/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright 1997-2002, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Tunneling_sample
*
* %I
* Written by: Kim Lefmann
* Date: 10.05.07
* Origin: Risoe
*
* A Double-cylinder shaped all-incoherent scatterer
* with elastic, quasielastic (Lorentzian), and tunneling (sharp)
* components.
*
* %D
* A Double-cylinder shaped all-incoherent scatterer
* with both elastic, quasielastic (Lorentzian), and tunneling (sharp)
* components. No multiple scattering. Absorbtion included.
* The shape of the sample may be a box with dimensions xwidth, yheight, zdepth.
* The area to scatter to is a disk of radius 'focus_r' situated at the target.
* This target area may also be rectangular if specified focus_xw and focus_yh
* or focus_aw and focus_ah, respectively in meters and degrees.
* The target itself is either situated according to given coordinates (x,y,z),
* or defined with the relative target_index of the component to focus
* to (next is +1).
* This target position will be set to its AT position. When targeting to
* centered components, such as spheres or cylinders, define an Arm component
* where to focus to.
*
* The outgoing polarization is calculated as for nuclear spin incoherence:
* P' = 1/3*P-2/3P = -1/3P
* As above multiple scattering is ignored .
*
* Example: Tunneling_sample(thickness=0.001,radius=0.01,yheight=0.02,focus_r=0.035,
* target_index=1)
*
* %P
* INPUT PARAMETERS:
* radius: [m] Outer radius of sample in (x,z) plane
* yheight: [m] vert. dimension of sample, as a height
* thickness: [m] Thickness of cylindrical sample in (x,z) plane
* focus_r: [m] Radius of disk containing target. Use 0 for full space
* target_index: [1] relative index of component to focus at, e.g. next is +1
* xwidth: horiz. dimension of sample, as a width [m]
* zdepth: depth of sample [m]
* focus_xw: horiz. dimension of a rectangular area [m]
* focus_yh: vert. dimension of a rectangular area [m]
* focus_aw: horiz. angular dimension of a rectangular area [deg]
* focus_ah: vert. angular dimension of a rectangular area [deg]
* sigma_abs:Absorbtion cross section pr. unit cell [barns]
* sigma_inc:Total incoherent scattering cross section pr. unit cell [barns]
* Vc: Unit cell volume [AA^3]
* p_interact: MC Probability for scattering the ray; otherwise transmit [1]
* f_QE: Fraction of quasielastic scattering [1]
* f_tun: Fraction of tunneling scattering (f_QE+f_tun < 1) [1]
* gamma: Lorentzian width of quasielastic broadening (HWHM) [meV]
* E_tun: Tunneling energy [meV]
* target_x: X-position of target to focus at [m]
* target_y: Y-position of target to focus at [m]
* target_z: Z-position of target to focus at [m]
*
* 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
* Test
* results (not up-to-date).
*
* %E
*******************************************************************************/
DEFINE COMPONENT Tunneling_sample
DEFINITION PARAMETERS ()
SETTING PARAMETERS (thickness=0, radius=0.01, focus_r = 0,
p_interact = 1, f_QE=0, f_tun=0, gamma=0, E_tun=0,
target_x = 0, target_y = 0, target_z = 0.235, focus_xw=0, focus_yh=0,
focus_aw=0, focus_ah=0, xwidth=0, yheight=0.05, zdepth=0, sigma_abs=5.08, sigma_inc=4.935, Vc=13.827, int target_index=0)
OUTPUT PARAMETERS()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
struct StructVarsV
{
double sigma_a; /* Absorption cross section per atom (barns) */
double sigma_i; /* Incoherent scattering cross section per atom (barns) */
double rho; /* Density of atoms (AA-3) */
double my_s;
double my_a_v;
char isrect; /* true when sample is a box */
double distance; /* when non zero, gives rect target distance */
double aw,ah; /* rectangular angular dimensions */
double xw,yh; /* rectangular metrical dimensions */
double tx,ty,tz; /* target coords */
};
%}
DECLARE
%{
struct StructVarsV VarsV;
double ftun;
double fQE;
%}
INITIALIZE
%{
if (!xwidth || !yheight || !zdepth) /* Cannot define a rectangle */
if (!radius || !yheight) /* Cannot define a cylinder either */
exit(fprintf(stderr,"V_sample: %s: sample has no volume (zero dimensions)\n", NAME_CURRENT_COMP));
else /* It is a cylinder */
VarsV.isrect=0;
else /* It is a rectangle */
VarsV.isrect=1;
VarsV.sigma_a=sigma_abs;
VarsV.sigma_i=sigma_inc;
VarsV.rho = (1/Vc);
VarsV.my_s=(VarsV.rho * 100 * VarsV.sigma_i);
VarsV.my_a_v=(VarsV.rho * 100 * VarsV.sigma_a);
/* now compute target coords if a component index is supplied */
VarsV.tx= VarsV.ty=VarsV.tz=0;
if (target_index)
{
Coords ToTarget;
ToTarget = coords_sub(POS_A_COMP_INDEX(INDEX_CURRENT_COMP+target_index),POS_A_CURRENT_COMP);
ToTarget = rot_apply(ROT_A_CURRENT_COMP, ToTarget);
coords_get(ToTarget, &VarsV.tx, &VarsV.ty, &VarsV.tz);
}
else
{ VarsV.tx = target_x; VarsV.ty = target_y; VarsV.tz = target_z; }
if (!(VarsV.tx || VarsV.ty || VarsV.tz))
printf("Tunneling_sample: %s: The target is not defined. Using direct beam (Z-axis).\n",
NAME_CURRENT_COMP);
VarsV.distance=sqrt(VarsV.tx*VarsV.tx+VarsV.ty*VarsV.ty+VarsV.tz*VarsV.tz);
/* different ways of setting rectangular area */
VarsV.aw = VarsV.ah = 0;
if (focus_xw) {
VarsV.xw = focus_xw;
}
if (focus_yh) {
VarsV.yh = focus_yh;
}
if (focus_aw) {
VarsV.aw = DEG2RAD*focus_aw;
}
if (focus_ah) {
VarsV.ah = DEG2RAD*focus_ah;
}
/* Check that probabilities are positive and do not exceed unity */
if (f_tun<0)
ftun=0;
else
ftun=f_tun;
if(f_QE<0)
fQE=0;
else
fQE=f_QE;
if ((ftun+fQE)>1) {
ftun=0;
printf("Tunneling_sample: Sum of inelastic probabilities > 1. Setting f_tun=0");
if (fQE>1) {
fQE=0;
printf("Tunneling_sample: Probability fQE > 1. Setting fQE=0.");
}
}
%}
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=0; /* Flight path lenght in/out for scattered neutron */
double my_a=0; /* Velocity-dependent attenuation factor */
double solid_angle=0; /* Solid angle of target as seen from scattering point */
double aim_x=0, aim_y=0, aim_z=1; /* Position of target relative to scattering point */
double v_i, v_f, E_i, E_f; /* initial and final energies and velocities */
double dE; /* Energy transfer */
double scatt_choice; /* Representing random choice of scattering type */
int intersect=0;
if (VarsV.isrect)
intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
else
intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight);
if(intersect)
{
if(t0 < 0) ABSORB; /* we already passed the sample; this is illegal */
/* Neutron enters at t=t0. */
if(VarsV.isrect)
t1 = t2 = t3;
else
if(!thickness || !cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius-thickness, yheight))
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 */
if (v) my_a = VarsV.my_a_v*(2200/v);
if (p_interact >= 1 || rand01() dt0)
dt += dt1; /* jump to 2nd side of cylinder */
PROP_DT(dt+t0); /* Point of scattering */
if ((VarsV.tx || VarsV.ty || VarsV.tz)) {
aim_x = VarsV.tx-x; /* Vector pointing at target (anal./det.) */
aim_y = VarsV.ty-y;
aim_z = VarsV.tz-z;
}
if(VarsV.aw && VarsV.ah) {
randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
aim_x, aim_y, aim_z, VarsV.aw, VarsV.ah, ROT_A_CURRENT_COMP);
} else if(VarsV.xw && VarsV.yh) {
randvec_target_rect(&vx, &vy, &vz, &solid_angle,
aim_x, aim_y, aim_z, VarsV.xw, VarsV.yh, ROT_A_CURRENT_COMP);
} else {
randvec_target_circle(&vx, &vy, &vz, &solid_angle, aim_x, aim_y, aim_z, focus_r);
}
NORM(vx, vy, vz);
scatt_choice = rand01(); /* chooses type of scattering */
v_i = v; /* Store initial velocity in case of inel. */
E_i = VS2E*v_i*v_i;
if (scatt_choice<(fQE+ftun)) /* Inelastic choices */
{
if (scatt_choice0)
dE = E_tun;
else
dE = -E_tun;
}
E_f = E_i + dE;
if (E_f <= 0)
ABSORB;
v_f = SE2V*sqrt(E_f);
v = v_f;
}
vx *= v;
vy *= v;
vz *= v;
if(!VarsV.isrect) {
if(!cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight))
{
/* ??? did not hit cylinder */
printf("FATAL ERROR: Did not hit cylinder from inside.\n");
exit(1);
}
dt = t3; /* outgoing point */
if(thickness && cylinder_intersect(&t1, &t2, x, y, z, vx, vy, vz, radius-thickness, yheight) &&
t2 > 0)
dt -= (t2-t1); /* Subtract hollow part */
}
else
{
if(!box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth))
{
/* ??? did not hit box */
printf("FATAL ERROR: Did not hit box from inside.\n");
exit(1);
}
dt = t3;
}
l_o = v*dt; /* trajectory after scattering point: absorption only */
p *= v/v_i*l_full*VarsV.my_s*exp(-my_a*(l_i+v_i/v*l_o)-VarsV.my_s*l_i);
/* We do not consider scattering from 2nd part (outgoing) */
p /= 4*PI/solid_angle;
p /= p_interact;
/* Polarisation part (1/3 NSF, 2/3 SF) */
sx *= -1.0/3.0;
sy *= -1.0/3.0;
sz *= -1.0/3.0;
SCATTER;
}
else /* Transmitting; always elastic */
{
p *= exp(-(my_a+VarsV.my_s)*l_full);
p /= (1-p_interact);
}
}
%}
MCDISPLAY
%{
if (!VarsV.isrect)
{
circle("xz", 0, yheight/2.0, 0, radius);
circle("xz", 0, -yheight/2.0, 0, radius);
line(-radius, -yheight/2.0, 0, -radius, +yheight/2.0, 0);
line(+radius, -yheight/2.0, 0, +radius, +yheight/2.0, 0);
line(0, -yheight/2.0, -radius, 0, +yheight/2.0, -radius);
line(0, -yheight/2.0, +radius, 0, +yheight/2.0, +radius);
if (thickness)
{
double radius_i=radius-thickness;
circle("xz", 0, yheight/2.0, 0, radius_i);
circle("xz", 0, -yheight/2.0, 0, radius_i);
line(-radius_i, -yheight/2.0, 0, -radius_i, +yheight/2.0, 0);
line(+radius_i, -yheight/2.0, 0, +radius_i, +yheight/2.0, 0);
line(0, -yheight/2.0, -radius_i, 0, +yheight/2.0, -radius_i);
line(0, -yheight/2.0, +radius_i, 0, +yheight/2.0, +radius_i);
}
}
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*zdepth;
double zmax = 0.5*zdepth;
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