/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Virtual_mcnp_ss_Guide
*
* %I
* Written by: Esben klinkby and Peter Willendrup
* Date: Marts 2012
* Origin: Risoe-DTU
*
* Neutron guide initiated using Virtual_mcnp_ss_input.comp, and replacing Virtual_mcnp_ss_output.comp  - see examples//Test_SSR_SSW_Guide.instr
*
* %D
* Based on Kristian Nielsens Guide.comp
* Models a rectangular guide tube centered on the Z axis. The entrance lies
* in the X-Y plane.
* The component must be initiated after Virtual_mcnp_ss_input.comp,
* and replaces Virtual_mcnp_ss_output.comp  - see examples/Test_SSR_SSW_Guide.instr.
* The basic idea is, that rather than discarding unreflected (i.e. absorbed)
* neutrons at the guide mirrors, these neutron states are stored on disk.
* Thus, after the McStas simulation a MCNP simulation can be performed based
*  on the un-reflected neutrons - intended for shielding studies.
* (details: we don't deal with actual neutrons, so what is transferred between
* simulations suites is neutron state parameters: pos,mom,time,weight. The latter is
* whatever remains after reflection.
*
* For details on the geometry calculation see the description in the McStas
* reference manual.
* The reflectivity profile may either use an analytical mode (see Component
* Manual) or a 2-columns reflectivity free text file with format
* [q(Angs-1) R(0-1)].
*
* Example: Virtual_mcnp_ss_Guide(w1=0.1, h1=0.1, w2=0.1, h2=0.1, l=2.0,
*           R0=0.99, Qc=0.021, alpha=6.07, m=2, W=0.003
*
* %VALIDATION
* Upcomming in 2012 based on ESS shielding
* Validated by: D. Ene & E. Klinkby
*
* %BUGS
* This component does not work with gravitation on. Use component Guide_gravity then
* (doesn't work with SSR/SSW unfortunately)
*
* %P
* INPUT PARAMETERS:
*
* w1: [m]         Width at the guide entry
* h1: [m]         Height at the guide entry
* w2: [m]         Width at the guide exit
* h2: [m]         Height at the guide exit
* l: [m]          length of guide
* R0: [1]         Low-angle reflectivity
* Qc: [AA-1]      Critical scattering vector
* alpha: [AA]     Slope of reflectivity
* m: [1]          m-value of material. Zero means completely absorbing.
* W: [AA-1]       Width of supermirror cut-off
* reflect: [str]  Reflectivity file name. Format <q(Angs-1) R(0-1)>
*
* %D
* Example values: m=4 Qc=0.0219 W=1/300 alpha=6.49 R0=1
*
* %E
*******************************************************************************/

DEFINE COMPONENT Virtual_mcnp_ss_Guide

SETTING PARAMETERS (string reflect=0, w1, h1, w2, h2, l, R0=0.99, Qc=0.0219, alpha=6.07, m=2, W=0.003)

NOACC

SHARE
%{
%include "read_table-lib"
%}

DECLARE
%{
t_Table      pTable;
double       from_mcstas[8];
unsigned int ntrk;
unsigned int nhis;
%}

INITIALIZE
%{
  writeheader_(&ntrk,&nhis); //Important that ntrk & nhis does not exceed actuall rssa file content

  if (mcgravitation)
    fprintf(stderr,"WARNING: Virtual_mcnp_ss_Guide: %s: "
      "This component produces wrong results with gravitation !\n"
      "Use Guide_gravity (doesn't work with SSR/SSW).\n",
      NAME_CURRENT_COMP);

  if (reflect && strlen(reflect))
  {
    if (Table_Read(&pTable, reflect, 1) <= 0) /* read 1st block data from file into pTable */
      exit(fprintf(stderr,"Virtual_mcnp_ss_Guide: %s: can not read file %s\n", NAME_CURRENT_COMP, reflect));
  }
  else
  {
    if (W < 0 || R0 < 0 || Qc < 0 || m < 0)
    {
      fprintf(stderr,"Virtual_mcnp_ss_Guide: %s: W R0 Qc must be >0.\n", NAME_CURRENT_COMP);
      exit(-1);
    }
    if (m < 1 && m != 0)
      fprintf(stderr,"WARNING: Virtual_mcnp_ss_Guide: %s: m < 1 behaves as if m=1.\n", NAME_CURRENT_COMP);
  }
%}

TRACE
%{
  double t1,t2;                                 /* Intersection times. */
  double av,ah,bv,bh,cv1,cv2,ch1,ch2,d;         /* Intermediate values */
  double weight;                                /* Internal probability weight */
  double vdotn_v1,vdotn_v2,vdotn_h1,vdotn_h2;   /* Dot products. */
  int i;                                        /* Which mirror hit? */
  double q;                                     /* Q [1/AA] of reflection */
  double nlen2;                                 /* Vector lengths squared */

  /* ToDo: These could be precalculated. */
  double ww = .5*(w2 - w1), hh = .5*(h2 - h1);
  double whalf = .5*w1, hhalf = .5*h1;

  double vx_org=0.,vy_org=0.,vz_org=0.;



  /* Propagate neutron to guide entrance. */
  PROP_Z0;
  /* Scatter here to ensure that fully transmitted neutrons will not be
     absorbed in a GROUP construction, e.g. all neutrons - even the
     later absorbed ones are scattered at the guide entry. */
  SCATTER;
  if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf)
    ABSORB;
  for(;;)
  {
    vx_org=vx;
    vy_org=vy;
    vz_org=vz;
    /* Compute the dot products of v and n for the four mirrors. */
    av = l*vx; bv = ww*vz;
    ah = l*vy; bh = hh*vz;
    vdotn_v1 = bv + av;         /* Left vertical */
    vdotn_v2 = bv - av;         /* Right vertical */
    vdotn_h1 = bh + ah;         /* Lower horizontal */
    vdotn_h2 = bh - ah;         /* Upper horizontal */
    /* Compute the dot products of (O - r) and n as c1+c2 and c1-c2 */
    cv1 = -whalf*l - z*ww; cv2 = x*l;
    ch1 = -hhalf*l - z*hh; ch2 = y*l;
    /* Compute intersection times. */
    t1 = (l - z)/vz;
    i = 0;
    if(vdotn_v1 < 0 && (t2 = (cv1 - cv2)/vdotn_v1) < t1)
    {
      t1 = t2;
      i = 1;
    }
    if(vdotn_v2 < 0 && (t2 = (cv1 + cv2)/vdotn_v2) < t1)
    {
      t1 = t2;
      i = 2;
    }
    if(vdotn_h1 < 0 && (t2 = (ch1 - ch2)/vdotn_h1) < t1)
    {
      t1 = t2;
      i = 3;
    }
    if(vdotn_h2 < 0 && (t2 = (ch1 + ch2)/vdotn_h2) < t1)
    {
      t1 = t2;
      i = 4;
    }
    if(i == 0)
      break;                    /* Neutron left guide. */
    PROP_DT(t1);
    switch(i)
    {
      case 1:                   /* Left vertical mirror */
        nlen2 = l*l + ww*ww;
        q = V2Q*(-2)*vdotn_v1/sqrt(nlen2);
        d = 2*vdotn_v1/nlen2;
        vx = vx - d*l;
        vz = vz - d*ww;
        break;
      case 2:                   /* Right vertical mirror */
        nlen2 = l*l + ww*ww;
        q = V2Q*(-2)*vdotn_v2/sqrt(nlen2);
        d = 2*vdotn_v2/nlen2;
        vx = vx + d*l;
        vz = vz - d*ww;
        break;
      case 3:                   /* Lower horizontal mirror */
        nlen2 = l*l + hh*hh;
        q = V2Q*(-2)*vdotn_h1/sqrt(nlen2);
        d = 2*vdotn_h1/nlen2;
        vy = vy - d*l;
        vz = vz - d*hh;
        break;
      case 4:                   /* Upper horizontal mirror */
        nlen2 = l*l + hh*hh;
        q = V2Q*(-2)*vdotn_h2/sqrt(nlen2);
        d = 2*vdotn_h2/nlen2;
        vy = vy + d*l;
        vz = vz - d*hh;
        break;
    }
    /* Now compute reflectivity. */
    weight = 1.0; /* Initial internal weight factor */
    if(m == 0)
      ABSORB;
    if (reflect && strlen(reflect))
      weight = Table_Value(pTable, q, 1);
    else if(q > Qc)
    {
      double arg = (q-m*Qc)/W;
      //      printf("   %e %e %e %e %e  \n ", arg, q, m, Qc, W );
      if(arg < 10)
        weight = .5*(1-tanh(arg))*(1-alpha*(q-Qc));
      else
        ABSORB;                               /* Cutoff ~ 1E-10 */
      weight *= R0;
    }
    else
    { /* q <= Qc */
      weight *= R0;
    }
    p *= weight;
    SCATTER;

    from_mcstas[0]=x;
    from_mcstas[1]=y;
    from_mcstas[2]=z;
    from_mcstas[3]=vx_org;
    from_mcstas[4]=vy_org;
    from_mcstas[5]=vz_org;
    from_mcstas[6]=p*(1.-weight); //Hmm... assuming that weigths has the same meening in MCNP and McStas.
    from_mcstas[7]=t;
    writeneutron_(&from_mcstas);
    //    double v2=vx_org*vx_org+vy_org*vy_org+vz_org*vz_org;
    //    printf("%e %e  %e   %e \n ", weight, 1-weight, from_mcstas[6], v2);
  }
%}

MCDISPLAY
%{

  multiline(5,
            -w1/2.0, -h1/2.0, 0.0,
             w1/2.0, -h1/2.0, 0.0,
             w1/2.0,  h1/2.0, 0.0,
            -w1/2.0,  h1/2.0, 0.0,
            -w1/2.0, -h1/2.0, 0.0);
  multiline(5,
            -w2/2.0, -h2/2.0, (double)l,
             w2/2.0, -h2/2.0, (double)l,
             w2/2.0,  h2/2.0, (double)l,
            -w2/2.0,  h2/2.0, (double)l,
            -w2/2.0, -h2/2.0, (double)l);
  line(-w1/2.0, -h1/2.0, 0, -w2/2.0, -h2/2.0, (double)l);
  line( w1/2.0, -h1/2.0, 0,  w2/2.0, -h2/2.0, (double)l);
  line( w1/2.0,  h1/2.0, 0,  w2/2.0,  h2/2.0, (double)l);
  line(-w1/2.0,  h1/2.0, 0, -w2/2.0,  h2/2.0, (double)l);
%}

END