/*******************************************************************************
*
* McStas, the neutron ray-tracing package
*         Maintained by Kristian Nielsen and Kim Lefmann,
*         Copyright 1997-2000 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: KL
* Date: November 22 2000
* Version: $Revision: 1.1.1.1 $
* Origin: McStas release, component Guide
* Modified by: 
*
* Neutron guide with gravity.
* NB! This is a temporary component to be abandones when information about
* gravity is included in the kernel.
*
* %D
* Models a rectangular guide tube centered on the Z axis. The entrance lies
* in the X-Y plane.
* A constant gravity is applied in the Y direction of the LOCAL coordinate
* system. No gravity is applied for the flight path before the guide. 
* For details on the geometry calculation see the description in the McStas
* reference manual.
*
* %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
*
* %D
* Example values: m=4 Qc=0.02 W=1/300 alpha=6.49 R0=1
*
* %E
*******************************************************************************/

DEFINE COMPONENT Gravity_guide
DEFINITION PARAMETERS ()
SETTING PARAMETERS (w1, h1, w2, h2, l, R0, Qc, alpha, m, W)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE 
%{
#define G 9.82
 double ww,hh,whalf,hhalf,lwhalf,lhhalf;
%}

INITIALIZE
%{
  ww = .5*(w2 - w1), hh = .5*(h2 - h1);
  whalf = .5*w1, hhalf = .5*h1;
  lwhalf = l*whalf, lhhalf = l*hhalf;
%}

TRACE
%{
  double t1,tv1,tv2,th1,th2;                                 /* Intersection times. */
  double av,ah,bv,bh,cv1,cv2,ch1,ch2,d;         /* Intermediate values */
  double vdotn_v1,vdotn_v2,vdotn_h1,vdotn_h2;   /* Dot products. */
  double D_1,b_1,D_2,b_2;        /* coeefficients for 2nd order equations */
  int i;                                        /* Which mirror hit? */
  double q;                                     /* Q [1/AA] of reflection */
  double vlen2,nlen2;                           /* Vector lengths squared */
 

  /* Propagate neutron to guide entrance. */
  /* NO GRAVITY HERE! */
  PROP_Z0;
  if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf)
    ABSORB;
  for(;;)
  {
    /* Compute the dot products of v and n for the two vertical mirrors. */
    av = l*vx; bv = ww*vz;
    vdotn_v1 = bv + av;         /* Left vertical */
    vdotn_v2 = bv - av;         /* Right vertical */
   /* Compute the dot products of (O - r) and n as c1+c2 and c1-c2 */
    cv1 = -whalf*l - z*ww; cv2 = x*l;
    /* Compute intersection times. WARNING, there may be division by zero */
    t1 = (l - z)/vz;
    tv1= (cv1 - cv2)/vdotn_v1;
    tv2 = (cv1 + cv2)/vdotn_v2;
    /* Now: solve the 2nd order equation for the "horizontal" surfaces*/
    b_1 = -vy+hh*vz/l;
    D_1 = b_1*b_1+2*G*(y-hhalf);
    if(D_1 < 0)
      th1= -1;
    else   /* Choose first intersection with upper surface */
      th1= (-b_1-sqrt(D_1))/G;

    b_2 = -vy-hh*vz/l;
    D_2 = b_2*b_2+2*G*(y+hhalf);
    if(D_2 < 0)
     {
      printf("Fatal error, neutron cannot intersect lower guide surface \n");
      exit(1);
     }
    else   /* Choose second intersection with lower surface */
      th2= (-b_2+sqrt(D_2))/G;

    i = 0;
    if(vdotn_v1 < 0 && tv1 < t1)
    {
      t1 = tv1;
      i = 1;
    }
    if(vdotn_v2 < 0 && tv2 < t1)
    {
      t1 = tv2;
      i = 2;
    }
    if(th1>0 && th1 < t1)
    {
      t1 = th1;
      i = 3;
    }
    if(th2>0 && th2 < t1)
    {
      t1 = th2;
      i = 4;
    }
    if(i == 0)
      break;                    /* Neutron left guide. */

    ah = l*vy; bh = hh*vz;
    vdotn_h1 = bh + ah;         /* Lower "horizontal" */
    vdotn_h2 = bh - ah;         /* Upper "horizontal" */
    ch1 = -hhalf*l - z*hh; ch2 = y*l;
      x+=vx*t1;
      z+=vz*t1;
      y+=vy*t1-G*t1*t1/2;
      vy-= G*t1;  /* TO BE REPLACED BY A MACRO: PROP_GRAV(dt) */
      t+=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. */
    if(m == 0)
      ABSORB;
    if(q > Qc)
    {
      double arg = (q-m*Qc)/W;
      if(arg < 10)
        p *= .5*(1-tanh(arg))*(1-alpha*(q-Qc));
      else
        ABSORB;                               /* Cutoff ~ 1E-10 */
    }
    p *= R0;
    SCATTER;

  }
  x+=vx*t1;
  z+=vz*t1;
  y+=vy*t1-G*t1*t1/2;
  vy-= G*t1;   /* TO BE REPLACED WITH A MACRO: PROP_GRAV(dt) */
  t+=t1;
%}

MCDISPLAY
%{
  double x;
  int i;

  magnify("xy");
  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