/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2008, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: SasView_onion * * %Identification * Written by: Jose Robledo * Based on sasmodels from SasView * Origin: FZJ / DTU / ESS DMSC * * * SasView onion model component as sample description. * * %Description * * SasView_onion component, generated from onion.c in sasmodels. * * Example: * SasView_onion(sld_core, radius_core, sld_solvent, n_shells, sld_in[n_shells], sld_out[n_shells], thickness[n_shells], A[n_shells], * model_scale=1.0, model_abs=0.0, xwidth=0.01, yheight=0.01, zdepth=0.005, R=0, * int target_index=1, target_x=0, target_y=0, target_z=1, * focus_xw=0.5, focus_yh=0.5, focus_aw=0, focus_ah=0, focus_r=0, * pd_radius_core=0.0, pd_thickness[n_shells]=0.0) * * %Parameters * INPUT PARAMETERS: * sld_core: [1e-6/Ang^2] ([-inf, inf]) Core scattering length density. * radius_core: [Ang] ([0, inf]) Radius of the core. * sld_solvent: [1e-6/Ang^2] ([-inf, inf]) Solvent scattering length density. * n_shells: [] ([0, 10]) number of shells (must be integer). * sld_in[n_shells]: [1e-6/Ang^2] ([-inf, inf]) scattering length density at the inner radius of shell k. * sld_out[n_shells]: [1e-6/Ang^2] ([-inf, inf]) scattering length density at the outer radius of shell k. * thickness[n_shells]: [Ang] ([0, inf]) Thickness of shell k. * A[n_shells]: [] ([-inf, inf]) Decay rate of shell k. * Optional parameters: * model_abs: [ ] Absorption cross section density at 2200 m/s. * 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. * xwidth: [m] ([-inf, inf]) Horiz. dimension of sample, as a width. * yheight: [m] ([-inf, inf]) vert . dimension of sample, as a height for cylinder/box * zdepth: [m] ([-inf, inf]) depth of sample * R: [m] Outer radius of sample in (x,z) plane for cylinder/sphere. * target_x: [m] relative focus target position. * target_y: [m] relative focus target position. * target_z: [m] relative focus target position. * target_index: [ ] Relative index of component to focus at, e.g. next is +1. * 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. * focus_r: [m] case of circular focusing, focusing radius. * pd_radius_core: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_thickness[n_shells]: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable * * %Link * %End *******************************************************************************/ DEFINE COMPONENT SasView_onion SETTING PARAMETERS ( sld_core=1.0, radius_core=200.0, sld_solvent=6.4, n_shells=1, vector sld_in[n_shells]={1.7}, vector sld_out[n_shells]={2.0}, vector thickness[n_shells]={40.0}, vector A[n_shells]={1.0}, model_scale=1.0, model_abs=0.0, xwidth=0.01, yheight=0.01, zdepth=0.005, R=0, target_x=0, target_y=0, target_z=1, int target_index=1, focus_xw=0.5, focus_yh=0.5, focus_aw=0, focus_ah=0, focus_r=0, pd_radius_core=0.0, pd_thickness[n_shells]=0.0) SHARE %{ %include "sas_kernel_header.c" /* BEGIN Required header for SASmodel onion */ #define HAS_FQ #define FORM_VOL #ifndef SAS_HAVE_sas_3j1x_x #define SAS_HAVE_sas_3j1x_x #line 1 "sas_3j1x_x" /** * Spherical Bessel function 3*j1(x)/x * * Used for low q to avoid cancellation error. * Note that the values differ from sasview ~ 5e-12 rather than 5e-14, but * in this case it is likely cancellation errors in the original expression * using double precision that are the source. */ double sas_3j1x_x(double q); // The choice of the number of terms in the series and the cutoff value for // switching between series and direct calculation depends on the numeric // precision. // // Point where direct calculation reaches machine precision: // // single machine precision eps 3e-8 at qr=1.1 ** // double machine precision eps 4e-16 at qr=1.1 // // Point where Taylor series reaches machine precision (eps), where taylor // series matches direct calculation (cross) and the error at that point: // // prec n eps cross error // single 3 0.28 0.4 6.2e-7 // single 4 0.68 0.7 2.3e-7 // single 5 1.18 1.2 7.5e-8 // double 3 0.01 0.03 2.3e-13 // double 4 0.06 0.1 3.1e-14 // double 5 0.16 0.2 5.0e-15 // // ** Note: relative error on single precision starts increase on the direct // method at qr=1.1, rising from 3e-8 to 5e-5 by qr=1e3. This should be // safe for the sans range, with objects of 100 nm supported to a q of 0.1 // while maintaining 5 digits of precision. For usans/sesans, the objects // are larger but the q is smaller, so again it should be fine. // // See explore/sph_j1c.py for code to explore these ranges. // Use 4th order series #if FLOAT_SIZE>4 #define SPH_J1C_CUTOFF 0.1 #else #define SPH_J1C_CUTOFF 0.7 #endif #pragma acc routine seq double sas_3j1x_x(double q) { // 2017-05-18 PAK - support negative q if (fabs(q) < SPH_J1C_CUTOFF) { const double q2 = q*q; return (1.0 + q2*(-3./30. + q2*(3./840. + q2*(-3./45360.))));// + q2*(3./3991680.))))); } else { double sin_q, cos_q; SINCOS(q, sin_q, cos_q); return 3.0*(sin_q/q - cos_q)/(q*q); } } #endif // SAS_HAVE_sas_3j1x_x #ifndef SAS_HAVE_onion #define SAS_HAVE_onion #line 1 "onion" static double f_exp(double q, double r, double sld_in, double sld_out, double thickness, double A, double side) { const double vol = M_4PI_3 * cube(r); const double qr = q * r; const double bes = sas_3j1x_x(qr); const double alpha = A * r/thickness; double result; if (qr == 0.0) { result = 1.0; } else if (fabs(A) > 0.0) { const double qrsq = qr * qr; const double alphasq = alpha * alpha; const double sumsq = alphasq + qrsq; double sinqr, cosqr; SINCOS(qr, sinqr, cosqr); const double t1 = (alphasq - qrsq)*sinqr/qr - 2.0*alpha*cosqr; const double t2 = alpha*sinqr/qr - cosqr; const double fun = -3.0*(t1/sumsq - t2)/sumsq; const double slope = (sld_out - sld_in)/expm1(A); const double contrast = slope*exp(A*side); result = contrast*fun + (sld_in-slope)*bes; } else { result = sld_in*bes; } return vol * result; } static double outer_radius(double radius_core, double n_shells, double thickness[]) { int n = (int)(n_shells+0.5); double r = radius_core; for (int i=0; i < n; i++) { r += thickness[i]; } return r; } static double form_volume_onion(double radius_core, double n_shells, double thickness[]) { return M_4PI_3*cube(outer_radius(radius_core, n_shells, thickness)); } static double radius_effective_onion(int mode, double radius_core, double n_shells, double thickness[]) { // case 1: outer radius return outer_radius(radius_core, n_shells, thickness); } static void Fq_onion(double q, double *F1, double *F2, double sld_core, double radius_core, double sld_solvent, double n_shells, double sld_in[], double sld_out[], double thickness[], double A[]) { int n = (int)(n_shells+0.5); double r_out = radius_core; double f = f_exp(q, r_out, sld_core, 0.0, 0.0, 0.0, 0.0); for (int i=0; i < n; i++){ const double r_in = r_out; r_out += thickness[i]; f -= f_exp(q, r_in, sld_in[i], sld_out[i], thickness[i], A[i], 0.0); f += f_exp(q, r_out, sld_in[i], sld_out[i], thickness[i], A[i], 1.0); } f -= f_exp(q, r_out, sld_solvent, 0.0, 0.0, 0.0, 0.0); *F1 = 1e-2 * f; *F2 = 1e-4 * f * f; } #endif // SAS_HAVE_onion /* END Required header for SASmodel onion */ %} DECLARE %{ double shape; double my_a_v; %} INITIALIZE %{ shape = -1; /* -1:no shape, 0:cyl, 1:box, 2:sphere */ if (xwidth && yheight && zdepth) shape = 1; else if (R > 0 && yheight) shape = 0; else if (R > 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, my_s; double aim_x = 0, aim_y = 0, aim_z = 1, axis_x, axis_y, axis_z; double arg, tmp_vx, tmp_vy, tmp_vz, vout_x, vout_y, vout_z; double f, solid_angle, vx_i, vy_i, vz_i, q, 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, R, 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, R); } 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); double trace_radius_core = radius_core; double trace_thickness[n_shells] = thickness[n_shells]; if (pd_radius_core != 0.0 || pd_thickness[n_shells] != 0.0) { trace_radius_core = (randnorm () * pd_radius_core + 1.0) * radius_core; trace_thickness[n_shells] = (randnorm () * pd_thickness[n_shells] + 1.0) * thickness[n_shells]; } // Sample dependent. Retrieved from SasView.///////////////////// float Iq_out; Iq_out = 1; double F1 = 0.0, F2 = 0.0; Fq_onion (q, &F1, &F2, sld_core, trace_radius_core, sld_solvent, n_shells, sld_in[n_shells], sld_out[n_shells], trace_thickness[n_shells], A[n_shells]); Iq_out = F2; float vol; vol = 1; // 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 %{ if (shape == 0) { /* cylinder */ circle ("xz", 0, yheight / 2.0, 0, R); circle ("xz", 0, -yheight / 2.0, 0, R); line (-R, -yheight / 2.0, 0, -R, +yheight / 2.0, 0); line (+R, -yheight / 2.0, 0, +R, +yheight / 2.0, 0); line (0, -yheight / 2.0, -R, 0, +yheight / 2.0, -R); line (0, -yheight / 2.0, +R, 0, +yheight / 2.0, +R); } 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, R); circle ("xz", 0, 0.0, 0, R); circle ("yz", 0, 0.0, 0, R); } %} END