/******************************************************************************* * * McXtrace, X-ray tracing package * Copyright, All rights reserved * DTU Physics, Kgs. Lyngby, Denmark * Synchrotron SOLEIL, Saint-Aubin, France * * Component: SasView_pearl_necklace * * %Identification * Written by: Jose Robledo * Based on sasmodels from SasView * Origin: FZJ / DTU / ESS DMSC * * * SasView pearl_necklace model component as sample description. * * %Description * * SasView_pearl_necklace component, generated from pearl_necklace.c in sasmodels. * * Example: * SasView_pearl_necklace(radius, edge_sep, thick_string, num_pearls, sld, sld_string, sld_solvent, * 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=0.0, pd_thick_string=0.0) * * %Parameters * INPUT PARAMETERS: * radius: [Ang] ([0, inf]) Mean radius of the chained spheres. * edge_sep: [Ang] ([0, inf]) Mean separation of chained particles. * thick_string: [Ang] ([0, inf]) Thickness of the chain linkage. * num_pearls: [none] ([1, inf]) Number of pearls in the necklace (must be integer). * sld: [1e-6/Ang^2] ([-inf, inf]) Scattering length density of the chained spheres. * sld_string: [1e-6/Ang^2] ([-inf, inf]) Scattering length density of the chain linkage. * sld_solvent: [1e-6/Ang^2] ([-inf, inf]) Scattering length density of the solvent. * 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: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_thick_string: [] (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_pearl_necklace SETTING PARAMETERS ( radius=80.0, edge_sep=350.0, thick_string=2.5, num_pearls=3, sld=1.0, sld_string=1.0, sld_solvent=6.3, 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=0.0, pd_thick_string=0.0) SHARE %{ %include "sas_kernel_header.c" /* BEGIN Required header for SASmodel pearl_necklace */ #define HAS_Iq #define FORM_VOL #ifndef SAS_HAVE_sas_Si #define SAS_HAVE_sas_Si #line 1 "sas_Si" // integral of sin(x)/x Taylor series approximated to w/i 0.1% double sas_Si(double x); #pragma acc routine seq double sas_Si(double x) { if (x >= M_PI*6.2/4.0) { const double xxinv = 1./(x*x); // Explicitly writing factorial values triples the speed of the calculation const double out_cos = (((-720.*xxinv + 24.)*xxinv - 2.)*xxinv + 1.)/x; const double out_sin = (((-5040.*xxinv + 120.)*xxinv - 6.)*xxinv + 1)*xxinv; double sin_x, cos_x; SINCOS(x, sin_x, cos_x); return M_PI_2 - cos_x*out_cos - sin_x*out_sin; } else { const double xx = x*x; // Explicitly writing factorial values triples the speed of the calculation return (((((-1./439084800.*xx + 1./3265920.)*xx - 1./35280.)*xx + 1./600.)*xx - 1./18.)*xx + 1.)*x; } } #endif // SAS_HAVE_sas_Si #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_pearl_necklace #define SAS_HAVE_pearl_necklace #line 1 "pearl_necklace" double form_volume_pearl_necklace(double radius, double edge_sep, double thick_string, double fp_num_pearls); double Iq_pearl_necklace(double q, double radius, double edge_sep, double thick_string, double fp_num_pearls, double sld, double string_sld, double solvent_sld); // From Igor library static double pearl_necklace_kernel(double q, double radius, double edge_sep, double thick_string, int num_pearls, double sld_pearl, double sld_string, double sld_solv) { // number of string segments const int num_strings = num_pearls - 1; //each masses: contrast * volume const double contrast_pearl = sld_pearl - sld_solv; const double contrast_string = sld_string - sld_solv; const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; const double pearl_vol = M_4PI_3 * radius * radius * radius; const double m_string = contrast_string * string_vol; const double m_pearl = contrast_pearl * pearl_vol; // center to center distance between the neighboring pearls const double A_s = edge_sep + 2.0 * radius; //sine functions of a pearl // Note: lim_(q->0) Si(q*a)/(q*b) = a/b // So therefore: // beta = q==0. ? 1.0 : (Si(q*(A_s-radius)) - Si(q*radius))/q_edge; // gamma = q==0. ? 1.0 : Si(q_edge)/q_edge; // But there is a 1/(1-sinc) term below which blows up so don't bother const double q_edge = q * edge_sep; const double beta = (sas_Si(q*(A_s-radius)) - sas_Si(q*radius)) / q_edge; const double gamma = sas_Si(q_edge) / q_edge; const double psi = sas_3j1x_x(q*radius); // Precomputed sinc terms const double si = sas_sinx_x(q*A_s); const double omsi = 1.0 - si; const double pow_si = pown(si, num_pearls); // form factor for num_pearls const double sss = 2.0*sas_square(m_pearl*psi) * ( - si * (1.0 - pow_si) / (omsi*omsi) + num_pearls / omsi - 0.5 * num_pearls ); // form factor for num_strings (like thin rods) const double srr = m_string * m_string * ( - 2.0 * (1.0 - pow_si/si)*beta*beta / (omsi*omsi) + 2.0 * num_strings*beta*beta / omsi + num_strings * (2.0*gamma - sas_square(sas_sinx_x(q_edge/2.0))) ); // form factor for correlations const double srs = 4.0 * m_string * m_pearl * beta * psi * ( - si * (1.0 - pow_si/si) / (omsi*omsi) + num_strings / omsi ); const double form = sss + srr + srs; return 1.0e-4 * form; } double form_volume_pearl_necklace(double radius, double edge_sep, double thick_string, double fp_num_pearls) { const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls const int num_strings = num_pearls - 1; const double string_vol = edge_sep * M_PI_4 * thick_string * thick_string; const double pearl_vol = M_4PI_3 * radius * radius * radius; const double volume = num_strings*string_vol + num_pearls*pearl_vol; return volume; } static double radius_from_volume_pearl_necklace(double radius, double edge_sep, double thick_string, double fp_num_pearls) { const double vol_tot = form_volume_pearl_necklace(radius, edge_sep, thick_string, fp_num_pearls); return cbrt(vol_tot/M_4PI_3); } static double radius_effective_pearl_necklace(int mode, double radius, double edge_sep, double thick_string, double fp_num_pearls) { switch (mode) { default: case 1: return radius_from_volume_pearl_necklace(radius, edge_sep, thick_string, fp_num_pearls); } } double Iq_pearl_necklace(double q, double radius, double edge_sep, double thick_string, double fp_num_pearls, double sld, double string_sld, double solvent_sld) { const int num_pearls = (int)(fp_num_pearls + 0.5); //Force integer number of pearls const double form = pearl_necklace_kernel(q, radius, edge_sep, thick_string, num_pearls, sld, string_sld, solvent_sld); return form; } #endif // SAS_HAVE_pearl_necklace /* END Required header for SASmodel pearl_necklace */ %} DECLARE %{ double shape; double my_a_k; %} 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; } /*TODO fix absorption*/ my_a_k = model_abs; /* assume absorption is given in 1/m */ %} TRACE %{ double l0, l1, k, l_full, l, dl, d_Phi; double aim_x = 0, aim_y = 0, aim_z = 1, axis_x, axis_y, axis_z; double f, solid_angle, kx_i, ky_i, kz_i, q, qx, qy, qz; char intersect = 0; /* Intersection photon trajectory / sample (sample surface) */ if (shape == 0) { intersect = cylinder_intersect (&l0, &l1, x, y, z, kx, ky, kz, R, yheight); } else if (shape == 1) { intersect = box_intersect (&l0, &l1, x, y, z, kx, ky, kz, xwidth, yheight, zdepth); } else if (shape == 2) { intersect = sphere_intersect (&l0, &l1, x, y, z, kx, ky, kz, R); } if (intersect) { if (l0 < 0) ABSORB; /* Photon enters at l0. */ k = sqrt (kx * kx + ky * ky + kz * kz); l_full = (l1 - l0); /* Length of full path through sample */ dl = rand01 () * (l1 - l0) + l0; /* Point of scattering */ PROP_DL (dl); /* Point of scattering */ l = (dl - l0); /* Penetration in sample */ kx_i = kx; ky_i = ky; kz_i = kz; 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 (&kx, &ky, &kz, &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 (&kx, &ky, &kz, &solid_angle, aim_x, aim_y, aim_z, focus_xw, focus_yh, ROT_A_CURRENT_COMP); } else { randvec_target_circle (&kx, &ky, &kz, &solid_angle, aim_x, aim_y, aim_z, focus_r); } NORM (kx, ky, kz); kx *= k; ky *= k; kz *= k; qx = (kx_i - kx); qy = (ky_i - ky); qz = (kz_i - kz); q = sqrt (qx * qx + qy * qy + qz * qz); double trace_radius = radius; double trace_thick_string = thick_string; if (pd_radius != 0.0 || pd_thick_string != 0.0) { trace_radius = (randnorm () * pd_radius + 1.0) * radius; trace_thick_string = (randnorm () * pd_thick_string + 1.0) * thick_string; } // Sample dependent. Retrieved from SasView.///////////////////// float Iq_out; Iq_out = 1; Iq_out = Iq_pearl_necklace (q, trace_radius, edge_sep, trace_thick_string, num_pearls, sld, sld_string, sld_solvent); float vol; vol = 1; // Scale by 1.0E2 [SasView: 1/cm -> McXtrace: 1/m] Iq_out = model_scale * Iq_out / vol * 1.0E2; p *= l_full * solid_angle / (4 * PI) * Iq_out * exp (-my_a_k * (l + l1)); SCATTER; } %} MCDISPLAY %{ if (shape == 0) { /* cylinder */ cylinder (0, 0, 0, R, yheight, 0, 0, 1, 0); } else if (shape == 1) { /* box */ box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0); } else if (shape == 2) { /* sphere */ sphere (0, 0, 0, R); } %} END