PhononSimple_process, Phonon_simple: scale the one-phonon cross-section to Squires 3.120

mcstas-comps/examples/Tests_samples/Samples_Phonon/README.md

@@ -18,8 +18,8 @@ Refer to the component documentation for further instructions.
 
 ## Examples
 
-- **Test: E=10 -n 1e5 focus_r=0.001 Detector: mon1_I=2.86265e-25**
-- **Test: E=10 -n 1e5 focus_a=0.1 Detector: mon1_I=2.86265e-25**
+- **Test: E=10 -n 1e5 focus_r=0.001 Detector: mon1_I=3.62846e-24**
+- **Test: E=10 -n 1e5 focus_a=0.1 Detector: mon1_I=3.62846e-24**
 
 ## Input parameters
 

mcstas-comps/examples/Tests_samples/Samples_Phonon/Samples_Phonon.instr

@@ -17,8 +17,8 @@
 * Simple test instrument for the Phonon_simple component.
 * Refer to the component documentation for further instructions.
 *
-* %Example: E=10 -n 1e5 focus_r=0.001 Detector: mon1_I=2.86265e-25
-* %Example: E=10 -n 1e5 focus_a=0.1 Detector: mon1_I=2.86265e-25
+* %Example: E=10 -n 1e5 focus_r=0.001 Detector: mon1_I=3.62846e-24
+* %Example: E=10 -n 1e5 focus_a=0.1 Detector: mon1_I=3.62846e-24
 *
 * %Parameters
 * E:    [meV] Mean energy at source   

mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/README.md

@@ -55,9 +55,9 @@ generate_dispersion_files.py (Python with numpy).
 
 ## Examples
 
-- **Test: h=1.2 l=1 dE=5 sample_choice=0 Detector: detector_I=1.46e-11**
-- **Test: h=1.2 l=1 dE=5 sample_choice=1 Detector: detector_I=1.41e-11**
-- **Test: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.46e-11**
+- **Test: h=1.2 l=1 dE=5 sample_choice=0 Detector: detector_I=1.85e-10**
+- **Test: h=1.2 l=1 dE=5 sample_choice=1 Detector: detector_I=1.79e-10**
+- **Test: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.85e-10**
 - **Test: h=1.15 l=0 dE=3 T=2 sample_choice=2 Detector: detector_I=8.91e-10**
 - **Test: h=1.15 l=0 dE=3 T=2 sample_choice=3 Detector: detector_I=9.41e-10**
 

mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.instr

@@ -47,16 +47,16 @@
 * instrument is compiled, if they do not exist yet, by
 * generate_dispersion_files.py (Python with numpy).
 *
-* %Example: h=1.2 l=1 dE=5 sample_choice=0 Detector: detector_I=1.46e-11
-* %Example: h=1.2 l=1 dE=5 sample_choice=1 Detector: detector_I=1.41e-11
-* %Example: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.46e-11
+* %Example: h=1.2 l=1 dE=5 sample_choice=0 Detector: detector_I=1.85e-10
+* %Example: h=1.2 l=1 dE=5 sample_choice=1 Detector: detector_I=1.79e-10
+* %Example: h=1.2 l=1 dE=5 sample_choice=0 primitive_cell=1 Detector: detector_I=1.85e-10
 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=2 Detector: detector_I=8.91e-10
 * %Example: h=1.15 l=0 dE=3 T=2 sample_choice=3 Detector: detector_I=9.41e-10
 * %Scan: mcrun Test_Dispersion_relation_TAS.instr -M -N 5,5 -n 2e5 h=0.8,1.2 l=1.1,1.3 dE=5 sample_choice=0 Detector: detector_I={
-*   8.964e-12,1.086e-11,8.771e-12,1.261e-11,1.151e-11,7.748e-12,1.057e-11,1.137e-11,
-*   1.199e-11,1.014e-11,6.126e-12,8.926e-12,1.006e-11,1.007e-11,8.424e-12,3.55e-12,
-*   7.184e-12,8.291e-12,8.243e-12,5.35e-12,9.55e-13,4.109e-12,5.612e-12,4.924e-12,
-*   1.74e-12 }
+*   1.136e-10,1.377e-10,1.112e-10,1.598e-10,1.459e-10,9.821e-11,1.34e-10,1.441e-10,
+*   1.52e-10,1.285e-10,7.765e-11,1.131e-10,1.275e-10,1.276e-10,1.068e-10,4.5e-11,
+*   9.106e-11,1.051e-10,1.045e-10,6.781e-11,1.21e-11,5.208e-11,7.113e-11,6.241e-11,
+*   2.205e-11 }
 *
 * %Parameters
 * h: [r.l.u.]       Scattering vector component along a* (sample x axis)

mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/Test_Dispersion_relation_TAS.md

@@ -75,7 +75,7 @@ of the two scanned parameters.
 
 `generate_phonon_dispersion.py` writes the `Phonon_simple` dispersion of Pb to
 a grid. The last column is the `Phonon_simple` cross section,
-b²/(4π) · (ħ²κ²/2m_n)/ω / M per unit cell, in the form `Dispersion_relation`
+b² · (ħ²κ²/2M)/ω per unit cell, M in amu, in the form `Dispersion_relation`
 reads (see the component header). It depends on the full scattering vector κ,
 while the file is read at q folded into one unit cell, so the file is written
 for the Brillouin zone around one reciprocal lattice point (`--zone`, default
@@ -95,7 +95,7 @@ broadened by the resolution of the spectrometer:
 
 | | `Dispersion_relation` | `Phonon_simple` |
 |---|---|---|
-| Total intensity over the scan (n/s) | 1.572e-9 | 1.585e-9 |
+| Total intensity over the scan (n/s) | 1.993e-8 | 2.009e-8 |
 | Ring radius from intensity (r.l.u.), contour 0.215 | 0.215 | 0.215 |
 
 The intensity maps have a correlation of 0.999, and the total intensities

mcstas-comps/examples/Tests_samples/Test_Dispersion_relation_TAS/generate_phonon_dispersion.py

@@ -10,7 +10,7 @@
 
 The last column is the Phonon_simple cross section written in the form that
 Dispersion_relation reads, in barn/sr per unit cell of side a/2:
-    I = (rho_atom/rho_cell) * b^2/(4 pi) * (hbar^2 kappa^2/2m_n)/omega / M.
+    I = (rho_atom/rho_cell) * b^2 * (hbar^2 kappa^2/2M)/omega, M in amu.
 It depends on the full scattering vector kappa, which the folded grid cannot
 hold, so it is written for the Brillouin zone around one reciprocal lattice
 point given by --zone: each folded coordinate is unfolded to lie within half
@@ -32,6 +32,7 @@
 import numpy as np
 
 HBAR2_OVER_2MN = 2.0721  # hbar^2/(2 m_n) [meV AA^2], K2V^2*VS2E in McStas
+NEUTRON_MASS_U = 1.00866  # m_n/u, so that M is in amu
 
 
 def phonon_simple_omega(qx, qy, qz, a, c):
@@ -46,15 +47,15 @@ def phonon_simple_omega(qx, qy, qz, a, c):
 def phonon_simple_cross_section(kappa2, omega, a, b, M, cell_volume=None):
     """Phonon_simple cross section in barn/sr per unit cell of the given volume.
 
-    Phonon_simple weights a ray by rho_atom/(4 pi) * (hbar^2 kappa^2/2m_n)/omega * b^2/M,
+    Phonon_simple weights a ray by rho_atom * (hbar^2 kappa^2/2m_n)/omega * b^2 * (m_n/u)/M,
     with b in fm, while Dispersion_relation uses rho_cell * I with I in barn/sr.
     The default cell is the cube of side a/2.
     """
     rho_atom = 4/a**3
     rho_cell = 1/(cell_volume if cell_volume is not None else (a/2)**3)
     b2_barn = b**2*0.01
     with np.errstate(divide="ignore", invalid="ignore"):
-        cross_section = (rho_atom/rho_cell)*b2_barn/(4*np.pi)*HBAR2_OVER_2MN*kappa2/omega/M
+        cross_section = (rho_atom/rho_cell)*b2_barn*HBAR2_OVER_2MN*kappa2/omega*NEUTRON_MASS_U/M
     # omega = 0 only at the reciprocal lattice points, where no phonon is excited
     return np.where(omega > 0, cross_section, 0)
 
@@ -88,7 +89,7 @@ def main():
     parser.add_argument("--a", type=float, default=4.95, help="fcc lattice constant [AA] (Pb)")
     parser.add_argument("--c", type=float, default=10.0, help="Velocity of sound [meV AA] (Pb)")
     parser.add_argument("--b", type=float, default=9.405, help="Scattering length [fm] (Pb)")
-    parser.add_argument("--M", type=float, default=207.2, help="Atomic mass [a.u.] (Pb)")
+    parser.add_argument("--M", type=float, default=207.2, help="Atomic mass [amu] (Pb)")
     parser.add_argument("--zone", type=float, nargs=3, default=[1, 0, 1], metavar=("H", "K", "L"),
                         help="Reciprocal lattice point, in units of 4 pi/a, whose zone the cross section is written for")
     parser.add_argument("--cell", choices=["cubic", "primitive"], default="cubic",

mcstas-comps/examples/Tests_union/Test_inhomogenous_process/README.md

@@ -54,7 +54,7 @@ gaussian function and the heavyside function.
 - **Test: sample=inc_linear Detector: lin_det_I=0.017457**
 - **Test: sample=inho_linear Detector: lin_det_I=0.017457**
 - **Test: sample=inc_and_inho Detector: lin_det_I=0.017457**
-- **Test: sample=phonon thick=0.001 Detector: lin_det_I=0.00505839**
+- **Test: sample=phonon thick=0.001 Detector: lin_det_I=0.00529178**
 - **Test: sample=inc_trans Detector: lin_det_I=0.081029**
 - **Test: sample=inho_trans Detector: lin_det_I=0.081040**
 - **Test: sample=gauss_heavy Detector: lin_det_I=0.0985658**

mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.instr

@@ -53,7 +53,7 @@
 * %Example: sample=inc_linear Detector: lin_det_I=0.017457
 * %Example: sample=inho_linear Detector: lin_det_I=0.017457
 * %Example: sample=inc_and_inho Detector: lin_det_I=0.017457
-* %Example: sample=phonon thick=0.001 Detector: lin_det_I=0.00505839
+* %Example: sample=phonon thick=0.001 Detector: lin_det_I=0.00529178
 * %Example: sample=inc_trans Detector: lin_det_I=0.081029
 * %Example: sample=inho_trans Detector: lin_det_I=0.081040
 * %Example: sample=gauss_heavy Detector: lin_det_I=0.0985658

mcstas-comps/samples/Phonon_simple.comp

@@ -37,7 +37,7 @@
 * sigma_inc: [barns]  Incoherent scattering cross section per atom
 * a: [AA]             fcc Lattice constant
 * b: [fm]             Scattering length
-* M: [a.u.]           Atomic mass
+* M: [amu]            Atomic mass
 * c: [meV/AA^(-1)]    Velocity of sound
 * DW: [1]             Debye-Waller factor
 * T: [K]              Temperature
@@ -484,10 +484,10 @@ TRACE
     my_a_f = V_my_a_v / v_f;
     bose_factor = nbose (omega, T);
     p1 = exp (-(V_my_s * (l_i + l_o) + my_a_i * l_i + my_a_f * l_o));                 /* Absorption factor */
-    p2 = nf * solid_angle * l_full * V_rho / (4 * PI);                                /* Focusing factors; assume random choice of n_f possibilities */
+    p2 = nf * solid_angle * l_full * V_rho;                                           /* Focusing factors; assume random choice of n_f possibilities */
     p3 = (v_f / v_i) * DW * (kappa2 * K2V * K2V * VS2E) / fabs (omega) * bose_factor; /* Cross section factor 1 */
     p4 = 2 * VS2E * v_f / J_factor;                                                   /* Jacobian of delta functions in cross section */
-    p5 = b * b / M;                                                                   /* Cross section factor 2 */
+    p5 = b * b * 1.00866 / M;                                                         /* Cross section factor 2; 1.00866 = m_n/u, M in amu */
     p *= p1 * p2 * p3 * p4 * p5;
   } /* else transmit: Neutron did not hit the sample */
   free (vf_list);

mcstas-comps/union/PhononSimple_process.comp

@@ -39,7 +39,7 @@
 * unit_cell_volume:  [AA^3]   Unit cell volume
 * a:                 [AA]     fcc lattice constant
 * c:                 [meV*AA] Velocity of sound
-* M:                 [units]  Nucleus atomic mass in units
+* M:                 [amu]    Atomic mass
 * b:                 [fm]     Scattring length
 * T:                 [K]      Temperature
 * longitudinal:      [0/1]    Simulate longitudinal branches
@@ -505,7 +505,7 @@ INITIALIZE
     exit (1);
   }
 
-  PhononSimple_storage.my_constant = packing_factor * 4.0 / (a * a * a) * DW * K2V * K2V * VS2E * 2 * VS2E * K2V * b * b / M;
+  PhononSimple_storage.my_constant = packing_factor * 4.0 / (a * a * a) * DW * K2V * K2V * VS2E * 2 * VS2E * K2V * b * b * 4 * PI * 1.00866 / M; /* 4 PI: mu = 4 PI rho dsigma/dOmega; 1.00866 = m_n/u, M in amu */
   // my_constant = V_rho*DW*K2V*K2V*VS2E*2*VS"E*K2V*b*b/M // done in initialize
 
   // First initialise This_process with default values: