/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2003, 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
* Origin: ILL (France) (http://www.ill.fr).
* Modified by: E. Farhi, uniformize parameter names (Jul 2008)
*
* 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.  Draws a true depiction
* of the guide, and trajectories.  Guide is not focusing.
*
* Example: Guide_curved(w1=0.1, h1=0.1, l=2.0, R0=0.99, Qc=0.021,
*                alpha=6.07, m=2, W=0.003, curvature=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:
*
* w1: [m]         Width at the guide entry
* h1: [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
* curvature: [m]  Radius of curvature of the guide
*
* %L
* Bender
*
* %E
*******************************************************************************/
DEFINE COMPONENT Guide_curved

SETTING PARAMETERS (w1, h1, l, R0=0.995, Qc=0.0218, alpha=4.38, m=2, W=0.003, curvature=2700)

/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "ref-lib"
%}

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*w1, hhalf = 0.5*h1;   /* half width and height of guide */
  double z_off  = curvature*sin(l/curvature);       /* z-component of total guide length */
  double R1     = curvature - whalf;        /* radius of curvature of inside mirror */
  double R2     = curvature + 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;
  SCATTER;
  for(;;)
  {
    double par[]={R0, Qc, alpha, m, W};
    /* Find itersection points of neutron with inside and outside guide walls */
    ii = cylinder_intersect(&t11, &t12 ,x - curvature, y, z, vx, vy, vz, R1, h1);
    ii = cylinder_intersect(&t21, &t22 ,x - curvature, y, z, vx, vy, vz, R2, h1);

    /* 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 || !R0) ABSORB;

    q = 4.0*PI*sin(theta)/lambda;
    StdReflecFunc(q, par, &R);
    if (R >= 0) p *= R; else ABSORB;
    SCATTER;
  }
%}

MCDISPLAY
%{
  double x1, x2, z1, z2;
  double xplot1[100], xplot2[100], zplot1[100], zplot2[100];
  int n = 100;
  int j = 1;
  double R1 = (curvature - 0.5*w1);    /* radius of inside arc */
  double R2 = (curvature + 0.5*w1);    /* radius of outside arc */

  

  for(j=0;j<n;j++) {
    z1 = ((double)j)*(R1*l/curvature)/(double)(n - 1);
    z2 = ((double)j)*(R2*l/curvature)/(double)(n - 1);
    x1 = curvature - sqrt(R1*R1 - z1*z1);
    x2 = curvature - sqrt(R2*R2 - z2*z2);
    xplot1[j] = x1;
    xplot2[j] = x2;
    zplot1[j] = z1;
    zplot2[j] = z2;
  }
  line(xplot1[0], 0.5*h1,zplot1[0],xplot2[0], 0.5*h1,zplot2[0]);
  line(xplot1[0], 0.5*h1,zplot1[0],xplot1[0],-0.5*h1,zplot1[0]);
  line(xplot2[0],-0.5*h1,zplot2[0],xplot2[0], 0.5*h1,zplot2[0]);
  line(xplot1[0],-0.5*h1,zplot1[0],xplot2[0],-0.5*h1,zplot2[0]);
  for(j=0;j<n-1;j++) {
    line(xplot1[j],  0.5*h1, zplot1[j], xplot1[j+1],  0.5*h1, zplot1[j+1]);
    line(xplot2[j],  0.5*h1, zplot2[j], xplot2[j+1],  0.5*h1, zplot2[j+1]);
    line(xplot1[j], -0.5*h1, zplot1[j], xplot1[j+1], -0.5*h1, zplot1[j+1]);
    line(xplot2[j], -0.5*h1, zplot2[j], xplot2[j+1], -0.5*h1, zplot2[j+1]);
  }
  line(xplot1[n-1], 0.5*h1,zplot1[n-1],xplot2[n-1], 0.5*h1,zplot2[n-1]);
  line(xplot1[n-1], 0.5*h1,zplot1[n-1],xplot1[n-1],-0.5*h1,zplot1[n-1]);
  line(xplot2[n-1],-0.5*h1,zplot2[n-1],xplot2[n-1], 0.5*h1,zplot2[n-1]);
  line(xplot1[n-1],-0.5*h1,zplot1[n-1],xplot2[n-1],-0.5*h1,zplot2[n-1]);
%}

END
