Add MCViNE sample kernels as Union processes (and standalone samples)

mcstas-comps/contrib/MCViNE_Broadened_E_Q.comp

@@ -0,0 +1,149 @@
+/*******************************************************************************
+*
+* McStas, neutron ray-tracing package
+*
+* Component: MCViNE_Broadened_E_Q
+*
+* %I
+* Written by: Fahima Islam (McStas port of the MCViNE Broadened_E_Q_Kernel; MCViNE by J. Y. Y. Lin et al.)
+* Date: 2026-09-29
+* Origin: MCViNE, https://github.com/mcvine/mcvine (mccomponents/lib/kernels/sample/Broadened_E_Q_Kernel.icc)
+*
+* Isotropic dispersion E(Q) broadened by a Gaussian of Q-dependent width.
+*
+* %D
+* S(Q,E) = S(Q) G(E - E(Q); sigma(Q)) with G a normalised Gaussian of standard
+* deviation sigma(Q).
+* Q is sampled uniformly in [Qmin,Qmax]; the energy offset is sampled from the
+* line shape (up to 100 attempts). Neutrons with Ei below min(E(Q)-3w(Q)) do not
+* scatter.
+* 
+* Functions are given as strings evaluated at run time (MCViNE uses fparser):
+* + - * / ^ (or **), % , comparisons, && ||, if(c,a,b), sin cos tan asin acos
+* atan sinh cosh tanh exp log log10 log2 sqrt abs floor ceil int sign cbrt,
+* pow atan2 min max hypot fmod, constants pi and e.
+*
+* <b>Transport</b> (same as MCViNE HomogeneousNeutronScatterer): on the first
+* passage the neutron is forced to scatter at a uniformly chosen depth, weighted
+* by path length, attenuation exp(-(mu+sigma)x) and the scattering coefficient;
+* on exit it is attenuated by exp(-(mu+sigma)L_out). With order>1 further
+* scatterings are sampled analogically (truncated exponential). With
+* p_transmit>0 that fraction of events is kept as the attenuated direct beam.
+* mu(v) = absorption_coefficient*2200/v. Events for which the kernel cannot
+* scatter (kinematically forbidden) are absorbed.
+*
+* Geometry: box (xwidth,yheight,zdepth), cylinder (radius,yheight), hollow
+* cylinder (+thickness) or sphere (radius only). All vectors (Q, targets,
+* reciprocal vectors, atom positions) are in the component's local frame.
+*
+* Uses the shared runtime share/mcvine-lib.h/.c.
+*
+* Example: MCViNE_Broadened_E_Q(E_Q="20*sin(Q*1.5)^2", S_Q="1", sigma_Q="0.5", Qmin=0, Qmax=10, scattering_coefficient=10, radius=0.01, yheight=0.05)
+*
+* %P
+* INPUT PARAMETERS:
+* E_Q: [str] E(Q) expression [meV], variable Q [AA^-1]
+* S_Q: [str] S(Q) expression
+* sigma_Q: [str] Gaussian sigma(Q) [meV]
+* Qmin: [AA^-1] Lower Q bound
+* Qmax: [AA^-1] Upper Q bound
+* unbiased: [1] 0: MCViNE sampling (retries without acceptance correction: over-weights when part of the Q range is forbidden). 1: single attempt, unbiased
+* absorption_coefficient: [m^-1] Absorption coefficient at 2200 m/s (MCViNE convention; scales as 1/v)
+* scattering_coefficient: [m^-1] Scattering coefficient
+* sigma_abs: [barn] Alternative: absorption cross section per unit cell at 2200 m/s (used when Vc>0)
+* sigma_scat: [barn] Alternative: scattering cross section per unit cell (used when Vc>0)
+* Vc: [AA^3] Unit cell volume; when >0, coefficients are computed from sigma_abs, sigma_scat
+* radius: [m] Outer radius of a cylinder (with yheight) or of a sphere (yheight=0)
+* xwidth: [m] Width of a box sample
+* yheight: [m] Height of a box or cylinder sample
+* zdepth: [m] Depth of a box sample
+* thickness: [m] Wall thickness of a hollow cylinder (0: filled)
+* pack: [1] Packing factor (scales absorption and scattering coefficients)
+* p_transmit: [1] Monte Carlo fraction of events kept as unscattered transmitted beam (0: always scatter)
+* order: [1] Maximum number of scattering events per neutron (1: single scattering)
+*
+* %L
+* MCViNE documentation: https://mcvine.github.io
+*
+* %E
+*******************************************************************************/
+
+DEFINE COMPONENT MCViNE_Broadened_E_Q
+
+SETTING PARAMETERS (string E_Q="10", string S_Q="1", string sigma_Q="1", Qmin=0, Qmax=10, int unbiased=0, absorption_coefficient=0, scattering_coefficient=0, sigma_abs=0, sigma_scat=0, Vc=0, radius=0, xwidth=0, yheight=0, zdepth=0, thickness=0, pack=1, p_transmit=0, int order=1)
+
+NOACC
+
+SHARE
+%{
+%include "read_table-lib"
+%include "mcvine-lib"
+%}
+
+DECLARE
+%{
+  mcvine_kernel_Broadened_E_Q kernel;
+  mcvine_scatterer scatterer;
+%}
+
+INITIALIZE
+%{
+  mcvine_shape shape;
+  double mu = 0, sig = 0;
+  memset (&kernel, 0, sizeof (kernel));
+  const char* vars[] = { "Q" };
+  if (mcvine_func_setup (&kernel.m_E_Q, E_Q, NULL, 0, 1, vars, NAME_CURRENT_COMP))
+    exit (-1);
+  if (mcvine_func_setup (&kernel.m_S_Q, S_Q, NULL, 0, 1, vars, NAME_CURRENT_COMP))
+    exit (-1);
+  if (mcvine_func_setup (&kernel.m_W_Q, sigma_Q, NULL, 0, 1, vars, NAME_CURRENT_COMP))
+    exit (-1);
+  if (Qmin < 0 || Qmin >= Qmax) {
+    fprintf (stderr, "%s: need 0 <= Qmin < Qmax\n", NAME_CURRENT_COMP);
+    exit (-1);
+  }
+  kernel.m_Qmin = Qmin;
+  kernel.m_Qmax = Qmax;
+  kernel.m_lorentzian = 0;
+  kernel.m_unbiased = unbiased;
+  mcvine_Broadened_E_Q_init (&kernel);
+
+  if (Vc > 0) {
+    mu = mcvine_xs2coeff (sigma_abs, Vc);
+    sig = mcvine_xs2coeff (sigma_scat, Vc);
+  } else {
+    mu = absorption_coefficient;
+    sig = scattering_coefficient;
+  }
+  if (!(sig > 0)) {
+    fprintf (stderr, "%s: give scattering_coefficient>0 (or sigma_scat>0 and Vc>0)\n", NAME_CURRENT_COMP);
+    exit (-1);
+  }
+  if (mcvine_shape_init (&shape, radius, xwidth, yheight, zdepth, thickness, NAME_CURRENT_COMP))
+    exit (-1);
+  mcvine_scatterer_init (&scatterer, &shape, mu, sig, pack, p_transmit, order);
+%}
+
+TRACE
+%{
+  int ret = mcvine_scatterer_interact (&scatterer, &kernel, mcvine_S_Broadened_E_Q, _particle);
+  if (ret < 0)
+    ABSORB;
+  if (ret > 0)
+    SCATTER;
+%}
+
+MCDISPLAY
+%{
+  if (scatterer.shape.type == MCVINE_SHAPE_BOX) {
+    box (0, 0, 0, xwidth, yheight, zdepth, 0, 0, 1, 0);
+  } else if (scatterer.shape.type == MCVINE_SHAPE_CYLINDER) {
+    cylinder (0, 0, 0, radius, yheight, 0, 0, 1, 0);
+    if (thickness > 0)
+      cylinder (0, 0, 0, radius - thickness, yheight, 0, 0, 1, 0);
+  } else {
+    sphere (0, 0, 0, radius);
+  }
+%}
+
+END

mcstas-comps/contrib/MCViNE_Broadened_E_Q_process.comp

@@ -0,0 +1,131 @@
+/*******************************************************************************
+*
+* McStas, neutron ray-tracing package
+*
+* Component: MCViNE_Broadened_E_Q_process
+*
+* %I
+* Written by: Fahima Islam (McStas port of the MCViNE Broadened_E_Q_Kernel; MCViNE by J. Y. Y. Lin et al.)
+* Date: 2026-09-29
+* Origin: MCViNE, https://github.com/mcvine/mcvine (mccomponents/lib/kernels/sample/Broadened_E_Q_Kernel.icc)
+*
+* Union process: Isotropic dispersion E(Q) broadened by a Gaussian of Q-dependent width.
+*
+* %D
+* S(Q,E) = S(Q) G(E - E(Q); sigma(Q)) with G a normalised Gaussian of standard
+* deviation sigma(Q).
+* Q is sampled uniformly in [Qmin,Qmax]; the energy offset is sampled from the
+* line shape (up to 100 attempts). Neutrons with Ei below min(E(Q)-3w(Q)) do not
+* scatter.
+* 
+* Functions are given as strings evaluated at run time (MCViNE uses fparser):
+* + - * / ^ (or **), % , comparisons, && ||, if(c,a,b), sin cos tan asin acos
+* atan sinh cosh tanh exp log log10 log2 sqrt abs floor ceil int sign cbrt,
+* pow atan2 min max hypot fmod, constants pi and e.
+*
+* <b>Union process.</b> Part of the Union components: define this process, collect it
+* into a material with Union_make_material (absorption is set there with
+* my_absorption, the absorption inverse penetration depth at 2200 m/s), assign
+* the material to Union_box/Union_cylinder/Union_sphere/Union_mesh geometries,
+* and add a Union_master after them. Geometry, attenuation and multiple
+* scattering are handled by Union_master; this component provides the MCViNE
+* kernel: the scattering coefficient and the final-state sampling (weight).
+* Isotropic process (powder/liquid-like); rotation is irrelevant.
+* The kernel code is shared with the standalone component MCViNE_Broadened_E_Q (mcvine-lib.c).
+* Uses share/mcvine-lib.h/.c and share/mcvine-union-lib.h/.c; the Union core
+* registers the process type MCViNE (share/union-lib.c, share/union-suffix.c).
+*
+* Example: MCViNE_Broadened_E_Q_process(E_Q="20*sin(Q*1.5)^2", S_Q="1", sigma_Q="0.5", Qmin=0, Qmax=10, scattering_coefficient=10)
+*
+* %P
+* INPUT PARAMETERS:
+* E_Q: [str] E(Q) expression [meV], variable Q [AA^-1]
+* S_Q: [str] S(Q) expression
+* sigma_Q: [str] Gaussian sigma(Q) [meV]
+* Qmin: [AA^-1] Lower Q bound
+* Qmax: [AA^-1] Upper Q bound
+* unbiased: [1] 0: MCViNE sampling (retries without acceptance correction: over-weights when part of the Q range is forbidden). 1: single attempt, unbiased
+* scattering_coefficient: [m^-1] Scattering coefficient (inverse penetration depth for scattering)
+* sigma_scat: [barn] Alternative: scattering cross section per unit cell (used when Vc>0)
+* Vc: [AA^3] Unit cell volume; when >0 the coefficient is sigma_scat/Vc
+* packing_factor: [1] Packing factor (scales the scattering coefficient)
+* interact_fraction: [1] Union: fraction of interactions forced to this process (-1: by cross section)
+*
+* %L
+* MCViNE documentation: https://mcvine.github.io
+*
+* %E
+*******************************************************************************/
+
+DEFINE COMPONENT MCViNE_Broadened_E_Q_process
+
+SETTING PARAMETERS (string E_Q="10", string S_Q="1", string sigma_Q="1", Qmin=0, Qmax=10, int unbiased=0, scattering_coefficient=0, sigma_scat=0, Vc=0, packing_factor=1, interact_fraction=-1)
+
+NOACC
+
+SHARE
+%{
+  #ifndef Union
+  #error "ERROR: MCViNE_Broadened_E_Q_process requires the Union library. Add a Union_master component to load the libraries."
+  #endif
+  %include "read_table-lib"
+  %include "mcvine-lib"
+  %include "mcvine-union-lib"
+  #ifndef PROCESS_DETECTOR
+  #define PROCESS_DETECTOR dummy
+  #endif
+  #ifndef PROCESS_MCVINE_DETECTOR
+  #define PROCESS_MCVINE_DETECTOR dummy
+  #endif
+%}
+
+DECLARE
+%{
+  mcvine_kernel_Broadened_E_Q kernel;
+  struct MCViNE_physics_storage_struct storage;
+  struct scattering_process_struct This_process;
+  struct global_process_element_struct global_process_element;
+%}
+
+INITIALIZE
+%{
+  double mu = 0, sig = 0;
+  memset (&kernel, 0, sizeof (kernel));
+  const char* vars[] = { "Q" };
+  if (mcvine_func_setup (&kernel.m_E_Q, E_Q, NULL, 0, 1, vars, NAME_CURRENT_COMP))
+    exit (-1);
+  if (mcvine_func_setup (&kernel.m_S_Q, S_Q, NULL, 0, 1, vars, NAME_CURRENT_COMP))
+    exit (-1);
+  if (mcvine_func_setup (&kernel.m_W_Q, sigma_Q, NULL, 0, 1, vars, NAME_CURRENT_COMP))
+    exit (-1);
+  if (Qmin < 0 || Qmin >= Qmax) {
+    fprintf (stderr, "%s: need 0 <= Qmin < Qmax\n", NAME_CURRENT_COMP);
+    exit (-1);
+  }
+  kernel.m_Qmin = Qmin;
+  kernel.m_Qmax = Qmax;
+  kernel.m_lorentzian = 0;
+  kernel.m_unbiased = unbiased;
+  mcvine_Broadened_E_Q_init (&kernel);
+
+  if (Vc > 0)
+    sig = mcvine_xs2coeff (sigma_scat, Vc);
+  else
+    sig = scattering_coefficient;
+  if (!(sig > 0)) {
+    fprintf (stderr, "%s: give scattering_coefficient>0 (or sigma_scat>0 and Vc>0)\n", NAME_CURRENT_COMP);
+    exit (-1);
+  }
+  storage.m_kernel = &kernel;
+  storage.m_S = mcvine_S_Broadened_E_Q;
+  storage.m_my_scattering = sig * packing_factor;
+  storage.m_kind = MCVINE_UNION_GENERIC;
+  mcvine_union_register (&This_process, &global_process_element, &storage, NAME_CURRENT_COMP, INDEX_CURRENT_COMP, interact_fraction, 0, ROT_A_CURRENT_COMP);
+%}
+
+TRACE
+%{
+  // the simulation is done in Union_master
+%}
+
+END