/*******************************************************************************
*
* McStas, neutron ray-tracing package
*         Copyright (C) 1997-2006, All rights reserved
*         Risoe National Laboratory, Roskilde, Denmark
*         Institut Laue Langevin, Grenoble, France
*
* Component: Powder1
*
* %I
* Written by: E.M.Lauridsen, N.B.Christensen, A.B.Abrahamsen
* Date: 4.2.98
* Version: $Revision: 1.29 $
* Origin: Risoe
* Release: McStas 1.9.1
*
* General powder sample with a single scattering vector.
*
* %D
* General powder sample with a single scattering vector. No multiple ,
* scattering, no incoherent scattering, no secondary extinction.
* The shape of the sample may be a cylinder of given radius or a box with
* dimensions xwidth, yheight, zthick.
* The efficient is highly improved when restricting the vertical scattering
* range on the Debye-Scherrer cone (with 'd_phi').
* You may use PowderN to use N scattering lines defined in a file.
*
* Example: Powder1(radius=0.015,yheight=0.05,q =1.8049,d_phi=0.07,pack=1,
*       j=6,DW=1,F2=56.8,Vc=85.0054,sigma_a=0.463)
*
* %P
*
* INPUT PARAMETERS
*
* d_phi:    Angle corresponding to the vertical angular range
*             to focus to, e.g. detector height. 0 for no focusing [deg,0-180]
* radius:  Radius of sample in (x,z) plane [m]
* yheight: Height of sample y direction [m]
* pack:    Packing factor [1]
* Vc:      Volume of unit cell [AA^3]
* sigma_a: Absorption cross section per unit cell at 2200 m/s [fm^2]
*
* q:       Scattering vector of reflection [AA^-1]
* d:       d-spacing for sample, overrides 'q' [AA]
* j:       Multiplicity of reflection [1]
* F2:      Structure factor of reflection [fm^2]
* DW:      Debye-Waller factor of reflection [1]
*
* Optional parameters:
* xwidth:  horiz. dimension of sample, as a width [m]
* zthick:  thickness of sample [m]
* h:       the same as yheight [m]
*
* Variables calculated in the component
*
* my_s: Attenuation factor due to scattering [m^-1]
* my_a: Attenuation factor due to absorbtion [m^-1]
*
* %L
* <A HREF="http://neutron.risoe.dk/mcstas/components/Powder1/">Old
* description</A>, with explanation of focusing.
* %L
* <A HREF="http://neutron.risoe.dk/mcstas/components/tests/powder/">
* Test results</A> (not up-to-date).
* %L
* See also: Powder1, Powder2 and PowderN
*
* %E
*******************************************************************************/

DEFINE COMPONENT Powder1
DEFINITION PARAMETERS ()
SETTING PARAMETERS (radius=0.01, yheight=0.05, q= 1.8049, d=0, d_phi= 0,
  pack= 1, j= 6, DW= 1, F2= 56.8, Vc= 85.0054, sigma_a= 0.463,
  xwidth=0, zthick=0, h=0)
OUTPUT PARAMETERS (my_s_v2, my_a_v, q_v, isrect)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)

DECLARE
%{
  double my_s_v2, my_a_v, q_v;
  char   isrect=0;
%}
INITIALIZE
%{
  if (h) yheight=h;
  if (!radius || !yheight) {
    if (!xwidth || !yheight || !zthick) exit(fprintf(stderr,"Powder1: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
    else isrect=1; }

  my_a_v = pack*sigma_a/Vc*2200*100;        /* Is not yet divided by v */
  my_s_v2 = 4*PI*PI*PI*pack*j*F2*DW/(Vc*Vc*V2K*V2K*q)*100;
  /* Is not yet divided by v^2. 100: convert from barns to fm^2 */
  /* Squires [3.103] */
  if (d) q=2*PI/d;
  q_v = q*K2V;
%}
TRACE
%{
  double t0, t1, v, l_full, l, l_1, dt, d_phi0, theta, my_s;
  double arg, tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z;
  char   intersect=0;

  if (isrect)
    intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zthick);
  else
    intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight);
  if(intersect)
  {
    if(t0 < 0)
      ABSORB;
    /* Neutron enters at t=t0. */
    v = sqrt(vx*vx + vy*vy + vz*vz);
    l_full = v * (t1 - t0);        /* Length of full path through sample */
    dt = rand01()*(t1 - t0);       /* Time of scattering */
    PROP_DT(dt+t0);                /* Point of scattering */
    l = v*dt;                      /* Penetration in sample */

    /* choose line theta */
    arg = q_v/(2.0*v);
    if(arg > 1)
      ABSORB;                   /* No bragg scattering possible*/
    theta = asin(arg);          /* Bragg scattering law */

/* Choose point on Debye-Scherrer cone */
      if (d_phi)
      { /* relate height of detector to the height on DS cone */
        arg = sin(d_phi*DEG2RAD/2)/sin(2*theta);
        if (arg < -1 || arg > 1) d_phi = 0;
        else d_phi = 2*asin(arg);
      }
      if (d_phi) {
        d_phi = fabs(d_phi);
        d_phi0= 2*rand01()*d_phi;
        if (d_phi0 > d_phi) arg = 1; else arg = 0;
        if (arg) {
          d_phi0=PI+(d_phi0-1.5*d_phi);
        } else {
          d_phi0=d_phi0-0.5*d_phi;
        }
        p *= d_phi/PI;
      }
      else
        d_phi0 = PI*randpm1();

    /* now find a nearly vertical rotation axis:
      *  (v along Z) x (X axis) -> nearly Y axis
      */
    vec_prod(tmp_vx,tmp_vy,tmp_vz, vx,vy,vz, 1,0,0);

    /* handle case where v and aim are parallel */
    if (!tmp_vx && !tmp_vy && !tmp_vz) { tmp_vx=tmp_vz=0; tmp_vy=1; }

    /* v_out = rotate 'v' by 2*theta around tmp_v: Bragg angle */
    rotate(vout_x,vout_y,vout_z, vx,vy,vz, 2*theta, tmp_vx,tmp_vy,tmp_vz);

    /* tmp_v = rotate v_out by d_phi0 around 'v' (Debye-Scherrer cone) */
    rotate(tmp_vx,tmp_vy,tmp_vz, vout_x,vout_y,vout_z, d_phi0, vx, vy, vz);
    vx = tmp_vx;
    vy = tmp_vy;
    vz = tmp_vz;

    arg=0;
    if (isrect && !box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zthick)) arg=1;
    else if(!isrect && !cylinder_intersect(&t0, &t1, x, y, z,
                          vx, vx, vx, radius, yheight)) arg=1;

    if (arg) {
      /* Strange error: did not hit cylinder */
      fprintf(stderr, "PowderN: FATAL ERROR: Did not hit sample from inside.\n");
      ABSORB;
    }
    l_1 = v*t1; /* go to exit */

    my_s = my_s_v2/(v*v);
    p *= l_full*my_s*exp(-(my_a_v/v+my_s)*(l+l_1));
    SCATTER;
  }
%}

MCDISPLAY
%{
  double h;
  h=yheight;
  magnify("xyz");
  if (!isrect) {
    circle("xz", 0,  h/2.0, 0, radius);
    circle("xz", 0, -h/2.0, 0, radius);
    line(-radius, -h/2.0, 0, -radius, +h/2.0, 0);
    line(+radius, -h/2.0, 0, +radius, +h/2.0, 0);
    line(0, -h/2.0, -radius, 0, +h/2.0, -radius);
    line(0, -h/2.0, +radius, 0, +h/2.0, +radius);
  } else {
    double xmin = -0.5*xwidth;
    double xmax =  0.5*xwidth;
    double ymin = -0.5*yheight;
    double ymax =  0.5*yheight;
    double zmin = -0.5*zthick;
    double zmax =  0.5*zthick;
    multiline(5, xmin, ymin, zmin,
                 xmax, ymin, zmin,
                 xmax, ymax, zmin,
                 xmin, ymax, zmin,
                 xmin, ymin, zmin);
    multiline(5, xmin, ymin, zmax,
                 xmax, ymin, zmax,
                 xmax, ymax, zmax,
                 xmin, ymax, zmax,
                 xmin, ymin, zmax);
    line(xmin, ymin, zmin, xmin, ymin, zmax);
    line(xmax, ymin, zmin, xmax, ymin, zmax);
    line(xmin, ymax, zmin, xmin, ymax, zmax);
    line(xmax, ymax, zmin, xmax, ymax, zmax);
  }
%}
END
