/*******************************************************************************
*
* McStas, the neutron ray-tracing package
* Maintained by Kristian Nielsen and Kim Lefmann,
* Copyright 1997-2000 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
* Old
* description, with explanation of focusing.
* %L
*
* Test results (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