/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, 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.6
*
* 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
* Old
* description, with explanation of focusing.
* %L
*
* Test results (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