/*******************************************************************************
*
* 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: <a href="http://www.ill.fr">ILL (France)</a>.
* 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
* <a href="../optics/Bender.comp.html">Bender</a>
*
* %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