/*******************************************************************************
*
* McStas, the neutron ray-tracing package
*         Maintained by Kristian Nielsen and Kim Lefmann,
*         Copyright (C) 1997-2006 Risoe National Laboratory, Roskilde, Denmark
*
* %I
* Written by: E.M.Lauridsen, N.B.Christensen, A.B.Abrahamsen
* Date: 4.2.98
* Version: $Revision: 1.1 $
* Origin: McStas release
* Modified by: KL, KN 20.3.98 (rewrite)
* Modified by: KL, 28.09.01   (two lines)
* Modified by: KL, 22.05.03   (background)
*
* General powder sample with two scattering vectors.
*
* %D
* General powder sample with two scattering vectors. 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.
*
* %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]
* sigma_inc:Incoherent cross section per unit cell [fm^2]
* frac:     Fraction of incoherently scattered neutron rays [1]
* q_1:      Scattering vector of reflection, Line 1 [AA^-1]
* d_1:      d-spacing for sample, overrides 'q_1' [AA]
* j_1:      Multiplicity of reflection, Line 1 [1]
* F2_1:     Structure factor of reflection, Line 1 [fm^2]
* w_1:      Intrinsic line broadening of reflection, Line 1 [AA^-1]
* q_2:      Scattering vector of reflection, Line 2 [AA^-1]
* d_2:      d-spacing for sample, overrides 'q_2' [AA]
* j_2:      Multiplicity of reflection, Line 2 [1]
* F2_2:     Structure factor of reflection, Line 2 [fm^2]
* w_2:      Intrinsic line broadening of reflection, Line 2 [AA^-1]
* 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 Powder2
DEFINITION PARAMETERS ()
SETTING PARAMETERS (d_phi=0, radius=0.01, xwidth=0, zthick=0, yheight=0.05,
  pack=1, Vc=85, sigma_a=0.4, sigma_inc=0, frac=0, h=0,
  j_1, q_1, d_1=0, F2_1, w_2=0,
  j_2, q_2, d_2=0, F2_2, w_1=0, DW)
OUTPUT PARAMETERS (my_s_v2, my_a_v, q_v)
STATE PARAMETERS (x,y,z,vx,vy,vz,t,s1,s2,p)
DECLARE
%{
  double my_s_v2_1, my_s_v2_2, my_a_v, my_inc, q_v1, q_v2, w_v1, w_v2, v, solid_angle;
  char   isrect=0;
%}
INITIALIZE
%{
  if (h) yheight=h;

  if (!radius || !yheight) {
    if (!xwidth || !yheight || !zthick) exit(fprintf(stderr,"Powder2: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
    else isrect=1; }
  /* 100: barns -> fm^2 */
  my_a_v = 100*pack*sigma_a/Vc*2200;           /* Is not yet divided by v */
  my_inc = 100*pack*sigma_inc/Vc;
  my_s_v2_1 = 100*4*PI*PI*PI*pack*DW/(Vc*Vc*V2K*V2K)*(j_1*F2_1/q_1);
  my_s_v2_2 = 100*4*PI*PI*PI*pack*DW/(Vc*Vc*V2K*V2K)*(j_2*F2_2/q_2);
                                      /* Is not yet divided by v^2 */
  /* Squires [3.103] */
  if (d_1) q_1=2*PI/d_1;
  if (d_2) q_2=2*PI/d_2;
  q_v1 = q_1*K2V;
  q_v2 = q_2*K2V;
  w_v1 = w_1*K2V;
  w_v2 = w_2*K2V;
%}
TRACE
%{
  double t0, t1, v, v1,l_full, l, l_1, dt, d_phi0, theta, my_s, my_s_n;
  double aim_x, aim_y, aim_z, axis_x, axis_y, axis_z;
  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)
      t0 = 0;
    /* 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 */

    my_s = (my_s_v2_1+my_s_v2_2)/(v*v)+my_inc;
        /* Total attenuation from scattering */

    if (!frac || rand01() >= frac)
     {   /* Make coherent scattering event */

      /* choose line */
      if (rand01() > 0.5)  /* Select between the two powder lines */
       {
        arg = (q_v1+w_v1*randnorm())/(2.0*v);
        my_s_n = my_s_v2_1/(v*v);
       }
      else
       {
        arg = (q_v2+w_v2*randnorm())/(2.0*v);
        my_s_n = my_s_v2_2/(v*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;

      p *= 2*l_full*my_s_n*exp(-(my_a_v/v+my_s)*(l+l_1))/(1-frac);
/* printf("Powder p: %g \n",p); */
      }  /* Coherent scattering event */
    else
     {  /* Make incoherent scattering event */
      v = sqrt(vx*vx+vy*vy+vz*vz);
      if(d_phi) {
        randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
            0, 0, 1,
            2*PI, d_phi*DEG2RAD, ROT_A_CURRENT_COMP);
      } else {
        randvec_target_circle(&vx, &vy, &vz,
          &solid_angle, 0, 0, 1, 0);
      }
      p *= solid_angle/4*PI;
      v1 = sqrt(vx*vx+vy*vy+vz*vz);
      vx *= v/v1;
      vy *= v/v1;
      vz *= v/v1;
      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;

      p *= l_full*my_inc*exp(-(my_a_v/v+my_s)*(l+l_1))/(frac);
      p *= solid_angle/(4*PI);
     }  /* Incoherent scattering event */
  } /* else transmit */
%}

MCDISPLAY
%{
  magnify("xyz");
  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);
%}
END
