/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright 1997-2002, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: V_selector
*
* %I
* Written by:  Kim Lefmann
* Date: Nov 25, 1998
* Origin: Risoe
*
* Velocity selector.
*
* %D
* Velocity selector consisting of rotating Soller-like blades
* defining a helically twisted passage.
* Geometry defined by two identical, centered apertures at 12 o'clock
* position, Origo is at the centre of the selector (input is at -zdepth/2).
* Transmission is analytical assuming a continuous source.
*
* Example: V_selector(xwidth=0.03, yheight=0.05, zdepth=0.30, radius=0.12, alpha=48.298, length=0.25, d=0.0004, nu=20000, nslit=72)
* These are values for the D11@ILL Dornier 'Dolores' Velocity Selector (NVS 023)
*
* %VALIDATION
* Jun 2005: extensive external test, no problems found
* Validated by: K. Lieutenant
*
* %P
* INPUT PARAMETERS:
*
* xwidth: [m]   Width of entry aperture
* yheight: [m]  Height of entry aperture
* zdepth: [m]   Distance between apertures, for housing containing the rotor
* radius: [m]   Height from aperture centre to rotation axis
* alpha: [deg]  Twist angle along the cylinder
* length: [m]   Length of cylinder/rotor (less than zdepth)
* d: [m]        Thickness of blades
* nu: [Hz]      Cylinder rotation speed, counter-clockwise, which is ideally 3956*alpha*DEG2RAD/2/PI/lambda/length
* nslit: [1]    Number of Soller blades
*
* %E
*******************************************************************************/

DEFINE COMPONENT V_selector



SETTING PARAMETERS (xwidth=0.03, yheight=0.05, zdepth=0.30, radius=0.12, alpha=48.298, length=0.25, d=0.0004, nu=300, nslit=72)


/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */

DECLARE
%{
  double omega;
  double alpha_rad;
%}

INITIALIZE
%{
  omega = nu * 2 * PI;
  alpha_rad = alpha * DEG2RAD;
  if (zdepth < length)
    zdepth = length;
%}

TRACE
%{
  double dt0, dt1;
  double r_i, r_f, r_mean;
  double theta_i, theta_f;
  double A;
  double d_s_alpha;

  if (vz == 0)
    ABSORB;
  dt1 = (-zdepth / 2.0 - z) / vz;
  PROP_DT (dt1); /* Propagate to the entry aperture */
  if (x < (-xwidth / 2.0) || x > (xwidth / 2.0) || y < (-yheight / 2.0) || y > (yheight / 2.0))
    ABSORB;

  dt0 = (zdepth - length) / (2.0 * vz); /* Propagate to the cylinder start */
  PROP_DT (dt0);
  r_i = sqrt (x * x + (y + radius) * (y + radius));
  theta_i = atan2 (x, y + radius);

  dt1 = length / vz; /* Propagate along the cylinder length */
  PROP_DT (dt1);
  r_f = sqrt (x * x + (y + radius) * (y + radius));
  theta_f = atan2 (x, y + radius);

  dt0 = (zdepth - length) / (2.0 * vz); /* Propagate to the exit aperture */
  PROP_DT (dt0);
  if (x < (-xwidth / 2.0) || x > (xwidth / 2.0) || y < (-yheight / 2.0) || y > (yheight / 2.0))
    ABSORB;

  /* Calculate analytical transmission assuming continuous source */

  r_mean = (r_i + r_f) / 2.0; /* Approximation using mean radius */
  d_s_alpha = theta_f - theta_i;
  A = nslit / (2 * PI) * (d / r_mean + fabs (alpha_rad + d_s_alpha - omega * length / vz));
  if (A >= 1)
    ABSORB;
  p *= (1 - A);
  SCATTER;
%}

MCDISPLAY
%{
  double r = radius + yheight;
  double x0 = -xwidth / 2.0;
  double x1 = xwidth / 2.0;
  double y0 = -yheight / 2.0;
  double y1 = yheight / 2.0;
  double z0 = -zdepth / 2.0;
  double z1 = -length / 2.0;
  double z2 = length / 2.0;
  double z3 = zdepth / 2.0;
  double a;
  double xw, yh;

  xw = xwidth / 2.0;
  yh = yheight / 2.0;
  /* Draw apertures */
  for (a = z0;;) {
    multiline (3, x0 - xw, (double)y1, a, (double)x0, (double)y1, a, (double)x0, y1 + yh, a);
    multiline (3, x1 + xw, (double)y1, a, (double)x1, (double)y1, a, (double)x1, y1 + yh, a);
    multiline (3, x0 - xw, (double)y0, a, (double)x0, (double)y0, a, (double)x0, y0 - yh, a);
    multiline (3, x1 + xw, (double)y0, a, (double)x1, (double)y0, a, (double)x1, y0 - yh, a);
    if (a == z3)
      break;
    else
      a = z3;
  }

  /* Draw cylinder. */
  circle ("xy", 0, -radius, z1, r);
  circle ("xy", 0, -radius, z2, r);
  line (0, -radius, z1, 0, -radius, z2);
  for (a = 0; a < 2 * PI; a += PI / 8) {
    multiline (4, 0.0, -radius, z1, r * cos (a), r * sin (a) - radius, z1, r * cos (a + DEG2RAD * alpha), r * sin (a + DEG2RAD * alpha) - radius, z2, 0.0,
               -radius, z2);
  }
%}

END