/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2011, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Guide_curved
*
* %I
* Written by: Ross Stewart
* Date: November 18 2003
*
* Non-focusing curved neutron guide.
*
* %D
* Models a rectangular curved guide tube with entrance centered on the Z axis.
* The entrance lies in the X-Y plane.  Unlike Bender.comp, draws a true depiction
* of the guide, and trajectories.  Guide is not focusing
*
* Example: Guide(w1=0.1, h1=0.1, l=2.0, R0=0.99, Qc=0.021,
*                alpha=6.07, m=2, W=0.003, r=2700)
*
* %BUGS
* This component does not work with gravitation on. Use component Guide_gravity then.
* Systematic error on transmitted flux is found to be about 10%.
*
* %P
* INPUT PARAMETERS:
*
* w:       [m]    Width at the guide entry
* h:       [m]    Height at the guide entry
* 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
* r:       [m]    Radius of curvature of the guide
*
* %L
* <a href="../optics/Bender.comp.html">Bender</a>
*
* %E
*******************************************************************************/

DEFINE COMPONENT Guide_curved
DEFINITION PARAMETERS ()
SETTING PARAMETERS (w, h, l, R0=0.99, Qc=0.0219, alpha=6.07, m=2, W=0.003, r=2700)
OUTPUT PARAMETERS ()
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

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

TRACE
%{
  double t11, t12, t21, t22, theta, alphaAng, endtime, phi;
  double time, time1, time2, q, R;
  int ii, i_bounce;

  double whalf  = 0.5*w, hhalf = 0.5*h;   /* half width and height of guide */
  double z_off  = r*sin(l/r);       /* z-component of total guide length */
  double R1     = r - whalf;        /* radius of curvature of inside mirror */
  double R2     = r + whalf;        /* radius of curvature of outside mirror */
  double vel    = sqrt(vx*vx + vy*vy + vz*vz);  /* neutron velocity */
  double vel_xz = sqrt(vx*vx + vz*vz);      /* in plane velocity */
  double K      = V2K*vel;        /* neutron wavevector */
  double lambda = 2.0*PI/K;       /* neutron wavelength */

/* Propagate neutron to guide entrance. */

  PROP_Z0;
  if(x <= -whalf || x >= whalf || y <= -hhalf || y >= hhalf)
    ABSORB;

  for(;;)
  {
    /* Find itersection points of neutron with inside and outside guide walls */
    ii = cylinder_intersect(&t11, &t12 ,x - r, y, z, vx, vy, vz, R1, h);
    ii = cylinder_intersect(&t21, &t22 ,x - r, y, z, vx, vy, vz, R2, h);

    /* Choose appropriate reflection time */
    time1 = (t11 < 1e-7) ? t12 : t11;
    time2 = (t21 < 1e-7) ? t22 : t21;
    time  = (time1 < 1e-7 || time2 < time1) ? time2 : time1;

    /* Has neutron left the guide? */
    endtime = (z_off - z)/vz;
    if (time > endtime || time <= 1e-7) break;

    PROP_DT(time);

    /* Find reflection surface */
    R = (time == time1) ? R1 : R2;
    i_bounce = (fabs(y - hhalf) < 1e-7 || fabs(y + hhalf) < 1e-7) ? 2 : 1;
    switch(i_bounce) {
  case 1:           /* Inside or Outside wall */
            phi   = atan(vx/vz);        /* angle of neutron trajectory */
    alphaAng = asin(z/R);      /* angle of guide wall */
    theta = fabs(phi - alphaAng);    /* angle of reflection */
            vz    = vel_xz*cos(2.0*alphaAng - phi);
    vx    = vel_xz*sin(2.0*alphaAng - phi);
    break;
  case 2:       /* Top or Bottom wall */
    theta = fabs(atan(vy/vz));
    vy    = -vy;
    break;
    }
    /* Now compute reflectivity. */
    if (m == 0) ABSORB;

    q = 4.0*PI*sin(theta)/lambda;
    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;
  }
%}

MCDISPLAY
%{
  double y1 = h/2.0;
  double y2 = -y1;
  double x1, x2, z1, z2;
  double xplot1[100], xplot2[100], zplot1[100], zplot2[100];
  int n = 100;
  int j = 0;
  double R1 = (r - 0.5*w);    /* radius of inside arc */
  double R2 = (r + 0.5*w);    /* radius of outside arc */
  
  magnify("xy");

  for(;;) {
    if(j == n) break;
      z1 = ((double)j - 1.0)*(R1*l/r)/(double)(n - 1);
      z2 = ((double)j - 1.0)*(R2*l/r)/(double)(n - 1);
    x1 = r - sqrt(R1*R1 - z1*z1);
    x2 = r - sqrt(R2*R2 - z2*z2);
            xplot1[j] = x1;
    xplot2[j] = x2;
    zplot1[j] = z1;
    zplot2[j] = z2;
    j++;
  }
  /* Draw opening aperture */
  j=0;
  line(xplot1[j],y1,zplot1[j],xplot2[j],y1,zplot2[j]);
  line(xplot1[j],y1,zplot1[j],xplot1[j],y2,zplot1[j]);
  line(xplot2[j],y2,zplot2[j],xplot2[j],y1,zplot2[j]);
  line(xplot1[j],y2,zplot1[j],xplot2[j],y2,zplot2[j]);
  
  for(j=0;j<n-1;j++) {
    /* Arc elements in xz-plane(s) */
    line(xplot1[j],y1,zplot1[j],xplot1[j+1],y1,zplot1[j+1]);
    line(xplot2[j],y1,zplot2[j],xplot2[j+1],y1,zplot2[j+1]);
    line(xplot1[j],y2,zplot1[j],xplot1[j+1],y2,zplot1[j+1]);
    line(xplot2[j],y2,zplot2[j],xplot2[j+1],y2,zplot2[j+1]);
  }
  /* Draw closing aperture */
  j=n-1;
  line(xplot1[j],y1,zplot1[j],xplot2[j],y1,zplot2[j]);
  line(xplot1[j],y1,zplot1[j],xplot1[j],y2,zplot1[j]);
  line(xplot2[j],y2,zplot2[j],xplot2[j],y1,zplot2[j]);
  line(xplot1[j],y2,zplot1[j],xplot2[j],y2,zplot2[j]);
%}

END