/*******************************************************************************
*
* 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
* Version: $Revision: 1.27 $
* Origin: Risoe
* Release: McStas 1.6
*
* 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 -l0/2).
* Transmission is analytical assuming a continuous source.
*
* Example: V_selector(width=0.03, height=0.05, l0=0.30, r0=0.12, phi=42.298,
*            l1=0.25, tb=0.0004, rot=20000, nb=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:
*
* width:     (m)   Width of entry aperture
* height:    (m)   Height of entry aperture
* l0:        (m)   Distance between apertures
* r0:        (m)   Height from aperture centre to rotation axis
* phi:       (deg) Twist angle along the cylinder
* l1:        (m)   Length of cylinder (less than l0)
* tb:        (m)   Thickness of blades
* rot:       (rpm) Cylinder rotation speed, counter-clockwise
* nb:        (1)   Number of Soller blades
*
* %E
*******************************************************************************/


DEFINE COMPONENT V_selector
DEFINITION PARAMETERS ()
SETTING PARAMETERS (width=0.03, height=0.05, l0=0.30, r0=0.12, phi=48.298, l1=0.25, tb=0.0004, rot=20000, nb=72)
OUTPUT PARAMETERS(RPM2OM, omega, phi_rad, dt0, dt1, r_i, r_f, r_mean, theta_i, theta_f, A, d_s_phi)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE
%{
  double RPM2OM, omega, phi_rad, dt0, dt1, r_i, r_f, r_mean, theta_i, theta_f, A, d_s_phi;
%}
INITIALIZE
%{
    RPM2OM = 2*PI/60.0;
    omega=rot*RPM2OM;
    phi_rad = phi*DEG2RAD;
%}
TRACE
%{
    if (vz == 0)
      ABSORB;
    dt1= (-l0/2.0 - z)/vz;
    PROP_DT(dt1); /* Propagate to the entry aperture */
    if (x<(-width/2.0) || x>(width/2.0) || y<(-height/2.0) || y>(height/2.0))
      ABSORB;

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

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

    dt0 = (l0-l1)/(2.0*vz); /* Propagate to the exit aperture */ 
    PROP_DT(dt0);
    if (x<(-width/2.0) || x>(width/2.0) || y<(-height/2.0) || y>(height/2.0))
      ABSORB;

    /* Calculate analytical transmission assuming continuous source */
    
    r_mean = (r_i + r_f)/2.0;          /* Approximation using mean radius */
    d_s_phi = theta_f-theta_i;
    A = nb/(2*PI)*( tb/r_mean + fabs(phi_rad+d_s_phi-omega*l1/vz) ); 
    if (A >= 1)
      ABSORB;
    p*= (1-A);
    SCATTER;
%}

MCDISPLAY
%{
  double r = r0 + height;
  double x0 = -width/2.0;
  double x1 =  width/2.0;
  double y0 = -height/2.0;
  double y1 =  height/2.0;
  double z0 = -l0/2.0;
  double z1 = -l1/2.0;
  double z2 =  l1/2.0;
  double z3 =  l0/2.0;
  double a;
  double xw, yh;
  
  magnify("xy");
  xw = width/2.0;
  yh = height/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, -r0, z1, r);
  circle("xy", 0, -r0, z2, r);
  line(0, -r0, z1, 0, -r0, z2);
  for(a = 0; a < 2*PI; a += PI/8)
  {
    multiline(4,
              0.0, -r0, z1,
              r*cos(a), r*sin(a) - r0, z1,
              r*cos(a + DEG2RAD*phi), r*sin(a + DEG2RAD*phi) - r0, z2,
              0.0, -r0, z2);
  }
%}

END