Skip to content
Permalink
Branch: master
Find file Copy path
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
295 lines (273 sloc) 11.3 KB
/*******************************************************************************
*
* Component: SasView_model
*
* %I
* Written by: Jakob Garde, Torben Nielsen, Peter Willendrup
* Date: 03.02.2016
* Version: 0.5
* Origin: SasView, DTU, European Spallation Source ERIC
* Release: McStas 2.3
*
* This SANS sample exposes <a href="http://www.sasview.org">SasView's</a> scattering kernels to McStas. In this way SasView's monodisperse scattering kernels can be call from McStas.
*
* %D
* Sample for use in SANS instruments. The models describe mono disperse particles in thin solution. The sample geometry may have the shape:
*
* Shape: - A filled box with dimensions xwidth, yheight and zdepth.
* - A cylinder with dimensions radius and yheight.
* - A filled sphere given by radius.
*
* These parameters are mutually exclusive.
*
* Example using spheres in thin solution, with radius=200 AA and a delta_sld=0.6 fm/AA^3:
* SasView_model(model_index=47, model_scale=1.0, model_pars={1, 7, 200}, model_abs=0.0,
xwidth=0.01, yheight=0.01, zdepth=0.005,)
*
*
* The algorithm of this component requires use of the ISO C standard c99. You could use MCSTAS_CFLAGS like MCSTAS_CFLAGS = -std=c99 -g -O2 -lm.
* Or, add the flag "-std=c99" to "C flags" in the mcgui or mcgui-py config.
*
* The list of scattering models in SasView is called sasmodels. The list of McStas available <a href="http://www.sasview.org/sasview/user/models/model_functions.html#model">SasView sasmodels</a> are found in the table below.
*
* A few models may require manual documentation lookup using the above link to the SasView site.
*
* McStas does currently not support multiplication of formfactor models with structure factor models.
*
* The 2D scattering scattering kernels are denoted by modelname_xy. I.e. to evalulate scattering from aligned cylinders use model_index=10 to use cylinder_xy.
*
* MDOC
* 0 - None
* 1 - barbell (sld, solvent_sld, bell_radius, radius, length)
* 2 - barbell_xy (sld, solvent_sld, bell_radius, radius, length, theta, phi)
* 3 - bcc_paracrystal (dnn, d_factor, radius, sld, solvent_sld)
* 4 - bcc_paracrystal_xy (dnn, d_factor, radius, sld, solvent_sld, theta, phi, psi)
* 5 - capped_cylinder (sld, solvent_sld, radius, cap_radius, length)
* 6 - capped_cylinder_xy (sld, solvent_sld, radius, cap_radius, length, theta, phi)
* 7 - core_shell_cylinder (core_sld, shell_sld, solvent_sld, radius, thickness, length)
* 8 - core_shell_cylinder_xy (core_sld, shell_sld, solvent_sld, radius, thickness, length, theta, phi)
* 9 - cylinder (sld, solvent_sld, radius, length)
* 10 - cylinder_xy (sld, solvent_sld, radius, length, theta, phi)
* 11 - dab (length)
* 12 - dab_xy (length)
* 13 - ellipsoid (sld, solvent_sld, rpolar, requatorial)
* 14 - ellipsoid_xy (sld, solvent_sld, rpolar, requatorial, theta, phi)
* 15 - fcc_paracrystal (dnn, d_factor, radius, sld, solvent_sld)
* 16 - fcc_paracrystal_xy (dnn, d_factor, radius, sld, solvent_sld, theta, phi, psi)
* 17 - gaussian_peak (q0, sigma)
* 18 - gaussian_peak_xy (q0, sigma)
* 19 - hardsphere (effect_radius, volfraction)
* 20 - hardsphere_xy (effect_radius, volfraction)
* 21 - HayterMSAsq (effect_radius, zz, VolFrac, Temp, csalt, dialec)
* 22 - HayterMSAsq_xy (effect_radius, charge, volfraction, temperature, saltconc, dielectconst)
* 23 - lamellar (sld, solvent_sld, thickness)
* 24 - lamellar_xy (sld, solvent_sld, thickness)
* 25 - lamellarCailleHG (tail_length, head_length, Nlayers, dd, Cp, tail_sld, head_sld, solvent_sld)
* 26 - lamellarCailleHG_xy (tail_length, head_length, Nlayers, spacing, Caille_parameter, sld, head_sld, solvent_sld)
* 27 - lamellarPC (th, Nlayers, davg, pd, sld, solvent_sld)
* 28 - lamellarPC_xy (thickness, Nlayers, spacing, spacing_polydisp, sld, solvent_sld)
* 29 - parallelepiped (sld, solvent_sld, a_side, b_side, c_side)
* 30 - parallelepiped_xy (sld, solvent_sld, a_side, b_side, c_side, theta, phi, psi)
* 31 - sphere (sld, solvent_sld, radius)
* 32 - sphere_xy (sld, solvent_sld, radius)
* 33 - stickyhardsphere (effect_radius, volfraction, perturb, stickiness)
* 34 - stickyhardsphere_xy (effect_radius, volfraction, perturb, stickiness)
* 35 - triaxial_ellipsoid (sld, solvent_sld, req_minor, req_major, rpolar)
* 36 - triaxial_ellipsoid_xy (sld, solvent_sld, req_minor, req_major, rpolar, theta, phi, psi)
* MDOC_END
*
* %P
* Definition parameters:
*
* model_index: Index of the applied sasview model. Recompile instrument for changes to take effect.
* model_scale: Global scale factor for scattering kernel. For systems without inter-particle interference, the form factors can be related
* to the scattering intensity by the particle volume fraction.
* model_pars: Model parameters are given as a set of comma-separated values enclosed by {}, e.g. {60} for the model index 21 (the guinier model).
* Consult the sasview docs for further info.
* model_abs: Absorption cross section density at 2200 m/s (1/m)
*
* Input parameters:
*
* radius: Outer radius of sample in (x,z) plane for cylinder/sphere [m]
* xwidth: horiz. dimension of sample, as a width (m)
* yheight: vert . dimension of sample, as a height for cylinder/box (m)
* zdepth: depth of sample (m)
* target_index: Relative index of component to focus at, e.g. next is +1 [1]
* 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]
* focus_r: Detector (disk-shaped) radius (m)
*
* Optional parameters:
*
* target_x: -\
* target_y: 3- relative focus target position [m]
* target_z: -/
*
*
* %Link
* http://www.sasview.org/sasview/user/models/model_functions.html
*
* %E
*******************************************************************************/
DEFINE COMPONENT SasView_model
DEFINITION PARAMETERS (
model_index=21,
model_scale= 1.0,
model_pars={60},
model_abs=0.5
)
SETTING PARAMETERS (
xwidth=0,
yheight=0,
zdepth=0,
radius=0,
target_x=0,
target_y=0,
target_z=6,
int target_index=0,
focus_xw=0,
focus_yh=0,
focus_aw=0,
focus_ah=0,
focus_r=0)
OUTPUT PARAMETERS (
my_a_v,
shape)
/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */
DEPENDENCY "-std=c99"
DECLARE
%{
%include "sasview_proxy.c"
float model_pars_ptr[] = model_pars;
double my_a_v, q, shape;
%}
INITIALIZE
%{
shape=-1; /* -1:no shape, 0:cyl, 1:box, 2:sphere */
if (xwidth && yheight && zdepth)
shape=1;
else if (radius > 0 && yheight)
shape=0;
else if (radius > 0 && !yheight)
shape=2;
if (shape < 0)
exit(fprintf(stderr, "SasView_model: %s: sample has invalid dimensions.\n"
"ERROR Please check parameter values.\n", NAME_CURRENT_COMP));
/* now compute target coords if a component index is supplied */
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, &target_x, &target_y, &target_z);
}
if (!(target_x || target_y || target_z)) {
printf("SasView_model: %s: The target is not defined. Using direct beam (Z-axis).\n",
NAME_CURRENT_COMP);
target_z=1;
}
my_a_v = model_abs*2200*100; /* Is not yet divided by v. 100: Convert barns -> fm^2 */
%}
TRACE
%{
double t0, t1, v, l_full, l, l_1, dt, d_phi, theta;
double aim_x=0, aim_y=0, aim_z=1, axis_x, axis_y, axis_z;
double arg, tmp_vx, tmp_vy, tmp_vz;
double f, solid_angle, vx_i, vy_i, vz_i, qx, qy, qz;
char intersect=0;
/* Intersection neutron trajectory / sample (sample surface) */
if (shape == 0)
intersect = cylinder_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius, yheight);
else if (shape == 1)
intersect = box_intersect(&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zdepth);
else if (shape == 2)
intersect = sphere_intersect(&t0, &t1, x, y, z, vx, vy, vz, radius);
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) + t0; /* Time of scattering */
PROP_DT(dt); /* Point of scattering */
l = v*(dt-t0); /* Penetration in sample */
vx_i=vx;
vy_i=vy;
vz_i=vz;
if ((target_x || target_y || target_z)) {
aim_x = target_x-x; /* Vector pointing at target (anal./det.) */
aim_y = target_y-y;
aim_z = target_z-z;
}
if(focus_aw && focus_ah) {
randvec_target_rect_angular(&vx, &vy, &vz, &solid_angle,
aim_x, aim_y, aim_z, focus_aw, focus_ah, ROT_A_CURRENT_COMP);
} else if(focus_xw && focus_yh) {
randvec_target_rect(&vx, &vy, &vz, &solid_angle,
aim_x, aim_y, aim_z, focus_xw, focus_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);
vx *= v;
vy *= v;
vz *= v;
qx = V2K*(vx_i-vx);
qy = V2K*(vy_i-vy);
qz = V2K*(vz_i-vz);
q = sqrt(qx*qx+qy*qy+qz*qz);
float Iq_out;
Iq_out = getIq(q, qx, qy, model_pars_ptr);
float vol;
vol=getFormVol(model_pars_ptr);
// Scale by 1.0E2 [SasView: 1/cm -> McStas: 1/m]
Iq_out = model_scale*Iq_out / vol * 1.0E2;
l_1 = v*t1;
p *= l_full*solid_angle/(4*PI)*Iq_out*exp(-my_a_v*(l+l_1)/v);
SCATTER;
}
%}
MCDISPLAY
%{
magnify("xyz");
if (shape == 0) { /* cylinder */
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);
}
else if (shape == 1) { /* box */
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);
}
else if (shape == 2) { /* sphere */
circle("xy", 0, 0.0, 0, radius);
circle("xz", 0, 0.0, 0, radius);
circle("yz", 0, 0.0, 0, radius);
}
%}
END
You can’t perform that action at this time.