/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Guide_wavy
*
* %I
* Written by: Kim Lefmann
* Origin: Risoe
*
* Neutron guide with gaussian waviness.
*
* %D
* Models a rectangular guide tube centered on the Z axis. The entrance lies
* in the X-Y plane.
* For details on the geometry calculation see the description in the McStas
* reference manual.
*
* Example: m=2 Qc=0.0218 (nat. Ni) W=1/300 alpha=4.38 R0=0.995 (given by Daniel Clemens, PSI)
*
* %BUGS
* This component does not work with gravitation on. Use component Guide_gravity then.
*
* %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 for all mirrors
*
* alpha1: [AA]    Slope of reflectivity left
* m1: [1]         m-value of material, left
* W1: [AA^-1]      Width of supermirror cut-off left
* alpha2: [AA]    Slope of reflectivity right
* m2: [1]         m-value of material, right.
* W2: [AA^-1]      Width of supermirror cut-off right
* alpha3: [AA]    Slope of reflectivity top
* m3: [1]         m-value of material, top.
* W3: [AA^-1]      Width of supermirror cut-off top
* alpha4: [AA]    Slope of reflectivity bottom
* m4: [1]         m-value of material, bottom.
* W4: [AA^-1]      Width of supermirror cut-off bottom
*
* wavy_z: [deg]   Waviness in the z-(flight-)direction
* wavy_xy: [deg]  Waviness in the transverse direction
*
* %E
******************************************************************************/

DEFINE COMPONENT Guide_wavy

SETTING PARAMETERS (w1, h1, w2=0, h2=0, l,
  R0=0.995, Qc=0.0218, alpha=0, m=0, W=0,
  alpha1=4.38, m1=2, W1=0.003,
  alpha2=4.38, m2=2, W2=0.003,
  alpha3=4.38, m3=2, W3=0.003,
  alpha4=4.38, m4=2, W4=0.003,
  wavy_z=0, wavy_xy=0)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE %{
  %include "ref-lib"
%}
DECLARE
%{
  double whalf;
  double hhalf;
  double lwhalf;
  double lhhalf;
  double norm_nv;
  double norm_nh;
  double f_h;
  double f_v;
  double eta_z;
  double eta_xy;
%}

INITIALIZE
%{
  if (m) {
    m1 = m2 = m3 = m4 = m;
  }
  if (alpha) {
    alpha1 = alpha2 = alpha3 = alpha4 = alpha;
  }
  if (W) {
    W1 = W2 = W3 = W4 = W;
  }

  if (!w2)
    w2 = w1;
  if (!h2)
    h2 = h1;
  f_h = 0.5 * (w2 - w1), f_v = 0.5 * (h2 - h1);
  whalf = 0.5 * w1, hhalf = 0.5 * h1;
  lwhalf = l * whalf, lhhalf = l * hhalf;
  norm_nv = sqrt (l * l + f_v * f_v);
  norm_nh = sqrt (l * l + f_h * f_h);
  eta_z = wavy_z / (sqrt (8 * log (2))); /* Convert from FWHM to Gaussian sigma */
  eta_xy = wavy_xy / (sqrt (8 * log (2)));

  if (mcgravitation)
    fprintf (stderr,
             "WARNING: Guide_wavy: %s: "
             "This component produces wrong results with gravitation !\n"
             "Use Guide_gravity.\n",
             NAME_CURRENT_COMP);
%}

TRACE
%{
  double t_min, t_tmp;                           /* 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 vdotn;
  double d_xy, d_z; /* Random angles */
  double m0, alpha0, w;
  int i;                        /* Which mirror hit? */
  double q;                     /* Q [1/AA] of reflection */
  double dvx, dvy, dvz;         /* Velocity change */
  double vlen2, nlen2, norm_n2; /* Vector lengths squared */
  double nperp, pz, nxy;        /* for dot products */
  double R;                     /* Reflectivity */

  /* 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 (;;) {
    /* Compute the dot products of v and n for the four mirrors. */
    av = -l * vx;
    bv = f_h * vz;
    ah = -l * vy;
    bh = f_v * 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 * f_h;
    cv2 = x * l;
    ch1 = -hhalf * l - z * f_v;
    ch2 = y * l;
    /* Compute intersection times. */
    t_min = (l - z) / vz;
    /*  printf(" (x,y,z)=(%g %g %g) (vx,vy,vz)=(%g %g %g) Exit time : %g \n",
         x,y,z,vx,vy,vz,t_min); */
    i = 0;
    if (vdotn_v1 < 0 && (t_tmp = (cv1 + cv2) / vdotn_v1) < t_min) {
      t_min = t_tmp;
      i = 1;
      /*      printf("Left vertical: t=%g \n",t_min); */
    }
    if (vdotn_v2 < 0 && (t_tmp = (cv1 - cv2) / vdotn_v2) < t_min) {
      t_min = t_tmp;
      i = 2;
      /*      printf("Right vertical: t=%g \n",t_min); */
    }
    if (vdotn_h1 < 0 && (t_tmp = (ch1 + ch2) / vdotn_h1) < t_min) {
      t_min = t_tmp;
      i = 3;
      /*      printf("Lower horizontal: t=%g \n",t_min); */
    }
    if (vdotn_h2 < 0 && (t_tmp = (ch1 - ch2) / vdotn_h2) < t_min) {
      t_min = t_tmp;
      i = 4;
      /*      printf("Upper horizontal: t=%g \n",t_min); */
    }
    if (i == 0)
      break; /* Neutron left guide. */
    PROP_DT (t_min);

    /* ******* Recalculate dot products ********* */

    d_z = DEG2RAD * eta_z * randnorm ();
    d_xy = DEG2RAD * eta_xy * randnorm ();

    /* Now the normal vector is rotated. To 1st order in waviness and 2nd order
      in f=(w2-w1)/l and (f times waviness) the rotation matrix for
      the left vertical mirror is:

      {  1     d_xy      d_z      }
      { -d_xy  1         d_xy f_h }   (for the right vertical mirror the )
      {  d_z  -d_xy f_h  1        }   (terms d_xy f_h changes sign )

    the left vertical normal vector is  { -l,  0,  f_h }

    giving the rotated normal vector  { -l + f_h d_z,  l d_xy,  f_h - l d_z }

    for the right vertical mirror the normal vector is { l,  0,  f_h }

    and the rotated right normal vector is { l + f_h d_z,  -l d_xy,  f_h + l d_z }


    The top horizontal mirror must be (something like)

      {  1        -d_xy  d_xy f_h }
      {  d_xy      1     d_z      }   (for the bottom horizontal mirror the )
      { -d_xy f_h  d_z   1        }   (terms d_xy f_h changes sign )

    The top horizontal normal vector is { 0,  -l,  f_v}

    giving the rotated normal vector  { l d_xy,  -l + f_v d_z,  f_v - l d_z }

    for the bottom mirror the normal vector is { 0,  l,  f_h }

    and the rotated bottom normal vector is { -l d_xy,  l + f_v d_z,  f_v + l d_z }

    */

    switch (i) {
    case 1: /* Left vertical mirror */
      m0 = m1;
      w = W1;
      alpha0 = alpha1;
      norm_n2 = (-l + d_z * f_h) * (-l + d_z * f_h) + (d_xy * l) * (d_xy * l) + (f_h - d_z * l) * (f_h - d_z * l); /* Square of length of n vector */
      vdotn = (vx * (-l + f_h * d_z) + vy * (l * d_xy) + vz * (f_h - l * d_z));
      q = 2 * V2Q * fabs (vdotn) / sqrt (norm_n2);
      dvx = -2 * (-l + f_h * d_z) * vdotn / norm_n2;
      dvy = -2 * (l * d_xy) * vdotn / norm_n2;
      dvz = -2 * (f_h - l * d_z) * vdotn / norm_n2;
      break;
    case 2: /* Right vertical mirror */
      m0 = m2;
      w = W2;
      alpha0 = alpha2;
      norm_n2 = (l + d_z * f_h) * (l + d_z * f_h) + (-d_xy * l) * (-d_xy * l) + (f_h + d_z * l) * (f_h + d_z * l); /* Square of length of n vector */
      vdotn = (vx * (l + f_h * d_z) + vy * (-l * d_xy) + vz * (f_h + l * d_z));
      q = 2 * V2Q * fabs (vdotn) / sqrt (norm_n2);
      dvx = -2 * (l + f_h * d_z) * vdotn / norm_n2;
      dvy = -2 * (-l * d_xy) * vdotn / norm_n2;
      dvz = -2 * (f_h + l * d_z) * vdotn / norm_n2;
      break;
    case 3: /* Lower horizontal mirror */
      m0 = m3;
      w = W3;
      alpha0 = alpha3;
      norm_n2 = (d_xy * l) * (d_xy * l) + (-l + d_z * f_v) * (-l + d_z * f_v) + (f_v - d_z * l) * (f_v - d_z * l); /* Square of length of n vector */
      vdotn = (vx * (l * d_xy) + vy * (-l + f_v * d_z) + vz * (f_v - l * d_z));
      q = 2 * V2Q * fabs (vdotn) / sqrt (norm_n2);
      dvx = -2 * (l * d_xy) * vdotn / norm_n2;
      dvy = -2 * (-l + f_v * d_z) * vdotn / norm_n2;
      dvz = -2 * (f_v - l * d_z) * vdotn / norm_n2;
      break;
    case 4: /* Upper horizontal mirror */
      m0 = m4;
      w = W4;
      alpha0 = alpha4;
      norm_n2 = (-d_xy * l) * (-d_xy * l) + (l + d_z * f_v) * (l + d_z * f_v) + (f_v + d_z * l) * (f_v + d_z * l); /* Square of length of n vector */
      vdotn = (vx * (-l * d_xy) + vy * (l + f_v * d_z) + vz * (f_v + l * d_z));
      q = 2 * V2Q * fabs (vdotn) / sqrt (norm_n2);
      dvx = -2 * (-l * d_xy) * vdotn / norm_n2;
      dvy = -2 * (l + f_v * d_z) * vdotn / norm_n2;
      dvz = -2 * (f_v + l * d_z) * vdotn / norm_n2;
      break;
    default:
      printf ("Fatal error: No guide wall hit");
      exit (1);
    }
    /* Now compute reflectivity. */
    {
      double par[] = { R0, Qc, alpha0, m0, w };
      StdReflecFunc (q, par, &R);
      if (R > 0)
        p *= R;
      else
        ABSORB;
    }
    vx += dvx;
    vy += dvy;
    vz += dvz;
    SCATTER;
  }
%}

MCDISPLAY
%{
  polygon (4, -w1 / 2.0, h1 / 2.0, 0.0, w1 / 2.0, h1 / 2.0, 0.0, w2 / 2.0, h2 / 2.0, (double)l, -w2 / 2.0, h2 / 2.0, (double)l);
  polygon (4, -w1 / 2.0, -h1 / 2.0, 0.0, w1 / 2.0, -h1 / 2.0, 0.0, w2 / 2.0, -h2 / 2.0, (double)l, -w2 / 2.0, -h2 / 2.0, (double)l);
  polygon (4, -w1 / 2.0, h1 / 2.0, 0.0, -w1 / 2.0, -h1 / 2.0, 0.0, -w2 / 2.0, -h2 / 2.0, (double)l, -w2 / 2.0, h2 / 2.0, (double)l);
  polygon (4, w1 / 2.0, h1 / 2.0, 0.0, w1 / 2.0, -h1 / 2.0, 0.0, w2 / 2.0, -h2 / 2.0, (double)l, w2 / 2.0, h2 / 2.0, (double)l);
%}

END