diff --git a/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/README.md b/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/README.md new file mode 100644 index 000000000..1a23822f9 --- /dev/null +++ b/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/README.md @@ -0,0 +1,62 @@ +# The Test_Sans_spheres Instrument + +*McStas Simple instrument used for ensuring that Sans_spheres and SANS_spheres2 remain consistent* + +--- + +## Identification + +- **Instrument:** `Test_Sans_spheres` +- **Author:** Daniel Lomholt Christensen +- **Date:** 18/08/2026 +- **Origin:** Niels Bohr Institute @ UCPH +- **Instrument Site:** `Test_samples` + +### Summary + +A small instrument used to compare the `Sans_spheres` and `SANS_spheres2` components. + + +## Description + +```text +This instrument is a simple SANS setup designed to simulate both the +`Sans_spheres` component and the `SANS_spheres2` component, allowing direct +verification that they produce identical results for this specific test case. +``` + +## Input Parameters + +| Name | Unit | Description | Default | +| ------------------- | ---- | -------------------------------------------------------------------------------- | ------- | +| E\_i | meV | Energy of the simulated neutrons. | - | +| use\_SANS\_spheres2 | 1 | Selects the sample component. `0` uses `Sans_spheres`, `1` uses `SANS_spheres2`. | 0 | +| improved\_res | 1 | Enables the instrument's high-resolution mode. | 0 | +| flux\_mult | n/s | Multiplicative factor applied to the source flux. | 1 | + +--- +## Links + +- [Source code](Test_Sans_spheres.instr) for `Test_SANS.instr`. + +--- + + +## Data treated results + +The results from this instrument have been integrated azimuthally using Pyfai and compared to an analytical +comparison. The below plots show the instrument in the normal (so called "low" resolution mode), +and in the improved\_resolution mode. + +First, the normal resolution mode: +![Normal resolution mode](low_res.png) + + +And secondly the high resolution mode: + +![High resolution mode](high_res.png) + + +Note, that in the high resolution mode, the neutrons will not move their complete distance through the sample, and therefore +when calculating the analytical values, the flight length that is attenuated was approximated to $0.48\cdot l_{\rm full}$ where +$l_{\rm full}$ is the zdepth of the sample diff --git a/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/Test_Sans_spheres.instr b/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/Test_Sans_spheres.instr new file mode 100755 index 000000000..5c258d7eb --- /dev/null +++ b/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/Test_Sans_spheres.instr @@ -0,0 +1,183 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: Test_Sans_spheres +* +* %Identification +* Written by: Daniel Lomholt Christensen +* Date: 18/08/2026 +* Origin: Niels Bohr Institute @ UCPH +* %INSTRUMENT_SITE: Test_samples +* +* A small instrument to compare the Sans_spheres and SANS_spheres2 components +* +* %Description +* This instrument is a small SANS instrument that simulates the Sans_spheres +* component, and the SANS_spheres2 component, in order to ensure that they +* give the same value in this specific case. +* +* %Example: use_SANS_spheres2=0 Detector: psd_det_I=9.73795e+06 +* %Example: use_SANS_spheres2=1 Detector: psd_det_I=9.73795e+06 +* +* +* %Parameters +* E_i: [meV] Energy of the simulated neutrons. +* use_SANS_spheres2: [1] Flag to select whether to use the Sans_spheres component or the SANS_spheres2 component. +* improved_res: [1] Flag to indicate whether the instrument should run in the insanely good resolution mode. +* flux_mult: [n/s] Factor that is multiplied onto the source component in order to adjust the flux. +* +* %Link +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT Test_Sans_spheres(E_i=5, int use_SANS_spheres2=0, int improved_res = 0, flux_mult=0.5512) + +DECLARE +%{ + double focus_w; + double focus_h; + double samp_w; + double samp_h; + double dE; + double src_rad; + double src_focus_w; + double src_focus_h; +%} + +INITIALIZE +%{ + focus_w = 1; + focus_h = 1; + samp_w = 0.01; + samp_h = 0.01; + dE = 0.475; + src_rad = 0.05; + src_focus_w = samp_w; + src_focus_h = samp_h; + if (improved_res){ + samp_w = 0.000001; + samp_h = 0.000001; + dE = 0.00000000000001; + src_rad = 0.0000000001; + src_focus_w = samp_w - samp_w * 0.01; + src_focus_h = samp_h - samp_h * 0.01; + } +%} + +TRACE + +COMPONENT Origin = Progress_bar() + AT (0,0,0) ABSOLUTE + +COMPONENT source = Source_simple( + radius = src_rad, + dist = 10.495, + focus_xw = src_focus_w, + focus_yh = src_focus_h, + E0 = E_i, + dE = dE, + flux = 1e14*flux_mult +) AT (0, 0, 0) RELATIVE Origin +ROTATED (0, 0, 0) RELATIVE Origin + + +COMPONENT post_src_det = PSD_monitor( + nx = 200, + ny = 200, + xwidth = src_rad*2, + yheight = src_rad*2, + restore_neutron = 1 +) AT (0, 0, 0.001) RELATIVE source + + + +COMPONENT post_source_slit = Slit( + radius = 0.005 +) AT (0, 0, 0.025) RELATIVE source + + +COMPONENT between_slits_det = PSD_monitor( + nx = 200, + ny = 200, + xwidth = 0.005, + yheight = 0.005, + restore_neutron = 1 +) AT (0, 0, 0.001) RELATIVE PREVIOUS + + +COMPONENT pre_samp_arm = Arm( + ) AT (0, 0, 10.495) RELATIVE source + + +COMPONENT pre_samp_slit = Slit( + radius = 0.005 +) AT (0, 0, 0) RELATIVE pre_samp_arm + + + +COMPONENT pre_samp_det = PSD_monitor( + nx = 200, + ny = 200, + filename="pre_det", + xwidth = samp_h, + yheight = samp_h, + restore_neutron = 1 +) AT (0, 0, 0.0005) RELATIVE pre_samp_arm + + + + +COMPONENT arm_sample = Arm( + +) AT (0, 0, 0.005) RELATIVE pre_samp_arm +ROTATED (0, 0, 0) RELATIVE pre_samp_arm + +COMPONENT sans_spheres = Sans_spheres( + R = 100, Phi = 1e-2, Delta_rho = 5, sigma_abs = 0, + xwidth=samp_w, yheight=samp_h, zdepth=0.001, + focus_xw=focus_w, + focus_yh=focus_h, + target_index=4 +) WHEN(use_SANS_spheres2 == 0) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT sans_spheres2 = SANS_spheres2( + xwidth = samp_w, yheight = samp_h, zthick = 0.001, + dsdw_inc = 0, singlesp = 1, + R = 100, phi = 1e-2, drho = 5e11, + Qmind = 0.001, Qmaxd = 100.25, + sc_aim =1, sans_aim = 1 +) WHEN (use_SANS_spheres2 == 1) +AT (0,0,0) RELATIVE PREVIOUS + + + +COMPONENT total_scat = PSD_monitor_4PI( + nx = 360, + ny = 360, + radius=0.1, + restore_neutron=1, + filename="total_scat" +) //WHEN (total_scattering==1) +AT (0, 0, 0) RELATIVE arm_sample + + +COMPONENT arm_det = Arm( +) AT (0, 0, 0) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +COMPONENT psd_det = PSD_monitor( + nx = 1000, + ny = 1000, + filename="post_det", + xwidth = 1, + yheight = 1, + restore_neutron = 1 +) AT (0, 0, 5) RELATIVE arm_det + + +FINALLY +%{ +%} +END + diff --git a/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/high_res.png b/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/high_res.png new file mode 100644 index 000000000..20c2ae2e2 Binary files /dev/null and b/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/high_res.png differ diff --git a/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/low_res.png b/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/low_res.png new file mode 100644 index 000000000..b6936c224 Binary files /dev/null and b/mcstas-comps/examples/Tests_samples/Test_Sans_spheres/low_res.png differ diff --git a/mcstas-comps/samples/SANS_spheres2.comp b/mcstas-comps/samples/SANS_spheres2.comp index 8dead0a0d..77b6c58c1 100644 --- a/mcstas-comps/samples/SANS_spheres2.comp +++ b/mcstas-comps/samples/SANS_spheres2.comp @@ -343,6 +343,9 @@ TRACE intersect = box_intersect (&t0, &t1, x, y, z, vx, vy, vz, xwidth, yheight, zthick); if (!intersect) ABSORB; + if (singlesp == 1) { + p *= exp (-v * t1 * Stot); // Attenuate beam if no multiple scattering is simulated + } PROP_DT (t1); } else { diff --git a/mcstas-comps/samples/Sans_spheres.comp b/mcstas-comps/samples/Sans_spheres.comp index 58aa21231..ccc5099e5 100644 --- a/mcstas-comps/samples/Sans_spheres.comp +++ b/mcstas-comps/samples/Sans_spheres.comp @@ -49,8 +49,8 @@ * * Variables calculated in the component * -* my_s: Attenuation factor due to scattering [m^-1] -* my_a: Attenuation factor due to absorbtion [m^-1] +* mu_s: Attenuation factor due to scattering [m^-1] +* mu_a: Attenuation factor due to absorbtion [m^-1] * * %Link * The test/example instrument SANS.instr. @@ -67,11 +67,60 @@ 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) /* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ +SHARE %{ + + #pragma acc routine seq + double + Mind (double A, double B) { + if (A < B) + return A; + else + return B; + }; + + #pragma acc routine seq + double + Maxd (double A, double B) { + if (A > B) + return A; + else + return B; + }; + + #pragma acc routine seq + int + MinI (int A, int B) { + if (A < B) + return A; + else + return B; + }; + + #pragma acc routine seq + int + MaxI (int A, int B) { + if (A > B) + return A; + else + return B; + }; + + + double dSigdOmega(double R, double Phi, double Delta_rho, double q){ + + double f = 3 * (sin (q * R) - q * R * cos (q * R)) / (q * R * q * R * q * R); + double V = 4 * PI * R * R * R / 3 ; + + return Phi * Delta_rho * Delta_rho * V * f*f; + } +%} DECLARE %{ - double my_s_pre; - double my_a_v; + double mu_a_v; double shape; + DArray1d Idsdw; + double Qminl; + double Qmaxl; %} INITIALIZE %{ @@ -103,14 +152,50 @@ INITIALIZE printf ("Sans_spheres: %s: The target is not defined. Using direct beam (Z-axis).\n", NAME_CURRENT_COMP); target_z = 1; } + if (focus_ah || focus_aw) { // Convert focus from degrees to radians + focus_aw *= PI / 180; + focus_ah *= PI / 180; + } + + mu_a_v = sigma_abs * 2200; /* Is not yet divided by v.*/ + + // Perform integral of dsigdomega over all of omega for use in adjusting the weight of the neutron. + int iii, kkk; + double Qmind = 0.0000001; + double Qmaxd = 100; + + Qminl = log10 (Qmind); + Qmaxl = log10 (Qmaxd); + + double q, Isq; + double qmin, qmax, step; + int istp; + + istp = floor ((Qmaxl - Qminl) * 300.0 + 0.5); - my_a_v = sigma_abs * 2200; /* Is not yet divided by v.*/ + Idsdw = create_darr1d (31); - my_s_pre = Phi * 4 * PI * R * R * R / 3 * Delta_rho * Delta_rho; + // Use the integration over Q for 30 different wavelengths, taken from SANS_spheres2.comp + /* By integration, calculate the coherent scattering cross-section for the relevant wavelength range */ + for (iii = 1; iii <= 30; iii++) { /* wavelength in AA, up to 30 */ + Idsdw[iii] = 0.0; + Isq = 0.0; + qmin = 0.0; + step = (log10 (Mind (Qmaxd, 4.0 * PI / iii)) - Qminl) / istp; + for (kkk = 0; kkk <= istp; kkk++) { + qmax = pow (10.0, Qminl + kkk * step); + q = 0.5 * (qmin + qmax); + Isq += dSigdOmega (R, Phi, Delta_rho, q) * q * (qmax - qmin); + qmin = qmax; + }; + Idsdw[iii] = Isq; + }; %} TRACE %{ - double t0, t1, v, l_full, l, l_1, dt, my_s; + double t0, t1, v, l_full, l_in, l_out, dt, mu_s, Scoh; + double k0, lambda; + int Ilam, Ilam2; double aim_x = 0, aim_y = 0, aim_z = 1; double f, solid_angle, vx_i, vy_i, vz_i, qx, qy, qz, q; char intersect = 0; @@ -130,7 +215,20 @@ TRACE 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 */ + l_in = v * (dt - t0); /* Penetration in sample */ + k0 = V2K * v; + lambda = 2.0 * PI / k0; + Ilam = MaxI (floor (lambda), 1); + Ilam2 = MinI (Ilam + 1, 30); + /* Coherent "SANS" scattering - in 3 intervals, asymptotic values at the low and high WL end */ + if (lambda <= 1.0) + Scoh = 2 * PI * Idsdw[1] / (k0 * k0); + else { + if (lambda >= 30.0) + Scoh = 2 * PI * Idsdw[30] / (k0 * k0); + else + Scoh = 2 * PI * ((Ilam2 - lambda) * Idsdw[Ilam] + (lambda - Ilam) * Idsdw[Ilam2]) / (k0 * k0); + }; vx_i = vx; vy_i = vy; @@ -155,10 +253,18 @@ TRACE qy = V2K * (vy_i - vy); qz = V2K * (vz_i - vz); q = sqrt (qx * qx + qy * qy + qz * qz); + 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); + double l_out = v * t1; + PROP_DT (t1); - f = 3 * (sin (q * R) - q * R * cos (q * R)) / (q * R * q * R * q * R); + mu_s = dSigdOmega (R, Phi, Delta_rho, q); - double pmul = l_full * solid_angle / (4 * PI) * my_s_pre * f * f * exp (-(my_a_v / v) * l_full); + double pmul = l_full * solid_angle * mu_s * exp (-(mu_a_v / v + Scoh) * (l_in + l_out)); p = p * pmul; SCATTER; }