/*******************************************************************************
*
* McStas, the neutron ray-tracing package
*         Maintained by Kristian Nielsen and Kim Lefmann,
*         Copyright 1997-2000 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: KL, KN
* Date: 15.4.98
* Version: $Revision: 1.1.1.1 $
* Origin: McStas release
*
* Vanadium sample.
*
* %D
* A Double-cylinder shaped incoherent scatterer (a V-sample)
* No multiple scattering. Absorbtion included.
*
* %P
* INPUT PARAMETERS:
*       
* radius_i  : Inner radius of sample in (x,z) plane (m)         
* radius_o  : Outer radius of sample in (x,z) plane (m)         
* h         : Height of sample y direction (m)
* pack      : Packing factor (1)
* focus_r   : Radius of sphere containing target. (m)
* target_x  :
* target_y  : position of target to focus at (m)
* target_z  :
*
* Variables calculated in the component
*
* V_my_s      : Attenuation factor due to scattering (m^-1)
* V_my_a      : Attenuation factor due to absorbtion (m^-1)
*
* %L
* <A HREF="http://neutron.risoe.dk/mcstas/components/tests/v_sample/">Test
* results</A> (not up-to-date).
*
* %E
*******************************************************************************/

DEFINE COMPONENT V_sample
DEFINITION PARAMETERS ()
SETTING PARAMETERS (radius_i,radius_o,h,pack,focus_r,
		    target_x, target_y, target_z)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
/* ToDo: Should be component local names. */
#define V_sigma_a 5.08     /* Absorption cross section per atom (barns) */
#define V_sigma_i 4.935    /* Incoherent scattering cross section per atom (barns) */
#define V_rho (2*pack/(3.024*3.024*3.024)) /* Density of atoms (AA-3) */
#define V_my_s (V_rho * 100 * V_sigma_i)
#define V_my_a_v (V_rho * 100 * V_sigma_a * 2200)
%}
INITIALIZE
%{
%}                                                                             
TRACE
%{
  double t0, t3;                /* Entry/exit time for outer cylinder */
  double t1, t2;                /* Entry/exit time for inner cylinder */
  double v;                     /* Neutron velocity */
  double dt0, dt1, dt2, dt;     /* Flight times through sample */
  double l_full;                /* Flight path length for non-scattered neutron */
  double l_i, l_o;              /* Flight path lenght in/out for scattered neutron */
  double my_a;                  /* Velocity-dependent attenuation factor */
  double solid_angle;           /* Solid angle of target as seen from scattering point */
  double aim_x, aim_y, aim_z;   /* Position of target relative to scattering point */

  if(cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius_o, h))
  {
    if(t0 < 0)
      ABSORB;
    /* Neutron enters at t=t0. */
    if(!cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius_i, h))
      t1 = t2 = t3;

    dt0 = t1-t0;                /* Time in sample, ingoing */
    dt1 = t2-t1;                /* Time in hole */
    dt2 = t3-t2;                /* Time in sample, outgoing */
    v = sqrt(vx*vx + vy*vy + vz*vz);
    l_full = v * (dt0 + dt2);   /* Length of full path through sample */
    dt = rand01()*(dt0+dt2);    /* Time of scattering (relative to t0) */
    l_i = v*dt;                 /* Penetration in sample */
    if (dt > dt0)
      dt += dt1;

    PROP_DT(dt+t0);             /* Point of scattering */
    
    aim_x = target_x-x;         /* Vector pointing at target (anal./det.) */
    aim_y = target_y-y;
    aim_z = target_z-z;
    randvec_target_sphere(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r);
    NORM(vx, vy, vz);
    vx *= v;
    vy *= v;
    vz *= v;

    if(!cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius_o, h))
    {
      /* ??? did not hit cylinder */
      printf("FATAL ERROR: Did not hit cylinder from inside.\n");
      exit(1);
    }
    dt = t3;
    if(cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius_i, h) &&
       t2 > 0)
      dt -= (t2-t1);            /* Subtract hollow part */
    l_o = v*dt;

    my_a = V_my_a_v/v;
    p *= l_full*V_my_s*exp(-(my_a+V_my_s)*(l_i+l_o));
    p /= 4*PI/solid_angle;     
  }
  else
    ABSORB;
%}

MCDISPLAY
%{
  magnify("xyz");
  circle("xz", 0,  h/2.0, 0, radius_i);
  circle("xz", 0,  h/2.0, 0, radius_o);
  circle("xz", 0, -h/2.0, 0, radius_i);
  circle("xz", 0, -h/2.0, 0, radius_o);
  line(-radius_i, -h/2.0, 0, -radius_i, +h/2.0, 0);
  line(+radius_i, -h/2.0, 0, +radius_i, +h/2.0, 0);
  line(0, -h/2.0, -radius_i, 0, +h/2.0, -radius_i);
  line(0, -h/2.0, +radius_i, 0, +h/2.0, +radius_i);
  line(-radius_o, -h/2.0, 0, -radius_o, +h/2.0, 0);
  line(+radius_o, -h/2.0, 0, +radius_o, +h/2.0, 0);
  line(0, -h/2.0, -radius_o, 0, +h/2.0, -radius_o);
  line(0, -h/2.0, +radius_o, 0, +h/2.0, +radius_o);
%}

END