/*******************************************************************************
*
* McStas, neutron ray-tracing package
* Copyright (C) 1997-2008, All rights reserved
* Risoe National Laboratory, Roskilde, Denmark
* Institut Laue Langevin, Grenoble, France
*
* Component: Incoherent
*
* %I
* Written by: Kim Lefmann and Kristian Nielsen
* Date: 15.4.98
* Origin: Risoe
* Modified by: Aziz Daoud-aladine, ISIS, 2007: Added option to handle a spherical sample shape
* Modified by: Peter Christiansen, Risoe: Added outgoing polarization: P' = 1/3*P-2/3P = -1/3P NB! As above multiple scattering is ignored .
* Modified by: Reynald Arnerin, ILL, 2008: Added option to handle a complex geometry (OFF files)
*
* Incoherent sample (such as Vanadium) sample, with quasielastic component OR or global energy transfer.
*
* %D
* A Double-cylinder shaped incoherent scatterer (like Vanadium)
* with both elastic and quasielastic (Lorentzian) components.
* No multiple scattering (but approximation available). Absorption included.
* Sample focusing:
* 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.
* Sample shape:
* Sample shape may be a cylinder, a sphere, a box or any other shape
* box/plate: xwidth x yheight x zdepth (thickness=0)
* hollow box/plate:xwidth x yheight x zdepth and thickness>0
* cylinder: radius x yheight (thickness=0)
* hollow cylinder: radius x yheight and thickness>0
* sphere: radius (yheight=0 thickness=0)
* hollow sphere: radius and thickness>0 (yheight=0)
* any shape: geometry=OFF file
*
* The complex geometry option handles any closed non-convex polyhedra.
* It computes the intersection points of the neutron ray with the object
* transparently, so that it can be used like a regular sample object.
* It supports the PLY, OFF and NOFF file format but not COFF (colored faces).
* Such files may be generated from XYZ data using:
* qhull < coordinates.xyz Qx Qv Tv o > geomview.off
* or
* powercrust coordinates.xyz
* and viewed with geomview or java -jar jroff.jar (see below).
* The default size of the object depends of the OFF file data, but its
* bounding box may be resized using xwidth,yheight and zdepth.
*
* Example: Incoherent(radius=0.05,focus_r=0.035, pack=1, target_index=1)
* Incoherent(geometry="socket.off",focus_r=0.035, pack=1, target_index=1)
*
* %P
* INPUT PARAMETERS:
* radius: [m] Outer radius of sample in (x,z) plane
* xwidth: [m] Horiz. dimension of sample (bounding box if off file), as a width
* yheight: [m] Vert. dimension of sample (bounding box if off file), as a height. A sphere shape is used when 0 and radius is set
* zdepth: [m] Depth of sample (bounding box if off file)
* thickness: [m] Thickness of hollow sample
* 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
*
* Optional parameters:
* pack: [1] Packing factor
* p_interact: [1] MC Probability for scattering the ray; otherwise transmit
* focus_xw: [m] horiz. dimension of a rectangular area
* focus_yh: [m] vert. dimension of a rectangular area
* focus_aw: [deg] horiz. angular dimension of a rectangular area
* focus_ah: [deg] vert. angular dimension of a rectangular area
* sigma_abs: [barns] Absorption cross section pr. unit cell at 2200 m/s
* sigma_inc: [barns] Incoherent scattering cross section pr. unit cell
* Vc: [AA^3] Unit cell volume
* f_QE: [1] Fraction of quasielastic scattering (rest is elastic)
* gamma: [1] Lorentzian width of quasielastic broadening (HWHM)
* Etrans: [meV] Global energy-transfer, for use in inelastic settings
* deltaE: [meV] Width in energy around Etrans, for use in inelastic settings
* geometry: [str] Name of an Object File Format (OFF) or PLY file for complex geometry. The OFF/PLY file may be generated from XYZ coordinates using qhull/powercrust
* concentric: [1] Indicate that this component has a hollow geometry and may contain other components. It should then be duplicated after the inside part (only for box, cylinder, sphere)
*
* order: [1] Limit multiple scattering up to given order 0 means all (default), 1 means single, 2 means double, ...
* target_x: []
* target_y: [m] position of target to focus at
* target_z: []
*
* %L
* Cross sections for single elements
* %L
* Web Elements
* %L
* Test
* results (not up-to-date).
* %L
* The test/example instrument vanadium_example.instr.
* %L
* The test/example instrument QENS_test.instr.
* %L
* Geomview and Object File Format (OFF)
* %L
* Java version of Geomview (display only) jroff.jar
* %L
* qhull
* %L
* powercrust
*
* %E
*******************************************************************************/
DEFINE COMPONENT Incoherent
DEFINITION PARAMETERS ()
SETTING PARAMETERS (string geometry=0, radius=0, xwidth=0, yheight=0, zdepth=0, thickness=0,
target_x = 0, target_y = 0, target_z = 0, focus_r = 0,
focus_xw=0, focus_yh=0, focus_aw=0, focus_ah=0, int target_index=0,
pack = 1, p_interact=1, f_QE=0, gamma=0, Etrans=0,deltaE=0,
sigma_abs=5.08, sigma_inc=5.08, Vc=13.827, concentric=0, order=0)
OUTPUT PARAMETERS()
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
SHARE
%{
%include "read_table-lib"
%include "interoff-lib"
struct StructVarsInc
{
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;
int shape; /* 0 cylinder, 1 box, 2 sphere, 3 OFF file */
double aw,ah; /* rectangular angular dimensions */
double xw,yh; /* rectangular metrical dimensions */
double tx,ty,tz; /* target coords */
};
%}
DECLARE
%{
struct StructVarsInc VarsInc;
off_struct offdata;
%}
INITIALIZE
%{
VarsInc.shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere, 3:any-shape */
if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) {
#ifndef USE_OFF
fprintf(stderr,"Error: You are attempting to use an OFF geometry without -DUSE_OFF. You will need to recompile with that define set!\n");
exit(-1);
#else
if (off_init(geometry, xwidth, yheight, zdepth, 0, &offdata)) {
VarsInc.shape=3; thickness=0; concentric=0;
}
#endif
}
else if (xwidth && yheight && zdepth) VarsInc.shape=1; /* box */
else if (radius > 0 && yheight) VarsInc.shape=0; /* cylinder */
else if (radius > 0 && !yheight) VarsInc.shape=2; /* sphere */
if (VarsInc.shape < 0)
exit(fprintf(stderr,"Incoherent: %s: sample has invalid dimensions.\n"
"ERROR Please check parameter values (xwidth, yheight, zdepth, radius).\n", NAME_CURRENT_COMP));
if (thickness) {
if (radius && (radius < thickness || ( yheight && (yheight < 2*thickness)))) {
fprintf(stderr,"Incoherent: %s: hollow sample thickness is larger than its volume (sphere/cylinder).\n"
"WARNING Please check parameter values. Using bulk sample (thickness=0).\n", NAME_CURRENT_COMP);
thickness=0;
}
else if (!radius && (xwidth < 2*thickness || yheight < 2*thickness || zdepth < 2*thickness)) {
fprintf(stderr,"Incoherent: %s: hollow sample thickness is larger than its volume (box).\n"
"WARNING Please check parameter values. Using bulk sample (thickness=0).\n", NAME_CURRENT_COMP);
thickness=0;
}
}
if (concentric && thickness<=0) {
printf("Incoherent: %s:Can not use concentric mode\n"
"WARNING on non hollow shape. Ignoring.\n",
NAME_CURRENT_COMP);
concentric=0;
}
VarsInc.sigma_a= sigma_abs;
VarsInc.sigma_i= sigma_inc;
VarsInc.rho = (pack/Vc);
VarsInc.my_s = (VarsInc.rho * 100 * VarsInc.sigma_i);
VarsInc.my_a_v = (VarsInc.rho * 100 * VarsInc.sigma_a);
/* now compute target coords if a component index is supplied */
VarsInc.tx= VarsInc.ty=VarsInc.tz=0;
if (!target_index && !target_x && !target_y && !target_z) target_index=1;
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, &VarsInc.tx, &VarsInc.ty, &VarsInc.tz);
}
else
{ VarsInc.tx = target_x; VarsInc.ty = target_y; VarsInc.tz = target_z; }
if (!(VarsInc.tx || VarsInc.ty || VarsInc.tz)) {
MPI_MASTER(
printf("Incoherent: %s: The target is not defined. Using direct beam (Z-axis).\n",
NAME_CURRENT_COMP);
);
VarsInc.tz=1;
}
/* different ways of setting rectangular area */
VarsInc.aw = VarsInc.ah = 0;
if (focus_xw) { VarsInc.xw = focus_xw; }
if (focus_yh) { VarsInc.yh = focus_yh; }
if (focus_aw) { VarsInc.aw = DEG2RAD*focus_aw; }
if (focus_ah) { VarsInc.ah = DEG2RAD*focus_ah; }
MPI_MASTER(
printf("Incoherent: %s: Vc=%g [Angs] sigma_abs=%g [barn] sigma_inc=%g [barn]\n",
NAME_CURRENT_COMP, Vc, VarsInc.sigma_a, VarsInc.sigma_i);
);
%}
TRACE
%{
double t0, t3; /* Entry/exit time for outer surface */
double t1, t2; /* Entry/exit time for inner surface */
double dt0, dt1, dt2, dt; /* Flight times through sample */
double v=0; /* Neutron velocity */
double d_path; /* 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,my_t=0; /* Velocity-dependent attenuation factor and total Xsec */
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 */
int intersect=0;
int flag_concentric=0;
int flag=0;
double mc_trans, p_trans, mc_scatt, p_scatt, ws;
double p_mult=1;
#ifdef OPENACC
#ifdef USE_OFF
off_struct thread_offdata = offdata;
#endif
#else
#define thread_offdata offdata
#endif
do { /* Main interaction loop. Ends with intersect=0 */
/* Intersection neutron trajectory / sample (sample surface) */
if (VarsInc.shape == 0)
intersect = cylinder_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius, yheight);
else if (VarsInc.shape == 1)
intersect = box_intersect(&t0, &t3, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
else if (VarsInc.shape == 2)
intersect = sphere_intersect(&t0, &t3, x, y, z, vx, vy, vz, radius);
#ifdef USE_OFF
else if (VarsInc.shape == 3)
intersect = off_intersect(&t0, &t3, NULL, NULL, x, y, z, vx, vy, vz, 0, 0, 0, thread_offdata );
#endif
if (intersect) {
int flag_ishollow = 0;
if (thickness>0) {
if (VarsInc.shape==0 && cylinder_intersect(&t1,&t2, x,y,z,vx,vy,vz, radius-thickness,yheight-2*thickness))
flag_ishollow=1;
else if (VarsInc.shape==2 && sphere_intersect (&t1,&t2, x,y,z,vx,vy,vz, radius-thickness))
flag_ishollow=1;
else if (VarsInc.shape==1 && box_intersect(&t1,&t2, x,y,z,vx,vy,vz, xwidth-2*thickness, yheight-2*thickness, zdepth-2*thickness))
flag_ishollow = 1;
}
if (!flag_ishollow) t1 = t2 = t3; /* no empty space inside */
dt0 = t1-t0; /* Time in sample, ingoing */
dt1 = t2-t1; /* Time in hole */
dt2 = t3-t2; /* Time in sample, outgoing */
if (t0 > 0) { /* we are before the sample */
PROP_DT(t0); /* propagates neutron to the entry of the sample */
} else if (t1 > 0 && t1 > t0) { /* we are inside first part of the sample */
/* no propagation, stay inside */
} else if (t2 > 0 && t2 > t1) { /* we are in the hole */
PROP_DT(t2); /* propagate to inner surface of 2nd part of sample */
} else if (t3 > 0 && t3 > t2) { /* we are in the 2nd part of sample */
/* no propagation, stay inside */
}
dt0=t1-(t0 > 0 ? t0 : 0); /* Time in first part of hollow/cylinder/box */
dt1=t2-(t1 > 0 ? t1 : 0); /* Time in hole */
dt2=t3-(t2 > 0 ? t2 : 0); /* Time in 2nd part of hollow cylinder */
if (dt0 < 0) dt0 = 0;
if (dt1 < 0) dt1 = 0;
if (dt2 < 0) dt2 = 0;
/* initialize concentric mode */
if (concentric && !flag_concentric && t0 >= 0
&& VarsInc.shape==0 && thickness>0) {
flag_concentric=1;
}
if (flag_concentric == 1) {
dt1=dt2=0; /* force exit when reaching hole/2nd part */
}
if (!dt0 && !dt2) {
intersect = 0; /* the sample was passed entirely */
break;
}
p_mult = 1;
if (!v) v = sqrt(vx*vx + vy*vy + vz*vz);
if (v) my_a = VarsInc.my_a_v*(2200/v);
else {
printf("Incoherent: %s: ERROR: Null velocity\n",NAME_CURRENT_COMP);
ABSORB; /* should never occur */
}
my_t = my_a + VarsInc.my_s; /* total scattering Xsect (tmp var) */
if (my_t <= 0) {
printf("Incoherent: %s: ERROR: Null total cross section %g. Removing event.\n",
NAME_CURRENT_COMP, my_t);
ABSORB; /* should never occur */
}
d_path = v * (dt0 + dt2); /* Length of full path through sample */
/* Proba of scattering vs absorption (integrating along the whole trajectory) */
ws = VarsInc.my_s/my_t; /* (inc+coh)/(inc+coh+abs) */
/* Proba of transmission along length d_path */
p_trans = exp(-my_t*d_path);
p_scatt = 1 - p_trans; /* portion of beam which scatters */
flag = 0; /* flag used for propagation to exit point before ending */
/* are we next to the exit ? probably no scattering (avoid rounding errors) */
if (VarsInc.my_s*d_path <= 4e-7) {
flag = 1; /* No interaction before the exit */
}
/* force a given fraction of the beam to scatter */
if (p_interact>0 && p_interact<=1) {
/* we force a portion of the beam to interact */
/* This is used to improve statistics on single scattering (and multiple) */
if (!SCATTERED) mc_trans = 1-p_interact;
else mc_trans = 1-p_interact/(4*SCATTERED+1); /* reduce effect on multi scatt */
} else {
mc_trans = p_trans; /* 1 - p_scatt */
}
mc_scatt = 1 - mc_trans; /* portion of beam to scatter (or force to) */
if (mc_scatt <= 0 || mc_scatt>1) flag=1;
/* MC choice: Interaction or transmission ? */
if (!flag && mc_scatt > 0 && (mc_scatt >= 1 || (rand01()) < mc_scatt)) { /* Interaction neutron/sample */
p_mult *= ws; /* Update weight ; account for absorption and retain scattered fraction */
if (!mc_scatt) ABSORB;
/* we have chosen portion mc_scatt of beam instead of p_scatt, so we compensate */
p_mult *= fabs(p_scatt/mc_scatt); /* lower than 1 */
} else {
flag = 1; /* Transmission : no interaction neutron/sample */
if (!mc_trans) ABSORB;
p_mult *= fabs(p_trans/mc_trans); /* attenuate beam by portion which is scattered (and left along) */
}
if (flag) { /* propagate to exit of sample and finish */
intersect = 0;
p *= p_mult; /* apply absorption correction */
PROP_DT(dt0+dt2);
break; /* exit main multi scatt while loop */
}
if (my_t*d_path < 1e-6)
/* For very weak scattering, use simple uniform sampling of scattering
point to avoid rounding errors. */
dt = rand0max(d_path); /* length */
else
dt = -log(1 - rand0max((1 - exp(-my_t*d_path)))) / my_t; /* length */
l_i = dt;/* Penetration in sample: scattering+abs */
dt /= v; /* Time from present position to scattering point */
/* If t0 is in hole, propagate to next part of the hollow cylinder */
if (dt1 > 0 && dt0 > 0 && dt > dt0) dt += dt1;
PROP_DT(dt); /* Point of scattering */
if ((VarsInc.tx || VarsInc.ty || VarsInc.tz)) {
aim_x = VarsInc.tx-x; /* Vector pointing at target (anal./det.) */
aim_y = VarsInc.ty-y;
aim_z = VarsInc.tz-z;
}
if(VarsInc.aw && VarsInc.ah) {
randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
aim_x, aim_y, aim_z, VarsInc.aw, VarsInc.ah, ROT_A_CURRENT_COMP);
} else if(VarsInc.xw && VarsInc.yh) {
randvec_target_rect(&vx, &vy, &vz, &solid_angle,
aim_x, aim_y, aim_z, VarsInc.xw, VarsInc.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);
v_i = v; /* Store initial velocity in case of quasielastic */
E_i = VS2E*v_i*v_i;
if (deltaE==0) {
if (rand01()= order) {
intersect=0; /* reached required number of SCATTERing */
break; /* finish multiple scattering loop */
}
} /* end if intersect */
} while (intersect); /* end do (intersect) (multiple scattering loop) */
%}
MCDISPLAY
%{
if (geometry && strlen(geometry) && strcmp(geometry, "NULL") && strcmp(geometry, "0")) { /* OFF file */
off_display(offdata);
}
else
if (radius > 0 && yheight) { /* cylinder */
cylinder(0,0,0,radius,yheight,0, 0,1,0);
cylinder(0,0,0,radius-thickness,yheight,0, 0,1,0);
}
else if (xwidth && yheight) { /* box/rectangle XY */
box(0,0,0,xwidth, yheight, zdepth);
}
else if (radius > 0 && !yheight) { /* sphere */
sphere(0,0,0,radius,0);
}
%}
END