/******************************************************************************* * * McStas, neutron ray-tracing package * Copyright (C) 1997-2008, All rights reserved * Risoe National Laboratory, Roskilde, Denmark * Institut Laue Langevin, Grenoble, France * * Component: SasView_core_shell_parallelepiped * * %Identification * Written by: Jose Robledo * Based on sasmodels from SasView * Origin: FZJ / DTU / ESS DMSC * * * SasView core_shell_parallelepiped model component as sample description. * * %Description * * SasView_core_shell_parallelepiped component, generated from core_shell_parallelepiped.c in sasmodels. * * Example: * SasView_core_shell_parallelepiped_aniso(sld_core, sld_a, sld_b, sld_c, sld_solvent, length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c, theta, phi, Psi, * 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_length_a=0.0, pd_length_b=0.0, pd_length_c=0.0, pd_thick_rim_a=0.0, pd_thick_rim_b=0.0, pd_thick_rim_c=0.0, pd_theta=0.0, pd_phi=0.0, pd_Psi=0.0) * * %Parameters * INPUT PARAMETERS: * sld_core: [1e-6/Ang^2] ([-inf, inf]) Parallelepiped core scattering length density. * sld_a: [1e-6/Ang^2] ([-inf, inf]) Parallelepiped A rim scattering length density. * sld_b: [1e-6/Ang^2] ([-inf, inf]) Parallelepiped B rim scattering length density. * sld_c: [1e-6/Ang^2] ([-inf, inf]) Parallelepiped C rim scattering length density. * sld_solvent: [1e-6/Ang^2] ([-inf, inf]) Solvent scattering length density. * length_a: [Ang] ([0, inf]) Shorter side of the parallelepiped. * length_b: [Ang] ([0, inf]) Second side of the parallelepiped. * length_c: [Ang] ([0, inf]) Larger side of the parallelepiped. * thick_rim_a: [Ang] ([0, inf]) Thickness of A rim. * thick_rim_b: [Ang] ([0, inf]) Thickness of B rim. * thick_rim_c: [Ang] ([0, inf]) Thickness of C rim. * 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_length_a: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_length_b: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_length_c: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_thick_rim_a: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_thick_rim_b: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_thick_rim_c: [] (0,inf) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_theta: [] (0,360) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_phi: [] (0,360) defined as (dx/x), where x is de mean value and dx the standard devition of the variable. * pd_Psi: [] (0,360) defined as (dx/x), where x is de mean value and dx the standard devition of the variable * * %Link * %End *******************************************************************************/ DEFINE COMPONENT SasView_core_shell_parallelepiped_aniso SETTING PARAMETERS ( sld_core=1, sld_a=2, sld_b=4, sld_c=2, sld_solvent=6, length_a=35, length_b=75, length_c=400, thick_rim_a=10, thick_rim_b=10, thick_rim_c=10, theta=0, phi=0, Psi=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_length_a=0.0, pd_length_b=0.0, pd_length_c=0.0, pd_thick_rim_a=0.0, pd_thick_rim_b=0.0, pd_thick_rim_c=0.0, pd_theta=0.0, pd_phi=0.0, pd_Psi=0.0) SHARE %{ %include "sas_kernel_header.c" /* BEGIN Required header for SASmodel core_shell_parallelepiped */ #define HAS_Iqabc #define HAS_FQ #define FORM_VOL #ifndef SAS_HAVE_gauss76 #define SAS_HAVE_gauss76 #line 1 "gauss76" // Created by Andrew Jackson on 4/23/07 #ifdef GAUSS_N # undef GAUSS_N # undef GAUSS_Z # undef GAUSS_W #endif #define GAUSS_N 76 #define GAUSS_Z Gauss76Z #define GAUSS_W Gauss76Wt // Gaussians constant double Gauss76Wt[76] = { .00126779163408536, //0 .00294910295364247, .00462793522803742, .00629918049732845, .00795984747723973, .00960710541471375, .0112381685696677, .0128502838475101, .0144407317482767, .0160068299122486, .0175459372914742, //10 .0190554584671906, .020532847967908, .0219756145344162, .0233813253070112, .0247476099206597, .026072164497986, .0273527555318275, .028587223650054, .029773487255905, .0309095460374916, //20 .0319934843404216, .0330234743977917, .0339977794120564, .0349147564835508, .0357728593807139, .0365706411473296, .0373067565423816, .0379799643084053, .0385891292645067, .0391332242205184, //30 .0396113317090621, .0400226455325968, .040366472122844, .0406422317102947, .0408494593018285, .040987805464794, .0410570369162294, .0410570369162294, .040987805464794, .0408494593018285, //40 .0406422317102947, .040366472122844, .0400226455325968, .0396113317090621, .0391332242205184, .0385891292645067, .0379799643084053, .0373067565423816, .0365706411473296, .0357728593807139, //50 .0349147564835508, .0339977794120564, .0330234743977917, .0319934843404216, .0309095460374916, .029773487255905, .028587223650054, .0273527555318275, .026072164497986, .0247476099206597, //60 .0233813253070112, .0219756145344162, .020532847967908, .0190554584671906, .0175459372914742, .0160068299122486, .0144407317482767, .0128502838475101, .0112381685696677, .00960710541471375, //70 .00795984747723973, .00629918049732845, .00462793522803742, .00294910295364247, .00126779163408536 //75 (indexed from 0) }; constant double Gauss76Z[76] = { -.999505948362153, //0 -.997397786355355, -.993608772723527, -.988144453359837, -.981013938975656, -.972229228520377, -.961805126758768, -.949759207710896, -.936111781934811, -.92088586125215, -.904107119545567, //10 -.885803849292083, -.866006913771982, -.844749694983342, -.822068037328975, -.7980001871612, -.77258672828181, -.74587051350361, -.717896592387704, -.688712135277641, -.658366353758143, //20 -.626910417672267, -.594397368836793, -.560882031601237, -.526420920401243, -.491072144462194, -.454895309813726, -.417951418780327, -.380302767117504, -.342012838966962, -.303146199807908, //30 -.263768387584994, -.223945802196474, -.183745593528914, -.143235548227268, -.102483975391227, -.0615595913906112, -.0205314039939986, .0205314039939986, .0615595913906112, .102483975391227, //40 .143235548227268, .183745593528914, .223945802196474, .263768387584994, .303146199807908, .342012838966962, .380302767117504, .417951418780327, .454895309813726, .491072144462194, //50 .526420920401243, .560882031601237, .594397368836793, .626910417672267, .658366353758143, .688712135277641, .717896592387704, .74587051350361, .77258672828181, .7980001871612, //60 .822068037328975, .844749694983342, .866006913771982, .885803849292083, .904107119545567, .92088586125215, .936111781934811, .949759207710896, .961805126758768, .972229228520377, //70 .981013938975656, .988144453359837, .993608772723527, .997397786355355, .999505948362153 //75 }; #pragma acc declare copyin(Gauss76Wt[0:76], Gauss76Z[0:76]) #endif // SAS_HAVE_gauss76 #ifndef SAS_HAVE_core_shell_parallelepiped #define SAS_HAVE_core_shell_parallelepiped #line 1 "core_shell_parallelepiped" // Set OVERLAPPING to 1 in order to fill in the edges of the box, with // c endcaps and b overlapping a. With the proper choice of parameters, // (setting rim slds to sld, core sld to solvent, rim thickness to thickness // and subtracting 2*thickness from length, this should match the hollow // rectangular prism.) Set it to 0 for the documented behaviour. #define OVERLAPPING 0 static double form_volume_core_shell_parallelepiped(double length_a, double length_b, double length_c, double thick_rim_a, double thick_rim_b, double thick_rim_c) { return #if OVERLAPPING // Hollow rectangular prism only includes the volume of the shell // so uncomment the next line when comparing. Solid rectangular // prism, or parallelepiped want filled cores, so comment when // comparing. //-length_a * length_b * length_c + (length_a + 2.0*thick_rim_a) * (length_b + 2.0*thick_rim_b) * (length_c + 2.0*thick_rim_c); #else length_a * length_b * length_c + 2.0 * thick_rim_a * length_b * length_c + 2.0 * length_a * thick_rim_b * length_c + 2.0 * length_a * length_b * thick_rim_c; #endif } static double radius_from_excluded_volume_core_shell_parallelepiped(double length_a, double length_b, double length_c, double thick_rim_a, double thick_rim_b, double thick_rim_c) { double r_equiv, length; double lengths[3] = {length_a+thick_rim_a, length_b+thick_rim_b, length_c+thick_rim_c}; double lengthmax = fmax(lengths[0],fmax(lengths[1],lengths[2])); double length_1 = lengthmax; double length_2 = lengthmax; double length_3 = lengthmax; for(int ilen=0; ilen<3; ilen++) { if (lengths[ilen] < length_1) { length_2 = length_1; length_1 = lengths[ilen]; } else { if (lengths[ilen] < length_2) { length_2 = lengths[ilen]; } } } if(length_2-length_1 > length_3-length_2) { r_equiv = sqrt(length_2*length_3/M_PI); length = length_1; } else { r_equiv = sqrt(length_1*length_2/M_PI); length = length_3; } return 0.5*cbrt(0.75*r_equiv*(2.0*r_equiv*length + (r_equiv + length)*(M_PI*r_equiv + length))); } static double radius_from_volume_core_shell_parallelepiped(double length_a, double length_b, double length_c, double thick_rim_a, double thick_rim_b, double thick_rim_c) { const double volume = form_volume_core_shell_parallelepiped(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); return cbrt(volume/M_4PI_3); } static double radius_from_crosssection_core_shell_parallelepiped(double length_a, double length_b, double thick_rim_a, double thick_rim_b) { const double area_xsec_paral = length_a*length_b + 2.0*thick_rim_a*length_b + 2.0*thick_rim_b*length_a; return sqrt(area_xsec_paral/M_PI); } static double radius_effective_core_shell_parallelepiped(int mode, double length_a, double length_b, double length_c, double thick_rim_a, double thick_rim_b, double thick_rim_c) { switch (mode) { default: case 1: // equivalent cylinder excluded volume return radius_from_excluded_volume_core_shell_parallelepiped(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); case 2: // equivalent volume sphere return radius_from_volume_core_shell_parallelepiped(length_a, length_b, length_c, thick_rim_a, thick_rim_b, thick_rim_c); case 3: // half outer length a return 0.5 * length_a + thick_rim_a; case 4: // half outer length b return 0.5 * length_b + thick_rim_b; case 5: // half outer length c return 0.5 * length_c + thick_rim_c; case 6: // equivalent circular cross-section return radius_from_crosssection_core_shell_parallelepiped(length_a, length_b, thick_rim_a, thick_rim_b); case 7: // half outer ab diagonal return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b)); case 8: // half outer diagonal return 0.5*sqrt(square(length_a+ 2.0*thick_rim_a) + square(length_b+ 2.0*thick_rim_b) + square(length_c+ 2.0*thick_rim_c)); } } static void Fq_core_shell_parallelepiped(double q, double *F1, double *F2, double core_sld, double arim_sld, double brim_sld, double crim_sld, double solvent_sld, double length_a, double length_b, double length_c, double thick_rim_a, double thick_rim_b, double thick_rim_c) { // Code converted from functions CSPPKernel and CSParallelepiped in libCylinder.c // Did not understand the code completely, it should be rechecked (Miguel Gonzalez) // Code is rewritten, the code is compliant with Diva Singh's thesis now (Dirk Honecker) // Code rewritten; cross checked against hollow rectangular prism and realspace (PAK) const double half_q = 0.5*q; const double tA = length_a + 2.0*thick_rim_a; const double tB = length_b + 2.0*thick_rim_b; const double tC = length_c + 2.0*thick_rim_c; // Scale factors const double dr0 = (core_sld-solvent_sld); const double drA = (arim_sld-solvent_sld); const double drB = (brim_sld-solvent_sld); const double drC = (crim_sld-solvent_sld); // outer integral (with gauss points), integration limits = 0, 1 // substitute d_cos_alpha for sin_alpha d_alpha double outer_sum_F1 = 0; //initialize integral double outer_sum_F2 = 0; //initialize integral for( int i=0; i 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_length_a = length_a; double trace_length_b = length_b; double trace_length_c = length_c; double trace_thick_rim_a = thick_rim_a; double trace_thick_rim_b = thick_rim_b; double trace_thick_rim_c = thick_rim_c; if (pd_length_a != 0.0 || pd_length_b != 0.0 || pd_length_c != 0.0 || pd_thick_rim_a != 0.0 || pd_thick_rim_b != 0.0 || pd_thick_rim_c != 0.0) { trace_length_a = (randnorm () * pd_length_a + 1.0) * length_a; trace_length_b = (randnorm () * pd_length_b + 1.0) * length_b; trace_length_c = (randnorm () * pd_length_c + 1.0) * length_c; trace_thick_rim_a = (randnorm () * pd_thick_rim_a + 1.0) * thick_rim_a; trace_thick_rim_b = (randnorm () * pd_thick_rim_b + 1.0) * thick_rim_b; trace_thick_rim_c = (randnorm () * pd_thick_rim_c + 1.0) * thick_rim_c; } double trace_theta = theta, dtheta = 0.0; double trace_phi = phi, dphi = 0.0; double trace_Psi = Psi, dPsi = 0.0; if (pd_theta != 0.0 || pd_phi != 0.0 || pd_Psi != 0.0) { trace_theta = ((rand01 () - 0.5) * pd_theta + 1.0) * theta; dtheta = trace_theta - theta; trace_phi = ((rand01 () - 0.5) * pd_phi + 1.0) * phi; dphi = trace_phi - phi; trace_Psi = ((rand01 () - 0.5) * pd_Psi + 1.0) * Psi; dPsi = trace_Psi - Psi; } // Sample dependent. Retrieved from SasView.///////////////////// float Iq_out; Iq_out = 1; double qa = 0.0, qb = 0.0, qc = 0.0; QABCRotation rotation; qabc_rotation (&rotation, trace_theta, trace_phi, trace_Psi, dtheta, dphi, dPsi); qabc_apply (&rotation, qx, qy, &qa, &qb, &qc); Iq_out = Iqabc_core_shell_parallelepiped (qa, qb, qc, sld_core, sld_a, sld_b, sld_c, sld_solvent, trace_length_a, trace_length_b, trace_length_c, trace_thick_rim_a, trace_thick_rim_b, trace_thick_rim_c); 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