diff --git a/.github/workflows/mcstas-basictest.yml b/.github/workflows/mcstas-basictest.yml index 3a11718e3f..3b96230e86 100644 --- a/.github/workflows/mcstas-basictest.yml +++ b/.github/workflows/mcstas-basictest.yml @@ -435,7 +435,7 @@ jobs: compindex=0 if [ "$NUMCHANGEDCOMPS" != "0" ]; then - if [ "$NUMCHANGEDCOMPS" -lt "5" ]; then + if [ "$NUMCHANGEDCOMPS" -lt "11" ]; then for comp in $CHANGEDCOMPS; do echo Finding tests including component $comp diff --git a/mcstas-comps/contrib/He3TubePack_detector.comp b/mcstas-comps/contrib/He3TubePack_detector.comp new file mode 100644 index 0000000000..f086cfc233 --- /dev/null +++ b/mcstas-comps/contrib/He3TubePack_detector.comp @@ -0,0 +1,253 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE detector kernels by Jiao Lin +* Date: 24.09.2026 +* Origin: DTU Physics +* +* Pack of 3He position sensitive detector tubes with absorption efficiency and MCViNE-compatible event mode output +* +* %D +* Port of the MCViNE detector kernels mccomponents/lib/kernels/detector/He3, He3Tube, +* Z2Channel, Tof2Channel and EventModeMCA (J.Y.Y. Lin et al., Caltech/ORNL), as used for +* the tube detector packs of the SNS direct geometry spectrometers in MCViNE. +* +* A pack of ntubes parallel, vertical (along y) cylindrical tubes of radius tube_radius +* and length tube_length, placed side by side along x with centre distance tube_spacing, +* centred on the component position. Each tube is divided into npixels pixels along its +* axis. For each tube crossed by the neutron, the absorption probability in the 3He gas is +* P = 1 - exp(-mu L), mu = 5333 barn * lambda/1.798 AA * n, n = pressure/(kB 300 K) +* with L the path length in the tube. The absorption point is sampled along the path and +* gives the pixel; the time of flight gives the tof channel. The ray then continues with +* the transmitted weight (tube walls are not modelled). +* +* Output: +* - a 2D monitor of intensity vs. tube and pixel, and a 1D time of flight monitor, +* - optionally (events_file) binary events in the MCViNE EventModeMCA format: records of +* {unsigned int pixelID; unsigned int tofChannel; double n} (16 bytes, native byte +* order), with pixelID = (pack_index*ntubes + tube)*npixels + pixel and +* tofChannel = floor((t - tofmin)/tofstep). Events outside the tof range are not written. +* With MPI, each node writes its own file with the suffix .. +* +* Example: He3TubePack_detector(ntubes=8, npixels=128, tube_radius=0.0127, tube_length=1, +* tube_spacing=0.0286, pressure=10, tofmin=0, tofmax=0.02, tofstep=1e-5, +* events_file="events.dat") +* +* %P +* INPUT PARAMETERS: +* ntubes: [1] Number of tubes in the pack +* npixels: [1] Number of pixels per tube +* tube_radius: [m] Inner radius of the tubes +* tube_length: [m] Length of the tubes (active length) +* tube_spacing: [m] Distance between tube axes. 0: 2*tube_radius +* pressure: [atm] 3He pressure +* tofmin: [s] Minimum time of flight of the tof channels +* tofmax: [s] Maximum time of flight of the tof channels +* tofstep: [s] Width of the tof channels +* pack_index: [1] Index of the pack, offsets the pixel IDs (several packs in an instrument) +* events_file: [string] Name of the binary event file (MCViNE format). Empty: no event file +* filename: [string] Base name of the histogram output files (default: component name) +* restore_neutron: [1] If set, the monitor does not influence the neutron state +* nowritefile: [1] If set, monitor will skip writing histograms to disk +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* +* %E +******************************************************************************/ + +DEFINE COMPONENT He3TubePack_detector + +SETTING PARAMETERS(int ntubes=8, int npixels=128, tube_radius=0.0127, tube_length=1, tube_spacing=0, pressure=10, + tofmin=0, tofmax=0.02, tofstep=1e-5, int pack_index=0, string events_file="", string filename="", + int restore_neutron=0, int nowritefile=0) + +NOACC + +SHARE +%{ + // one event in the MCViNE EventModeMCA binary format + struct he3tube_event { + unsigned int pixelID; + unsigned int tofChannelNo; + double n; + }; +%} + +DECLARE +%{ + double mu_ref; // absorption coefficient at 1.798 AA [1/m] + double pitch; + int ntof; + DArray2d Npix; + DArray2d Ppix; + DArray2d P2pix; + DArray1d Ntof; + DArray1d Ptof; + DArray1d P2tof; + FILE* evfp; + long nevents; +%} + +INITIALIZE +%{ + if (ntubes < 1 || npixels < 1 || tube_radius <= 0 || tube_length <= 0 || pressure <= 0 || tofmax <= tofmin || tofstep <= 0) { + fprintf (stderr, "He3TubePack_detector:%s: ERROR: invalid parameters\n", NAME_CURRENT_COMP); + exit (-1); + } + pitch = tube_spacing > 0 ? tube_spacing : 2 * tube_radius; + if (pitch < 2 * tube_radius) { + fprintf (stderr, "He3TubePack_detector:%s: ERROR: tube_spacing is smaller than the tube diameter\n", NAME_CURRENT_COMP); + exit (-1); + } + // MCViNE He3 kernel: sigma_abs = 5333 barn at 1.798 AA, n = P[Pa] * 2.414e20 m^-3 (300 K) + mu_ref = 5333.0e-28 * pressure * 101325.0 * 2.414e20; + ntof = (int)ceil ((tofmax - tofmin) / tofstep); + Npix = create_darr2d (ntubes, npixels); + Ppix = create_darr2d (ntubes, npixels); + P2pix = create_darr2d (ntubes, npixels); + Ntof = create_darr1d (ntof); + Ptof = create_darr1d (ntof); + P2tof = create_darr1d (ntof); + if (!strcmp (filename, "\0")) + sprintf (filename, "%s", NAME_CURRENT_COMP); + + evfp = NULL; + nevents = 0; + if (strlen (events_file) && strcmp (events_file, "NULL") && strcmp (events_file, "0")) { + char name[1024]; + #ifdef USE_MPI + if (mpi_node_count > 1) + snprintf (name, 1024, "%.1000s.%d", events_file, mpi_node_rank); + else + #endif + snprintf (name, 1024, "%.1000s", events_file); + char* path = mcfull_file (name, NULL); + evfp = fopen (path, "wb"); + if (!evfp) { + fprintf (stderr, "He3TubePack_detector:%s: ERROR: cannot open event file %s\n", NAME_CURRENT_COMP, path); + exit (-1); + } + free (path); + } +%} + +TRACE +%{ + // tubes crossed by the ray, sorted by entry time + int i, j, n_hit = 0, tube[64], order[64]; + double tin[64], tout[64], v = sqrt (vx * vx + vy * vy + vz * vz); + double x0 = -(ntubes - 1) * pitch / 2; + double p_det = 0; + if (v > 0) { + for (i = 0; i < ntubes && n_hit < 64; i++) { + double t0, t1; + if (cylinder_intersect (&t0, &t1, x - (x0 + i * pitch), y, z, vx, vy, vz, tube_radius, tube_length) && t1 > 0) { + if (t0 < 0) + t0 = 0; + tube[n_hit] = i; + tin[n_hit] = t0; + tout[n_hit] = t1; + order[n_hit] = n_hit; + n_hit++; + } + } + // insertion sort on entry time + for (i = 1; i < n_hit; i++) + for (j = i; j > 0 && tin[order[j]] < tin[order[j - 1]]; j--) { + int tmp = order[j]; + order[j] = order[j - 1]; + order[j - 1] = tmp; + } + double lambda = 2 * PI / (V2K * v), mu = mu_ref * lambda / 1.798, pw = p; + for (i = 0; i < n_hit; i++) { + int h = order[i]; + double L = (tout[h] - tin[h]) * v; + double pabs = 1 - exp (-mu * L); + if (pabs <= 0) + continue; + // absorption point: exponential distribution truncated to the path in the tube + double s = -log (1 - rand01 () * pabs) / mu; + double ta = tin[h] + s / v; + double ya = y + vy * ta, ta_abs = t + ta; + double w = pw * pabs; + int pix = (int)floor ((ya + tube_length / 2) / tube_length * npixels); + if (pix >= 0 && pix < npixels) { + double w2 = w * w; + Npix[tube[h]][pix] += 1; + Ppix[tube[h]][pix] += w; + P2pix[tube[h]][pix] += w2; + if (ta_abs >= tofmin && ta_abs < tofmax) { + int it = (int)floor ((ta_abs - tofmin) / tofstep); + if (it >= 0 && it < ntof) { + Ntof[it] += 1; + Ptof[it] += w; + P2tof[it] += w2; + if (evfp) { + struct he3tube_event ev; + ev.pixelID = (unsigned int)(((long)pack_index * ntubes + tube[h]) * npixels + pix); + ev.tofChannelNo = (unsigned int)it; + ev.n = w; + fwrite (&ev, sizeof (ev), 1, evfp); + nevents++; + } + } + } + } + p_det += w; + pw *= 1 - pabs; + } + if (n_hit > 0) { + if (!restore_neutron) { + // continue with the transmitted weight + p = pw; + SCATTER; + } + } + } + if (restore_neutron) { + RESTORE_NEUTRON (INDEX_CURRENT_COMP, x, y, z, vx, vy, vz, t, sx, sy, sz, p); + } +%} + +SAVE +%{ + if (!nowritefile) { + char fn[1024]; + snprintf (fn, 1024, "%.1000s_pixels", filename); + DETECTOR_OUT_2D ("He3 tube pack", "Tube index", "Pixel index", -0.5, ntubes - 0.5, -0.5, npixels - 0.5, ntubes, npixels, &Npix[0][0], &Ppix[0][0], + &P2pix[0][0], fn); + snprintf (fn, 1024, "%.1000s_tof", filename); + DETECTOR_OUT_1D ("He3 tube pack time of flight", "Time of flight [s]", "Intensity", "t", tofmin, tofmin + ntof * tofstep, ntof, &Ntof[0], &Ptof[0], + &P2tof[0], fn); + } + if (evfp) + fflush (evfp); +%} + +FINALLY +%{ + if (evfp) { + fclose (evfp); + MPI_MASTER (printf ("He3TubePack_detector:%s: %ld events written to %s (MCViNE EventModeMCA format)\n", NAME_CURRENT_COMP, nevents, events_file);); + } + destroy_darr2d (Npix); + destroy_darr2d (Ppix); + destroy_darr2d (P2pix); + destroy_darr1d (Ntof); + destroy_darr1d (Ptof); + destroy_darr1d (P2tof); +%} + +MCDISPLAY +%{ + int i; + double x0 = -(ntubes - 1) * pitch / 2; + for (i = 0; i < ntubes; i++) + cylinder (x0 + i * pitch, 0, 0, tube_radius, tube_length, 0, 0, 1, 0); +%} + +END diff --git a/mcstas-comps/data/fccNi-phonons/DOS b/mcstas-comps/data/fccNi-phonons/DOS new file mode 100644 index 0000000000..0ce4981214 Binary files /dev/null and b/mcstas-comps/data/fccNi-phonons/DOS differ diff --git a/mcstas-comps/data/fccNi-phonons/Ni.xyz b/mcstas-comps/data/fccNi-phonons/Ni.xyz new file mode 100644 index 0000000000..911c9284fa --- /dev/null +++ b/mcstas-comps/data/fccNi-phonons/Ni.xyz @@ -0,0 +1,3 @@ +1 +1.76 1.76 0 1.76 0 1.76 0 -1.76 -1.76 +Ni 0 0 0 10.3 58.6934 diff --git a/mcstas-comps/data/fccNi-phonons/Omega2 b/mcstas-comps/data/fccNi-phonons/Omega2 new file mode 100644 index 0000000000..07ac86e4b1 Binary files /dev/null and b/mcstas-comps/data/fccNi-phonons/Omega2 differ diff --git a/mcstas-comps/data/fccNi-phonons/Polarizations b/mcstas-comps/data/fccNi-phonons/Polarizations new file mode 100644 index 0000000000..b833c6038a Binary files /dev/null and b/mcstas-comps/data/fccNi-phonons/Polarizations differ diff --git a/mcstas-comps/data/fccNi-phonons/Qgridinfo b/mcstas-comps/data/fccNi-phonons/Qgridinfo new file mode 100644 index 0000000000..fb7e0c6467 --- /dev/null +++ b/mcstas-comps/data/fccNi-phonons/Qgridinfo @@ -0,0 +1,7 @@ +import math +twopi = 2*math.pi +b = 0.28409090909090912*twopi +b1=(b, -b, b) +b2=(b, b, -b) +b3=(-b, b, b) +n1=n2=n3=21 diff --git a/mcstas-comps/data/fccNi-phonons/README b/mcstas-comps/data/fccNi-phonons/README new file mode 100644 index 0000000000..0237f62b88 --- /dev/null +++ b/mcstas-comps/data/fccNi-phonons/README @@ -0,0 +1,14 @@ +Phonon dispersion of fcc Ni in the MCViNE "IDF" format, for the Union processes +CoherentPhononSingleXtal_process and CoherentPhononPowder_process. + +Qgridinfo reciprocal cell spanned by the grid (primitive bcc reciprocal cell of fcc + Ni, a=3.52 AA) and number of grid points (21x21x21, both ends included) +Omega2 squared phonon angular frequencies [rad^2/s^2] +Polarizations phonon polarization vectors +DOS phonon density of states (frequency in THz), used for the Debye-Waller factor +Ni.xyz crystal structure: number of atoms, lattice vectors a1 a2 a3 [AA], + then "Symbol x y z b_coh[fm] mass[amu]" with fractional coordinates + +Source: MCViNE test data (mcvine/packages/mccomponents/tests/mccomponents/sample/ +phonon/kernels/CoherentInelastic_PolyXtal/phonon-dispersion-fccNi-primitive-reciprocal-unitcell), +https://github.com/mcvine/mcvine diff --git a/mcstas-comps/examples/Tests_monitors/He3TubePack_test/He3TubePack_test.instr b/mcstas-comps/examples/Tests_monitors/He3TubePack_test/He3TubePack_test.instr new file mode 100644 index 0000000000..460c163ad1 --- /dev/null +++ b/mcstas-comps/examples/Tests_monitors/He3TubePack_test/He3TubePack_test.instr @@ -0,0 +1,82 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: He3TubePack_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 24.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_monitors +* +* Test of He3TubePack_detector: an 8-pack of 3He tubes behind a vanadium sample +* +* %Description +* A pulsed white beam hits a vanadium cylinder (Incoherent) which scatters towards +* an 8-pack of 3He position sensitive tubes (He3TubePack_detector, port of the MCViNE +* He3Tube / EventModeMCA detector kernels) at distance L and angle tth. A reference +* PSD monitor in front of the pack records the incoming intensity, so that the +* wavelength dependent detection efficiency can be compared between the two. +* With events=1 the pack writes the MCViNE event mode file events.dat. +* +* %Example: pressure=10 Detector: pack_I=101.7 +* +* %Parameters +* lambda0: [AA] Centre wavelength of the source +* dlambda: [AA] Half width of the wavelength band +* L: [m] Sample-pack distance +* tth: [deg] Scattering angle of the pack +* pressure: [atm] 3He pressure of the tubes +* events: [1] If set, write the binary event file events.dat +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT He3TubePack_test(lambda0=2, dlambda=1, L=2, tth=60, pressure=10, int events=0) + +DECLARE +%{ + double pitch; + double pack_w; +%} + +INITIALIZE +%{ + pitch = 0.0286; + pack_w = 8 * pitch; +%} + +TRACE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) RELATIVE ABSOLUTE + +COMPONENT source = Source_simple(radius=0.01, dist=10, focus_xw=0.01, focus_yh=0.03, + lambda0=lambda0, dlambda=dlambda, flux=1e10) +AT (0, 0, 0) RELATIVE origin +EXTEND %{ + t = 1e-5 * randpm1 (); +%} + +COMPONENT sample_pos = Arm() +AT (0, 0, 10) RELATIVE origin + +COMPONENT vanadium = Incoherent(radius=0.005, yheight=0.03, sigma_abs=5.08, sigma_inc=5.08, Vc=13.827, + target_index=2, focus_xw=pack_w, focus_yh=1) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT pack_arm = Arm() +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, tth, 0) RELATIVE sample_pos + +COMPONENT pack_ref = PSD_monitor(xwidth=pack_w, yheight=1, nx=8, ny=128, filename="pack_ref", restore_neutron=1) +AT (0, 0, L - 0.02) RELATIVE pack_arm + +COMPONENT pack = He3TubePack_detector(ntubes=8, npixels=128, tube_radius=0.0127, tube_length=1, tube_spacing=pitch, + pressure=pressure, tofmin=0.002, tofmax=0.01, tofstep=1e-5, + events_file=events ? "events.dat" : "NULL") +AT (0, 0, L) RELATIVE pack_arm + +END diff --git a/mcstas-comps/examples/Tests_union/CoherentPhononPowder_test/CoherentPhononPowder_test.instr b/mcstas-comps/examples/Tests_union/CoherentPhononPowder_test/CoherentPhononPowder_test.instr new file mode 100644 index 0000000000..f4d5d9d0e9 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/CoherentPhononPowder_test/CoherentPhononPowder_test.instr @@ -0,0 +1,97 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: CoherentPhononPowder_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 24.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of CoherentPhononPowder_process: coherent one-phonon scattering from fcc Ni powder +* +* %Description +* A monochromatic beam hits a Ni powder cylinder made from the Union process +* CoherentPhononPowder_process (port of the MCViNE CoherentInelastic_PolyXtal +* kernel), with the fcc Ni phonon dispersion from MCViNE (data/fccNi-phonons). +* +* method=0 samples phonon wave vectors as MCViNE does and scatters into 4pi. The +* banana monitor then shows scattering angle vs. final energy, i.e. the powder +* averaged one-phonon S(Q,E) seen through the kinematic constraints. +* method=1 picks the final direction with the Union focusing of the geometry; with +* focus=1 the scattering is focused onto the energy monitor at angle tth. +* +* %Example: method=1 focus=1 Detector: Edet_I=3.07e-10 +* +* %Parameters +* Ei: [meV] Incident energy +* T: [K] Sample temperature +* tth: [deg] Scattering angle of the energy monitor +* method: [1] Sampling method of CoherentPhononPowder_process, 0: MCViNE (4pi), 1: focusing +* focus: [1] With method=1: 1 focus scattering onto the energy monitor, 0: 4pi +* radius: [m] Radius of the Ni cylinder +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT CoherentPhononPowder_test(Ei=60, T=300, tth=60, int method=0, int focus=0, radius=0.002) + +DECLARE +%{ + int target; + double fw; +%} + +INITIALIZE +%{ + target = (focus && method == 1) ? 3 : 0; + fw = (focus && method == 1) ? 0.05 : 0; +%} + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.005, dist=1, focus_xw=2*radius, focus_yh=0.02, E0=Ei, dE=0, flux=1) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +COMPONENT Ni_phonons = CoherentPhononPowder_process(dispersion_dir="fccNi-phonons", xyz_file="fccNi-phonons/Ni.xyz", T=T, method=method) +AT (0, 0, 0) RELATIVE sample_pos + +// Ni: sigma_abs = 4.49 barn, V_uc = 10.9036 AA^3 +COMPONENT Ni = Union_make_material(my_absorption=100*4.49/10.9036, process_string="Ni_phonons") +AT (0, 0, 0) ABSOLUTE + +COMPONENT det_arm = Arm() +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, tth, 0) RELATIVE sample_pos + +COMPONENT sample = Union_cylinder(radius=radius, yheight=0.02, material_string="Ni", priority=1, + target_index=target, focus_xw=fw, focus_xh=fw) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Etheta = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, + options="banana, theta limits=[5 135] bins=130, energy limits=[0 100] bins=100", filename="Etheta.dat") +WHEN (!focus) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Edet = E_monitor(xwidth=0.05, yheight=0.05, nE=200, Emin=0, Emax=100, filename="Edet.dat", restore_neutron=1) +AT (0, 0, 1) RELATIVE det_arm + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/examples/Tests_union/CoherentPhononPowder_test/fccNi-phonons b/mcstas-comps/examples/Tests_union/CoherentPhononPowder_test/fccNi-phonons new file mode 120000 index 0000000000..f7c3c73a80 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/CoherentPhononPowder_test/fccNi-phonons @@ -0,0 +1 @@ +../../../data/fccNi-phonons/ \ No newline at end of file diff --git a/mcstas-comps/examples/Tests_union/CoherentPhononSingleXtal_test/CoherentPhononSingleXtal_test.instr b/mcstas-comps/examples/Tests_union/CoherentPhononSingleXtal_test/CoherentPhononSingleXtal_test.instr new file mode 100644 index 0000000000..93e6f518ac --- /dev/null +++ b/mcstas-comps/examples/Tests_union/CoherentPhononSingleXtal_test/CoherentPhononSingleXtal_test.instr @@ -0,0 +1,105 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: CoherentPhononSingleXtal_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 23.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of CoherentPhononSingleXtal_process: coherent one-phonon scattering from fcc Ni +* +* %Description +* A monochromatic beam hits a 2 mm thick fcc Ni single crystal plate, made from the +* Union process CoherentPhononSingleXtal_process (port of the MCViNE +* CoherentInelastic_SingleXtal kernel). The phonon dispersion (energies and +* polarizations on a 21x21x21 grid in the primitive reciprocal cell), the DOS +* used for the Debye-Waller factor and the crystal structure (Ni.xyz) are read +* from data/fccNi-phonons in the MCViNE IDF format (fcc Ni test data from MCViNE). +* +* With focus=0 the scattered neutrons go into 4pi and are recorded by a banana +* monitor as a map of scattering angle and final energy, where the phonon +* branches are seen as sharp lines. With focus=1 scattering is focused onto +* the small energy monitor at scattering angle tth, giving the energies where +* the dispersion fulfils energy and momentum conservation for that single final +* direction. +* +* The crystal can be rotated with A3 around the vertical axis. +* +* %Example: focus=1 Detector: Edet_I=3.78e-12 +* +* %Parameters +* Ei: [meV] Incident energy +* T: [K] Sample temperature +* A3: [deg] Rotation of the crystal around the vertical axis +* tth: [deg] Scattering angle of the energy monitor +* focus: [1] 0: scatter into 4pi, 1: focus scattering onto the energy monitor +* thick: [m] Thickness of the crystal plate +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT CoherentPhononSingleXtal_test(Ei=60, T=300, A3=0, tth=60, int focus=0, thick=0.002) + +DECLARE +%{ + int target; + double fw; +%} + +INITIALIZE +%{ + target = focus ? 3 : 0; + fw = focus ? 0.02 : 0; +%} + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.002, dist=1, focus_xw=0.01, focus_yh=0.01, E0=Ei, dE=0, flux=1) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +// Lattice, grid and polarization vectors are defined in the frame of the process component +COMPONENT Ni_phonons = CoherentPhononSingleXtal_process(dispersion_dir="fccNi-phonons", xyz_file="fccNi-phonons/Ni.xyz", T=T) +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, A3, 0) RELATIVE sample_pos + +// Ni: sigma_abs = 4.49 barn, V_uc = 10.9036 AA^3 +COMPONENT Ni = Union_make_material(my_absorption=100*4.49/10.9036, process_string="Ni_phonons") +AT (0, 0, 0) ABSOLUTE + +COMPONENT det_arm = Arm() +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, tth, 0) RELATIVE sample_pos + +COMPONENT crystal = Union_box(xwidth=0.02, yheight=0.02, zdepth=thick, material_string="Ni", priority=1, + target_index=target, focus_xw=fw, focus_xh=fw) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Etheta = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, + options="banana, theta limits=[5 135] bins=130, energy limits=[0 100] bins=200", filename="Etheta.dat") +WHEN (!focus) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Edet = E_monitor(xwidth=0.02, yheight=0.02, nE=500, Emin=0, Emax=100, filename="Edet.dat", restore_neutron=1) +AT (0, 0, 1) RELATIVE det_arm + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/examples/Tests_union/CoherentPhononSingleXtal_test/fccNi-phonons b/mcstas-comps/examples/Tests_union/CoherentPhononSingleXtal_test/fccNi-phonons new file mode 120000 index 0000000000..f7c3c73a80 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/CoherentPhononSingleXtal_test/fccNi-phonons @@ -0,0 +1 @@ +../../../data/fccNi-phonons/ \ No newline at end of file diff --git a/mcstas-comps/examples/Tests_union/DispersionPowder_test/DispersionPowder_test.instr b/mcstas-comps/examples/Tests_union/DispersionPowder_test/DispersionPowder_test.instr new file mode 100644 index 0000000000..82c55eff5f --- /dev/null +++ b/mcstas-comps/examples/Tests_union/DispersionPowder_test/DispersionPowder_test.instr @@ -0,0 +1,98 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: DispersionPowder_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 24.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of DispersionPowder_process: powder excitation with analytic, broadened dispersion +* +* %Description +* A monochromatic beam hits a powder sample made from the Union process +* DispersionPowder_process (port of the MCViNE E_Q_Kernel, Broadened_E_Q_Kernel and +* LorentzianBroadened_E_Q_Kernel). The example uses a roton-like dispersion +* E(Q) = sqrt((c*Q)^2 + (Q^2/2m)^2) - dip, with a minimum near Q0, +* written as an expression, with a Q dependent Gaussian or Lorentzian broadening and +* detailed balance at temperature T. +* method=0 samples Q as MCViNE (4pi), method=1 uses Union focusing (focus=1: onto the +* energy monitor at scattering angle tth). +* +* %Example: method=1 focus=1 Detector: Edet_I=0.130 +* +* %Parameters +* Ei: [meV] Incident energy +* tth: [deg] Scattering angle of the energy monitor +* broadening: [1] 0: none, 1: Gaussian, 2: Lorentzian +* width: [meV] Broadening at Q=0 (increases with Q) +* T: [K] Temperature +* method: [1] 0: MCViNE sampling (4pi), 1: focusing +* focus: [1] With method=1: 0 4pi, 1 focus onto the energy monitor +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT DispersionPowder_test(Ei=20, tth=40, int broadening=2, width=0.1, T=10, int method=0, int focus=0) + +DECLARE +%{ + int target; + double fw; +%} + +INITIALIZE +%{ + target = (focus && method == 1) ? 3 : 0; + fw = (focus && method == 1) ? 0.02 : 0; +%} + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.002, dist=1, focus_xw=0.01, focus_yh=0.01, E0=Ei, dE=0, flux=1e10) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +// phonon-roton like curve: linear at small Q, minimum near Q=2 AA^-1 +COMPONENT excitation = DispersionPowder_process(E_Q="abs(2.2*sin(1.3*Q))+0.3*Q+0.05", S_Q="1", W_Q="p1*(1+0.5*Q)", + broadening=broadening, sigma=1.34, unit_cell_volume=45, Qmin=0.05, Qmax=4, Emax=15, T=T, method=method, p1=width) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT powder_mat = Union_make_material(my_absorption=0, process_string="excitation") +AT (0, 0, 0) ABSOLUTE + +COMPONENT det_arm = Arm() +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, tth, 0) RELATIVE sample_pos + +COMPONENT sample = Union_cylinder(radius=0.005, yheight=0.02, material_string="powder_mat", priority=1, + target_index=target, focus_xw=fw, focus_xh=fw) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Etheta = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, + options="banana, theta limits=[5 135] bins=130, energy limits=[0 40] bins=160", filename="Etheta.dat") +WHEN (!focus) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Edet = E_monitor(xwidth=0.02, yheight=0.02, nE=200, Emin=0, Emax=40, filename="Edet.dat", restore_neutron=1) +AT (0, 0, 1) RELATIVE det_arm + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/examples/Tests_union/DispersionSingleXtal_test/DispersionSingleXtal_test.instr b/mcstas-comps/examples/Tests_union/DispersionSingleXtal_test/DispersionSingleXtal_test.instr new file mode 100644 index 0000000000..dc2c003591 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/DispersionSingleXtal_test/DispersionSingleXtal_test.instr @@ -0,0 +1,98 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: DispersionSingleXtal_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 24.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of DispersionSingleXtal_process: single crystal excitation with analytic dispersion +* +* %Description +* A monochromatic beam hits a single crystal plate made from the Union process +* DispersionSingleXtal_process (port of the MCViNE E_vQ_Kernel). The example uses a +* nearest neighbour spin wave in a simple cubic lattice (a=4 AA), +* E(Q) = D*(3 - cos(2 pi h) - cos(2 pi k) - cos(2 pi l)) + gap, +* with D=p1 and gap=p2, intensity S=1, and detailed balance at temperature T. +* The crystal is rotated by A3 around the vertical axis. +* With focus=0 a banana monitor records scattering angle vs. final energy (4pi), +* with focus=1 scattering is focused onto the energy monitor at scattering angle tth. +* +* %Example: focus=1 Detector: Edet_I=0.0914 +* +* %Parameters +* Ei: [meV] Incident energy +* A3: [deg] Crystal rotation around the vertical axis +* tth: [deg] Scattering angle of the energy monitor +* D: [meV] Spin wave stiffness parameter +* gap: [meV] Spin wave gap +* T: [K] Temperature +* focus: [1] 0: 4pi, 1: focus onto the energy monitor +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT DispersionSingleXtal_test(Ei=30, A3=0, tth=40, D=3, gap=1, T=50, int focus=0) + +DECLARE +%{ + int target; + double fw; +%} + +INITIALIZE +%{ + target = focus ? 3 : 0; + fw = focus ? 0.02 : 0; +%} + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.002, dist=1, focus_xw=0.01, focus_yh=0.01, E0=Ei, dE=0, flux=1e10) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +COMPONENT spinwave = DispersionSingleXtal_process(E_Q="p1*(3-cos(2*pi*h)-cos(2*pi*k)-cos(2*pi*l))+p2", S_Q="1", + sigma=5, unit_cell_volume=64, Emax=40, T=T, ax=4, by=4, cz=4, p1=D, p2=gap) +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, A3, 0) RELATIVE sample_pos + +COMPONENT crystal_mat = Union_make_material(my_absorption=0, process_string="spinwave") +AT (0, 0, 0) ABSOLUTE + +COMPONENT det_arm = Arm() +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, tth, 0) RELATIVE sample_pos + +COMPONENT crystal = Union_box(xwidth=0.02, yheight=0.02, zdepth=0.005, material_string="crystal_mat", priority=1, + target_index=target, focus_xw=fw, focus_xh=fw) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Etheta = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, + options="banana, theta limits=[5 135] bins=130, energy limits=[0 60] bins=120", filename="Etheta.dat") +WHEN (!focus) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Edet = E_monitor(xwidth=0.02, yheight=0.02, nE=300, Emin=0, Emax=60, filename="Edet.dat", restore_neutron=1) +AT (0, 0, 1) RELATIVE det_arm + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/examples/Tests_union/ElasticSQ_test/ElasticSQ_test.instr b/mcstas-comps/examples/Tests_union/ElasticSQ_test/ElasticSQ_test.instr new file mode 100644 index 0000000000..1cdf9fb6ed --- /dev/null +++ b/mcstas-comps/examples/Tests_union/ElasticSQ_test/ElasticSQ_test.instr @@ -0,0 +1,83 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: ElasticSQ_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 24.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of ElasticSQ_process: elastic scattering from S(|Q|) or single crystal S(Q) +* +* %Description +* A monochromatic beam hits a sample made from the Union process ElasticSQ_process +* (port of the MCViNE SQkernel and SvQkernel). +* crystal=0: liquid-like structure factor S(Q) with a first peak at Q0, given as an +* expression (1D table files "Q S" can be used via SQ_file). +* crystal=1: single crystal diffuse scattering, here rods of intensity along l through +* integer h and k of a simple cubic lattice (a=4 AA), rotated by A3. +* The banana monitor records the scattered intensity vs. scattering angle, the PSD +* monitor behind the sample shows the diffuse pattern for crystal=1. +* +* %Example: crystal=0 Detector: theta_I=1520 +* +* %Parameters +* Ei: [meV] Incident energy +* crystal: [1] 0: S(|Q|) liquid, 1: single crystal diffuse rods +* Q0: [AA^-1] Position of the first S(Q) peak (crystal=0) +* A3: [deg] Crystal rotation around the vertical axis (crystal=1) +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT ElasticSQ_test(Ei=20, int crystal=0, Q0=2, A3=10) + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.002, dist=1, focus_xw=0.01, focus_yh=0.02, E0=Ei, dE=0, flux=1e10) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +COMPONENT liquid = ElasticSQ_process(SQ_expr="1+1.5*exp(-(Q-p1)^2/0.08)-0.4*exp(-(Q-1.6*p1)^2/0.3)", crystal=0, + sigma=5, unit_cell_volume=30, p1=Q0) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT diffuse = ElasticSQ_process(SQ_expr="exp(-(sin(pi*h)^2+sin(pi*k)^2)/0.02)", crystal=1, + sigma=5, unit_cell_volume=64, ax=4, by=4, cz=4) +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, A3, 0) RELATIVE sample_pos + +COMPONENT mat_liquid = Union_make_material(my_absorption=0, process_string="liquid") +AT (0, 0, 0) ABSOLUTE + +COMPONENT mat_crystal = Union_make_material(my_absorption=0, process_string="diffuse") +AT (0, 0, 0) ABSOLUTE + +COMPONENT sample = Union_cylinder(radius=0.005, yheight=0.02, material_string=crystal ? "mat_crystal" : "mat_liquid", priority=1) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT theta = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, options="banana, theta limits=[5 175] bins=170", filename="theta.dat") +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT psd = PSD_monitor(xwidth=1, yheight=1, nx=200, ny=200, filename="psd.dat", restore_neutron=1) +AT (0, 0, 0.5) RELATIVE sample_pos + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/examples/Tests_union/IncoherentElastic_test/IncoherentElastic_test.instr b/mcstas-comps/examples/Tests_union/IncoherentElastic_test/IncoherentElastic_test.instr new file mode 100644 index 0000000000..048a2bd102 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/IncoherentElastic_test/IncoherentElastic_test.instr @@ -0,0 +1,71 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: IncoherentElastic_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 24.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of IncoherentElastic_process: incoherent elastic scattering with Debye-Waller factor +* +* %Description +* A monochromatic beam hits a vanadium cylinder made from the Union process +* IncoherentElastic_process (port of the MCViNE IncoherentElastic kernel). +* The banana monitor shows the angular distribution, which falls off as +* exp(-DW_core*Q^2) with Q = 2k sin(theta). With DW_core=0 the scattering is isotropic, +* as for Incoherent_process. The default DW_core is the value for vanadium at 300 K +* (from the V-dos.dat DOS used in IncoherentOnePhonon_test). +* +* %Example: DW_core=0.00669 Detector: theta_mon_I=1.29e+04 +* +* %Parameters +* Ei: [meV] Incident energy +* DW_core: [AA^2] Debye-Waller core, 2W = DW_core*Q^2 +* radius: [m] Radius of the vanadium cylinder +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT IncoherentElastic_test(Ei=100, DW_core=0.00669, radius=0.005) + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.005, dist=1, focus_xw=2*radius, focus_yh=0.03, E0=Ei, dE=0, flux=1e10) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +// V: sigma_inc = 5.08 barn, V = 13.827 AA^3 per atom +COMPONENT V_elastic = IncoherentElastic_process(sigma_inc=5.08, unit_cell_volume=13.827, DW_core=DW_core) +AT (0, 0, 0) RELATIVE sample_pos + +// V: sigma_abs = 5.08 barn +COMPONENT V = Union_make_material(my_absorption=100*5.08/13.827, process_string="V_elastic") +AT (0, 0, 0) ABSOLUTE + +COMPONENT sample = Union_cylinder(radius=radius, yheight=0.03, material_string="V", priority=1) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT theta_mon = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, + options="banana, theta limits=[5 175] bins=170", filename="theta.dat") +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/examples/Tests_union/IncoherentOnePhonon_test/IncoherentOnePhonon_test.instr b/mcstas-comps/examples/Tests_union/IncoherentOnePhonon_test/IncoherentOnePhonon_test.instr new file mode 100644 index 0000000000..ef1019e745 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/IncoherentOnePhonon_test/IncoherentOnePhonon_test.instr @@ -0,0 +1,100 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: IncoherentOnePhonon_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 23.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of IncoherentOnePhonon_process: incoherent one-phonon scattering from a DOS +* +* %Description +* A monochromatic beam hits a vanadium cylinder made from the Union process +* IncoherentOnePhonon_process (port of the MCViNE IncoherentInelastic kernel), +* using the vanadium phonon DOS from MCViNE (V-dos.dat, frequencies in THz). +* The Debye-Waller factor is calculated from the DOS. +* +* A banana monitor records scattering angle vs. final energy, and an energy +* monitor at scattering angle tth records the energy spectrum. With focus=1 the +* scattering is focused onto the energy monitor. With dEf_focus > 0 only final +* energies in [Ef_focus-dEf_focus/2, Ef_focus+dEf_focus/2] are generated. +* +* %Example: focus=1 Detector: Edet_I=9.91 +* +* %Parameters +* Ei: [meV] Incident energy +* T: [K] Sample temperature +* tth: [deg] Scattering angle of the energy monitor +* focus: [1] 0: scatter into 4pi, 1: focus scattering onto the energy monitor +* Ef_focus: [meV] Centre of the final energy window (energy focusing) +* dEf_focus: [meV] Width of the final energy window, 0: no energy focusing +* radius: [m] Radius of the vanadium cylinder +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT IncoherentOnePhonon_test(Ei=50, T=300, tth=60, int focus=0, Ef_focus=0, dEf_focus=0, radius=0.005) + +DECLARE +%{ + int target; + double fw; +%} + +INITIALIZE +%{ + target = focus ? 3 : 0; + fw = focus ? 0.05 : 0; +%} + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.005, dist=1, focus_xw=0.012, focus_yh=0.03, E0=Ei, dE=0, flux=1e10) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +// V: sigma_inc = 5.08 barn, V = 13.827 AA^3 per atom +COMPONENT V_phonons = IncoherentOnePhonon_process(DOS_file="V-dos.dat", sigma_inc=5.08, unit_cell_volume=13.827, mass=50.94, T=T, + Ef_focus=Ef_focus, dEf_focus=dEf_focus) +AT (0, 0, 0) RELATIVE sample_pos + +// V: sigma_abs = 5.08 barn +COMPONENT V = Union_make_material(my_absorption=100*5.08/13.827, process_string="V_phonons") +AT (0, 0, 0) ABSOLUTE + +COMPONENT det_arm = Arm() +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, tth, 0) RELATIVE sample_pos + +COMPONENT sample = Union_cylinder(radius=radius, yheight=0.03, material_string="V", priority=1, + target_index=target, focus_xw=fw, focus_xh=fw) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Etheta = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, + options="banana, theta limits=[5 170] bins=165, energy limits=[0 100] bins=200", filename="Etheta.dat") +WHEN (!focus) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Edet = E_monitor(xwidth=0.05, yheight=0.05, nE=200, Emin=0, Emax=100, filename="Edet.dat", restore_neutron=1) +AT (0, 0, 1) RELATIVE det_arm + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/examples/Tests_union/IncoherentOnePhonon_test/V-dos.dat b/mcstas-comps/examples/Tests_union/IncoherentOnePhonon_test/V-dos.dat new file mode 100644 index 0000000000..0eccf7d71f --- /dev/null +++ b/mcstas-comps/examples/Tests_union/IncoherentOnePhonon_test/V-dos.dat @@ -0,0 +1,384 @@ +# Vanadium phonon density of states. +# x axis: frequency in terahz (note: not angular frequency) +# Frequency(TeraHz) DOS +0.0 0.0 +0.02 0.0 +0.04 0.0 +0.06 0.0 +0.08 0.0 +0.1 0.0 +0.12 0.0 +0.14 0.0 +0.16 0.0 +0.18 0.0 +0.2 0.0 +0.22 0.000520833333333 +0.24 0.000520833333333 +0.26 0.00182291666667 +0.28 0.00104166666667 +0.3 0.000260416666667 +0.32 0.0 +0.34 0.00104166666667 +0.36 0.00078125 +0.38 0.00182291666667 +0.4 0.000260416666667 +0.42 0.000260416666667 +0.44 0.0015625 +0.46 0.00208333333333 +0.48 0.00104166666667 +0.5 0.00208333333333 +0.52 0.00104166666667 +0.54 0.00208333333333 +0.56 0.00260416666667 +0.58 0.00182291666667 +0.6 0.00390625 +0.62 0.00182291666667 +0.64 0.00286458333333 +0.66 0.00338541666667 +0.68 0.00286458333333 +0.7 0.00442708333333 +0.72 0.00520833333333 +0.74 0.00364583333333 +0.76 0.00442708333333 +0.78 0.00338541666667 +0.8 0.0046875 +0.82 0.00364583333333 +0.84 0.00442708333333 +0.86 0.00520833333333 +0.88 0.00390625 +0.9 0.00546875 +0.92 0.00416666666667 +0.94 0.00364583333333 +0.96 0.00442708333333 +0.98 0.00677083333333 +1.0 0.003125 +1.02 0.00572916666667 +1.04 0.00442708333333 +1.06 0.0078125 +1.08 0.00442708333333 +1.1 0.00598958333333 +1.12 0.00494791666667 +1.14 0.0104166666667 +1.16 0.00963541666667 +1.18 0.01171875 +1.2 0.00651041666667 +1.22 0.00833333333333 +1.24 0.0119791666667 +1.26 0.0106770833333 +1.28 0.009375 +1.3 0.01015625 +1.32 0.009375 +1.34 0.009375 +1.36 0.00807291666667 +1.38 0.00963541666667 +1.4 0.0130208333333 +1.42 0.0130208333333 +1.44 0.0122395833333 +1.46 0.0122395833333 +1.48 0.0135416666667 +1.5 0.015625 +1.52 0.0153645833333 +1.54 0.0153645833333 +1.56 0.015625 +1.58 0.01640625 +1.6 0.0143229166667 +1.62 0.0184895833333 +1.64 0.0203125 +1.66 0.0153645833333 +1.68 0.0143229166667 +1.7 0.015625 +1.72 0.0174479166667 +1.74 0.01796875 +1.76 0.0145833333333 +1.78 0.0200520833333 +1.8 0.0213541666667 +1.82 0.0200520833333 +1.84 0.0135416666667 +1.86 0.0174479166667 +1.88 0.0200520833333 +1.9 0.0208333333333 +1.92 0.0197916666667 +1.94 0.0184895833333 +1.96 0.02109375 +1.98 0.0234375 +2.0 0.0192708333333 +2.02 0.0205729166667 +2.04 0.0213541666667 +2.06 0.0239583333333 +2.08 0.0265625 +2.1 0.02734375 +2.12 0.0236979166667 +2.14 0.02578125 +2.16 0.0268229166667 +2.18 0.0286458333333 +2.2 0.0296875 +2.22 0.0265625 +2.24 0.0252604166667 +2.26 0.03203125 +2.28 0.03046875 +2.3 0.0286458333333 +2.32 0.028125 +2.34 0.0341145833333 +2.36 0.0333333333333 +2.38 0.0341145833333 +2.4 0.0330729166667 +2.42 0.0346354166667 +2.44 0.0309895833333 +2.46 0.0375 +2.48 0.0354166666667 +2.5 0.0346354166667 +2.52 0.0411458333333 +2.54 0.03984375 +2.56 0.0401041666667 +2.58 0.0356770833333 +2.6 0.0424479166667 +2.62 0.0354166666667 +2.64 0.0375 +2.66 0.03984375 +2.68 0.0427083333333 +2.7 0.0424479166667 +2.72 0.0432291666667 +2.74 0.0401041666667 +2.76 0.0447916666667 +2.78 0.0395833333333 +2.8 0.0484375 +2.82 0.0419270833333 +2.84 0.0455729166667 +2.86 0.0466145833333 +2.88 0.0486979166667 +2.9 0.0473958333333 +2.92 0.0486979166667 +2.94 0.0510416666667 +2.96 0.0434895833333 +2.98 0.0518229166667 +3.0 0.0559895833333 +3.02 0.0591145833333 +3.04 0.0526041666667 +3.06 0.04765625 +3.08 0.0513020833333 +3.1 0.0557291666667 +3.12 0.0614583333333 +3.14 0.0572916666667 +3.16 0.0598958333333 +3.18 0.0609375 +3.2 0.0565104166667 +3.22 0.0528645833333 +3.24 0.0572916666667 +3.26 0.05703125 +3.28 0.0552083333333 +3.3 0.06328125 +3.32 0.0578125 +3.34 0.0606770833333 +3.36 0.0708333333333 +3.38 0.0638020833333 +3.4 0.0666666666667 +3.42 0.06875 +3.44 0.071875 +3.46 0.0721354166667 +3.48 0.0645833333333 +3.5 0.0690104166667 +3.52 0.0713541666667 +3.54 0.0734375 +3.56 0.0765625 +3.58 0.0760416666667 +3.6 0.0684895833333 +3.62 0.0713541666667 +3.64 0.0783854166667 +3.66 0.0739583333333 +3.68 0.0872395833333 +3.7 0.0877604166667 +3.72 0.0817708333333 +3.74 0.0815104166667 +3.76 0.0880208333333 +3.78 0.0861979166667 +3.8 0.08828125 +3.82 0.090625 +3.84 0.096875 +3.86 0.1 +3.88 0.11015625 +3.9 0.105989583333 +3.92 0.0942708333333 +3.94 0.105989583333 +3.96 0.0971354166667 +3.98 0.10078125 +4.0 0.106770833333 +4.02 0.10625 +4.04 0.11875 +4.06 0.10234375 +4.08 0.124479166667 +4.1 0.130208333333 +4.12 0.16171875 +4.14 0.161458333333 +4.16 0.158333333333 +4.18 0.158854166667 +4.2 0.182552083333 +4.22 0.1890625 +4.24 0.196354166667 +4.26 0.203125 +4.28 0.232291666667 +4.3 0.24296875 +4.32 0.2609375 +4.34 0.276822916667 +4.36 0.270833333333 +4.38 0.2734375 +4.4 0.28671875 +4.42 0.286458333333 +4.44 0.27421875 +4.46 0.2828125 +4.48 0.2828125 +4.5 0.276302083333 +4.52 0.290104166667 +4.54 0.293489583333 +4.56 0.30546875 +4.58 0.306510416667 +4.6 0.280729166667 +4.62 0.299739583333 +4.64 0.311979166667 +4.66 0.3203125 +4.68 0.321354166667 +4.7 0.32890625 +4.72 0.322135416667 +4.74 0.28828125 +4.76 0.292447916667 +4.78 0.275 +4.8 0.265625 +4.82 0.271614583333 +4.84 0.263802083333 +4.86 0.251822916667 +4.88 0.251822916667 +4.9 0.245833333333 +4.92 0.253645833333 +4.94 0.253645833333 +4.96 0.23984375 +4.98 0.236979166667 +5.0 0.226041666667 +5.02 0.22265625 +5.04 0.238802083333 +5.06 0.24375 +5.08 0.22734375 +5.1 0.2265625 +5.12 0.236197916667 +5.14 0.2265625 +5.16 0.20703125 +5.18 0.21640625 +5.2 0.221875 +5.22 0.19453125 +5.24 0.204427083333 +5.26 0.19375 +5.28 0.173697916667 +5.3 0.182291666667 +5.32 0.179166666667 +5.34 0.182552083333 +5.36 0.17734375 +5.38 0.19296875 +5.4 0.18671875 +5.42 0.18828125 +5.44 0.184895833333 +5.46 0.192447916667 +5.48 0.19609375 +5.5 0.182552083333 +5.52 0.1953125 +5.54 0.195052083333 +5.56 0.196875 +5.58 0.20234375 +5.6 0.200520833333 +5.62 0.210677083333 +5.64 0.21015625 +5.66 0.215625 +5.68 0.22734375 +5.7 0.220833333333 +5.72 0.249479166667 +5.74 0.238802083333 +5.76 0.25078125 +5.78 0.236979166667 +5.8 0.251822916667 +5.82 0.238020833333 +5.84 0.256770833333 +5.86 0.27421875 +5.88 0.308333333333 +5.9 0.290885416667 +5.92 0.282291666667 +5.94 0.277604166667 +5.96 0.277083333333 +5.98 0.270833333333 +6.0 0.26640625 +6.02 0.262239583333 +6.04 0.265104166667 +6.06 0.261979166667 +6.08 0.267708333333 +6.1 0.269010416667 +6.12 0.28125 +6.14 0.331770833333 +6.16 0.358854166667 +6.18 0.363020833333 +6.2 0.345052083333 +6.22 0.315885416667 +6.24 0.29765625 +6.26 0.286458333333 +6.28 0.2984375 +6.3 0.280989583333 +6.32 0.28359375 +6.34 0.264583333333 +6.36 0.29140625 +6.38 0.305729166667 +6.4 0.27890625 +6.42 0.310416666667 +6.44 0.294791666667 +6.46 0.329166666667 +6.48 0.332291666667 +6.5 0.37421875 +6.52 0.390104166667 +6.54 0.406770833333 +6.56 0.390364583333 +6.58 0.373697916667 +6.6 0.3515625 +6.62 0.357291666667 +6.64 0.3609375 +6.66 0.359375 +6.68 0.366927083333 +6.7 0.363020833333 +6.72 0.348177083333 +6.74 0.349479166667 +6.76 0.385677083333 +6.78 0.370052083333 +6.8 0.35390625 +6.82 0.35546875 +6.84 0.352083333333 +6.86 0.372395833333 +6.88 0.35546875 +6.9 0.355729166667 +6.92 0.368489583333 +6.94 0.371354166667 +6.96 0.37109375 +6.98 0.371614583333 +7.0 0.361458333333 +7.02 0.336458333333 +7.04 0.276041666667 +7.06 0.247395833333 +7.08 0.236979166667 +7.1 0.211197916667 +7.12 0.184895833333 +7.14 0.16953125 +7.16 0.175520833333 +7.18 0.17421875 +7.2 0.173958333333 +7.22 0.160416666667 +7.24 0.168229166667 +7.26 0.148958333333 +7.28 0.149739583333 +7.3 0.144791666667 +7.32 0.144791666667 +7.34 0.129947916667 +7.36 0.112760416667 +7.38 0.133333333333 +7.4 0.0786458333333 +7.42 0.0328125 +7.44 0.00104166666667 +7.46 0.0 +7.48 0.0 +7.5 0.0 +7.52 0.0 +7.54 0.0 +7.56 0.0 +7.58 0.0 +7.6 0.0 \ No newline at end of file diff --git a/mcstas-comps/examples/Tests_union/IsotropicSqw_test/IsotropicSqw_test.instr b/mcstas-comps/examples/Tests_union/IsotropicSqw_test/IsotropicSqw_test.instr new file mode 100644 index 0000000000..97e1d2a8a7 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/IsotropicSqw_test/IsotropicSqw_test.instr @@ -0,0 +1,99 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: IsotropicSqw_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 24.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of IsotropicSqw_process: isotropic S(q,w) from a McStas .sqw file (liquid He4) +* +* %Description +* A monochromatic beam hits a cylinder of liquid He4 made from the Union process +* IsotropicSqw_process (port of the MCViNE SQEkernel / GridSQE), reading +* data/He4_liq_coh.sqw. Cross section, density and temperature are taken from the +* file header, and S is normalised with the sum rule as in Isotropic_Sqw (norm=-1). +* With comp_select=1 the standalone Isotropic_Sqw component is used +* instead for comparison. method=0 samples (Q,E) as MCViNE, method=1 uses Union +* focusing (focus=1: onto the energy monitor at angle tth). +* +* %Example: method=1 focus=1 Detector: Edet_I=2.81 +* +* %Parameters +* Ei: [meV] Incident energy +* tth: [deg] Scattering angle of the energy monitor +* method: [1] 0: MCViNE sampling (4pi), 1: focusing +* focus: [1] With method=1: 0 4pi, 1 focus onto the energy monitor +* comp_select: [1] 0: IsotropicSqw_process (Union), 1: Isotropic_Sqw +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT IsotropicSqw_test(Ei=4, tth=60, int method=1, int focus=0, int comp_select=0) + +DECLARE +%{ + int target; + double fw; +%} + +INITIALIZE +%{ + target = (focus && method == 1) ? 4 : 0; + fw = (focus && method == 1) ? 0.05 : 0; +%} + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.002, dist=1, focus_xw=0.01, focus_yh=0.02, E0=Ei, dE=0, flux=1e10) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +COMPONENT He4 = IsotropicSqw_process(Sqw_file="He4_liq_coh.sqw", method=method, norm=-1) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT He4_mat = Union_make_material(my_absorption=100*0.00747*0.072, process_string="He4") +AT (0, 0, 0) ABSOLUTE + +COMPONENT det_arm = Arm() +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, tth, 0) RELATIVE sample_pos + +COMPONENT sample = Union_cylinder(radius=0.005, yheight=0.02, material_string="He4_mat", priority=1, + target_index=target, focus_xw=fw, focus_xh=fw) +WHEN (comp_select == 0) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +WHEN (comp_select == 0) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT standalone = Isotropic_Sqw(Sqw_coh="He4_liq_coh.sqw", radius=0.005, yheight=0.02) +WHEN (comp_select == 1) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Etheta = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, + options="banana, theta limits=[5 135] bins=130, energy limits=[2 6] bins=100", filename="Etheta.dat") +WHEN (!focus) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Edet = E_monitor(xwidth=0.05, yheight=0.05, nE=200, Emin=2, Emax=6, filename="Edet.dat", restore_neutron=1) +AT (0, 0, 1) RELATIVE det_arm + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/examples/Tests_union/Resolution_test/Resolution_test.instr b/mcstas-comps/examples/Tests_union/Resolution_test/Resolution_test.instr new file mode 100644 index 0000000000..78186d59cc --- /dev/null +++ b/mcstas-comps/examples/Tests_union/Resolution_test/Resolution_test.instr @@ -0,0 +1,106 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: Resolution_test +* +* %Identification +* Written by: Peter Willendrup +* Date: 24.09.2026 +* Origin: DTU Physics +* %INSTRUMENT_SITE: Tests_union +* +* Test of Resolution_process: resolution / test scattering kernels in a Union sample +* +* %Description +* A monochromatic beam hits a sample made from the Union process Resolution_process +* (port of the MCViNE ConstantEnergyTransfer, ConstantQE, ConstantvQE and DGSSXRes +* kernels), inside an aluminium container to show a resolution calculation including +* the sample environment. +* mode=0: constant energy transfer E, isotropic (or focused) +* mode=1: constant |Q| and E, i.e. a cone of scattering (powder line) +* mode=2: constant Q vector (along x, Q=Qx) and Gaussian distribution of E +* mode=3: flat time of flight window [tof-dtof/2, tof+dtof/2] at the focusing target +* (the energy monitor) +* +* %Example: mode=0 Detector: Etheta_I=712 +* +* %Parameters +* Ei: [meV] Incident energy +* mode: [1] Resolution_process mode 0-3 +* E: [meV] Energy transfer (modes 0-2) +* Q: [AA^-1] Momentum transfer (mode 1: |Q|, mode 2: Qx) +* tth: [deg] Scattering angle of the energy monitor (focusing target for mode 3) +* +* %Link +* MCViNE: https://github.com/mcvine/mcvine +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT Resolution_test(Ei=30, int mode=0, E=10, Q=3, tth=60) + +DECLARE +%{ + int target; + double fw; + double tof_target; +%} + +INITIALIZE +%{ + target = (mode == 3) ? 4 : 0; + fw = (mode == 3) ? 0.05 : 0; + // time of flight for elastic scattering to the energy monitor (2 m flight path) + tof_target = 2.0 / (SE2V * sqrt (Ei)); +%} + +TRACE + +COMPONENT init = Union_init() +AT (0, 0, 0) ABSOLUTE + +COMPONENT origin = Progress_bar() +AT (0, 0, 0) ABSOLUTE + +COMPONENT source = Source_simple(radius=0.002, dist=1, focus_xw=0.01, focus_yh=0.02, E0=Ei, dE=0, flux=1e10) +AT (0, 0, 0) RELATIVE origin + +COMPONENT sample_pos = Arm() +AT (0, 0, 1) RELATIVE origin + +COMPONENT res = Resolution_process(mode=mode, E=E, Q=Q, Qx=Q, Qy=0, Qz=0, dE=0.5, tof=tof_target, dtof=0.1*tof_target, sigma=5, unit_cell_volume=40) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Al_inc = Incoherent_process(sigma=0.0082, unit_cell_volume=66.4) +AT (0, 0, 0) ABSOLUTE + +COMPONENT sample_mat = Union_make_material(my_absorption=0, process_string="res") +AT (0, 0, 0) ABSOLUTE + +COMPONENT Al = Union_make_material(my_absorption=100*4*0.231/66.4, process_string="Al_inc") +AT (0, 0, 0) ABSOLUTE + +COMPONENT det_arm = Arm() +AT (0, 0, 0) RELATIVE sample_pos +ROTATED (0, tth, 0) RELATIVE sample_pos + +COMPONENT sample = Union_cylinder(radius=0.004, yheight=0.02, material_string="sample_mat", priority=2, + target_index=target, focus_xw=fw, focus_xh=fw) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT can = Union_cylinder(radius=0.005, yheight=0.03, material_string="Al", priority=1) +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT master = Union_master() +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Etheta = Monitor_nD(radius=1, yheight=0.2, restore_neutron=1, + options="banana, theta limits=[5 175] bins=170, energy limits=[0 60] bins=120", filename="Etheta.dat") +AT (0, 0, 0) RELATIVE sample_pos + +COMPONENT Edet = Monitor_nD(xwidth=0.05, yheight=0.05, restore_neutron=1, options="energy limits=[0 60] bins=240", filename="Edet.dat") +AT (0, 0, 1) RELATIVE det_arm + +COMPONENT stop = Union_stop() +AT (0, 0, 0) ABSOLUTE + +END diff --git a/mcstas-comps/share/phonon-dispersion-lib.c b/mcstas-comps/share/phonon-dispersion-lib.c new file mode 100644 index 0000000000..a41a095e9b --- /dev/null +++ b/mcstas-comps/share/phonon-dispersion-lib.c @@ -0,0 +1,775 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* Library: share/phonon-dispersion-lib.c +* +* Written by: Peter Willendrup, port of MCViNE (J.Y.Y. Lin et al.) phonon kernel code +* Date: 23.09.2026 +* Origin: DTU Physics +* +* Shared code for the Union processes CoherentPhononSingleXtal_process and +* CoherentPhononPowder_process: +* - readers for phonon dispersions in the MCViNE IDF format (Qgridinfo, Omega2, +* Polarizations, DOS) and MCViNE/diffpy xyz crystal structure files +* - periodic, trilinearly interpolated phonon energies and polarizations +* (PeriodicDispersion_3D, ChangeCoordinateSystem_forDispersion_3D and +* LinearlyInterpolatedDispersionOnGrid_3D in MCViNE) +* - one-phonon coherent structure factor (scattering_length.icc) +* - Debye-Waller factor from the DOS (DWFromDOS.icc) +* - root finding for omega(Q) = |Ei-Ef| (Omega_minus_deltaE, FindRootsEvenly, ZRidd) +* +* Usage: %include "read_table-lib" and then %include "phonon-dispersion-lib" in SHARE. Error messages are prefixed +* with the string "comp" passed by the caller. +* +*******************************************************************************/ + +#ifndef PHONON_DISPERSION_LIB_H +#include "phonon-dispersion-lib.h" +#endif + +#ifndef PHONON_DISPERSION_LIB_C +#define PHONON_DISPERSION_LIB_C + +// --------------------------------------------------------------------------- +// Grid lookup with periodic folding (PeriodicDispersion_3D + +// ChangeCoordinateSystem + LinearlyInterpolatedDispersionOnGrid_3D in MCViNE) +// --------------------------------------------------------------------------- +void +phdisp_grid_coords (struct phdisp_struct* s, double* Q, int* i0, double* r) { + int d; + for (d = 0; d < 3; d++) { + double f = s->binv[d][0] * Q[0] + s->binv[d][1] * Q[1] + s->binv[d][2] * Q[2]; + f -= floor (f); + if (f >= 1.0 || f < 0.0) + f = 0.0; + double x = f * (s->n[d] - 1); + int i = (int)floor (x); + double rem = x - i; + if (i >= s->n[d] - 1) { + i = s->n[d] - 2; + rem = 1.0; + } + i0[d] = i; + r[d] = rem; + } +} + +// trilinear interpolation of a quantity stored at data[q*stride + offset] +double +phdisp_interp (struct phdisp_struct* s, double* data, long stride, long offset, int* i0, double* r) { + double res = 0; + int dx, dy, dz; + for (dx = 0; dx < 2; dx++) + for (dy = 0; dy < 2; dy++) + for (dz = 0; dz < 2; dz++) { + double w = (dx ? r[0] : 1 - r[0]) * (dy ? r[1] : 1 - r[1]) * (dz ? r[2] : 1 - r[2]); + if (w == 0) + continue; + long q = ((long)(i0[0] + dx) * s->n[1] + (i0[1] + dy)) * s->n[2] + (i0[2] + dz); + res += w * data[q * stride + offset]; + } + return res; +} + +double +phdisp_energy (struct phdisp_struct* s, int branch, double* Q) { + int i0[3]; + double r[3]; + phdisp_grid_coords (s, Q, i0, r); + return phdisp_interp (s, s->E, s->n_branches, branch, i0, r); +} + +void +phdisp_polarization (struct phdisp_struct* s, int branch, int atom, double* Q, double* re, double* im) { + int i0[3], d; + double r[3]; + long stride = (long)s->n_branches * s->n_atoms * 6; + long base = ((long)branch * s->n_atoms + atom) * 6; + phdisp_grid_coords (s, Q, i0, r); + for (d = 0; d < 3; d++) { + re[d] = phdisp_interp (s, s->pol, stride, base + 2 * d, i0, r); + im[d] = phdisp_interp (s, s->pol, stride, base + 2 * d + 1, i0, r); + } +} + +// omega(Q) - |Ei-Ef| as a function of the final speed (Omega_minus_deltaE.cc) +double +phdisp_omega_minus_dE (double vf, struct phdisp_omega_ctx* c) { + double Q[3]; + int d; + for (d = 0; d < 3; d++) + Q[d] = V2K * (c->vi[d] - vf * c->vf_dir[d]); + if (c->powder) { + double ql = sqrt (Q[0] * Q[0] + Q[1] * Q[1] + Q[2] * Q[2]); + for (d = 0; d < 3; d++) + Q[d] = ql * c->qhat[d]; + } + return phdisp_energy (c->s, c->branch, Q) - VS2E * fabs (c->vi_l * c->vi_l - vf * vf); +} + +// Ridder's method (mccomponents/math/rootfinding.cc). Returns 1 on success. +// fl, fh: function values at x1, x2 (already evaluated by the caller) +int +phdisp_zridd (struct phdisp_omega_ctx* c, double x1, double x2, double fl, double fh, double xacc, double* root) { + const double UNUSED = -1.11e30; + double ans, fm, fnew, s, xh, xl, xm, xnew, tmp; + int j; + if (fl * fh >= 0) { + if (fl == 0) { + *root = x1; + return 1; + } + if (fh == 0) { + *root = x2; + return 1; + } + return 0; + } + int converged = 0; + xl = x1; + xh = x2; + ans = UNUSED; + for (j = 1; j < 60 && !converged; j++) { + xm = 0.5 * (xl + xh); + fm = phdisp_omega_minus_dE (xm, c); + tmp = fm * fm - fh * fl; + if (tmp >= 0) { + s = sqrt (tmp); + if (s == 0.0) { + converged = 1; + break; + } + xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); + } else { + xnew = xm; + } + if (fabs (xnew - ans) <= xacc) { + converged = 1; + break; + } + ans = xnew; + fnew = phdisp_omega_minus_dE (ans, c); + if (fnew == 0.0) { + converged = 1; + break; + } + if ((fnew >= 0 ? fabs (fm) : -fabs (fm)) != fm) { + xl = xm; + fl = fm; + xh = ans; + fh = fnew; + } else if ((fnew >= 0 ? fabs (fl) : -fabs (fl)) != fl) { + xh = ans; + fh = fnew; + } else if ((fnew >= 0 ? fabs (fh) : -fabs (fh)) != fh) { + xl = ans; + fl = fnew; + } else { + return 0; + } + if (fabs (xh - xl) <= xacc) + converged = 1; + } + if (!converged || ans == UNUSED) + return 0; + *root = ans; + return 1; +} + +// FindRootsEvenly: search roots in nsteps equal intervals of [x1,x2]. +// The function is evaluated once per interval boundary and shared by neighbouring intervals. +int +phdisp_find_roots (struct phdisp_omega_ctx* c, double x1, double x2, int nsteps, double xacc, double* roots) { + int i, nf = 0; + double step = (x2 - x1) / nsteps, root; + double fl = phdisp_omega_minus_dE (x1, c), fh; + for (i = 0; i < nsteps; i++) { + fh = phdisp_omega_minus_dE (x1 + step * (i + 1), c); + if (phdisp_zridd (c, x1 + step * i, x1 + step * (i + 1), fl, fh, xacc, &root)) + roots[nf++] = root; + fl = fh; + } + return nf; +} + +// Bose factor: n+1 for phonon creation (omega>0), n for annihilation +double +phdisp_bose (double omega, double T) { + if (omega == 0.0) + return 1.0; + double n = 1.0 / (exp (fabs (omega) / (T * PHDISP_T2E)) - 1.0); + return omega > 0 ? 1.0 + n : n; +} + +// --------------------------------------------------------------------------- +// File readers +// --------------------------------------------------------------------------- +FILE* +phdisp_open (const char* dir, const char* name, const char* mode, const char* comp) { + char path[1024]; + FILE* fp; + if (dir && strlen (dir)) + snprintf (path, 1024, "%s%c%s", dir, MC_PATHSEP_C, name); + else + snprintf (path, 1024, "%s", name); + fp = Open_File (path, mode, NULL); + if (!fp) { + fprintf (stderr, "%s: ERROR: cannot open file %s\n", comp, path); + exit (-1); + } + return fp; +} + +// IDF header: 64 char filetype, int version, 1024 char comment +void +phdisp_idf_header (FILE* fp, const char* expected, const char* comp) { + char filetype[65]; + char comment[1024]; + int version; + if (fread (filetype, 1, 64, fp) != 64 || fread (&version, sizeof (int), 1, fp) != 1 || fread (comment, 1, 1024, fp) != 1024) { + fprintf (stderr, "%s: ERROR: could not read IDF header of %s file\n", comp, expected); + exit (-1); + } + filetype[64] = '\0'; + if (strncmp (filetype, expected, strlen (expected))) { + fprintf (stderr, "%s: ERROR: file type is '%s', expected '%s'\n", comp, filetype, expected); + exit (-1); + } +} + +// --- tiny evaluator for the python-like Qgridinfo file --- + +void +phdisp_skipws (struct phdisp_parser* P) { + while (*P->p == ' ' || *P->p == '\t' || *P->p == '\r') + P->p++; +} + +struct phdisp_val phdisp_parse_list (struct phdisp_parser* P); +struct phdisp_val phdisp_parse_expr (struct phdisp_parser* P); + +struct phdisp_val +phdisp_scalar (double x) { + struct phdisp_val r; + r.n = 1; + r.v[0] = x; + r.v[1] = r.v[2] = 0; + return r; +} + +struct phdisp_val +phdisp_binop (struct phdisp_parser* P, struct phdisp_val a, struct phdisp_val b, char op) { + struct phdisp_val r; + int i, n = a.n > b.n ? a.n : b.n; + if (a.n != b.n && a.n != 1 && b.n != 1) { + P->error = 1; + return a; + } + r.n = n; + for (i = 0; i < 3; i++) { + double x = a.v[a.n == 1 ? 0 : i], y = b.v[b.n == 1 ? 0 : i]; + switch (op) { + case '+': + r.v[i] = x + y; + break; + case '-': + r.v[i] = x - y; + break; + case '*': + r.v[i] = x * y; + break; + case '/': + r.v[i] = x / y; + break; + case '^': + r.v[i] = pow (x, y); + break; + } + } + return r; +} + +struct phdisp_val +phdisp_parse_primary (struct phdisp_parser* P) { + struct phdisp_val r = phdisp_scalar (0); + phdisp_skipws (P); + if (*P->p == '(' || *P->p == '[') { + char close = (*P->p == '(') ? ')' : ']'; + P->p++; + r = phdisp_parse_list (P); + phdisp_skipws (P); + if (*P->p != close) + P->error = 1; + else + P->p++; + return r; + } + if (isdigit (*P->p) || *P->p == '.') { + char* end; + r = phdisp_scalar (strtod (P->p, &end)); + P->p = end; + return r; + } + if (isalpha (*P->p) || *P->p == '_') { + char name[64]; + int len = 0, i; + while ((isalnum (*P->p) || *P->p == '_' || *P->p == '.') && len < 63) + name[len++] = *P->p++; + name[len] = '\0'; + const char* fname = strncmp (name, "math.", 5) == 0 ? name + 5 : name; + phdisp_skipws (P); + if (*P->p == '(') { // function call + P->p++; + struct phdisp_val a = phdisp_parse_expr (P); + phdisp_skipws (P); + if (*P->p == ')') + P->p++; + else + P->error = 1; + double x = a.v[0]; + if (!strcmp (fname, "sqrt")) + return phdisp_scalar (sqrt (x)); + if (!strcmp (fname, "sin")) + return phdisp_scalar (sin (x)); + if (!strcmp (fname, "cos")) + return phdisp_scalar (cos (x)); + if (!strcmp (fname, "tan")) + return phdisp_scalar (tan (x)); + if (!strcmp (fname, "abs") || !strcmp (fname, "fabs")) + return phdisp_scalar (fabs (x)); + if (!strcmp (fname, "float")) + return phdisp_scalar (x); + P->error = 1; + return a; + } + if (!strcmp (fname, "pi")) + return phdisp_scalar (PI); + for (i = P->nvars - 1; i >= 0; i--) + if (!strcmp (P->names[i], name)) + return P->vals[i]; + fprintf (stderr, "phonon-dispersion-lib: Qgridinfo: unknown name '%s'\n", name); + P->error = 1; + return r; + } + P->error = 1; + return r; +} + +struct phdisp_val phdisp_parse_unary (struct phdisp_parser* P); + +struct phdisp_val +phdisp_parse_power (struct phdisp_parser* P) { + struct phdisp_val a = phdisp_parse_primary (P); + phdisp_skipws (P); + if (P->p[0] == '*' && P->p[1] == '*') { + P->p += 2; + a = phdisp_binop (P, a, phdisp_parse_unary (P), '^'); + } + return a; +} + +struct phdisp_val +phdisp_parse_unary (struct phdisp_parser* P) { + phdisp_skipws (P); + if (*P->p == '-') { + P->p++; + return phdisp_binop (P, phdisp_scalar (-1), phdisp_parse_unary (P), '*'); + } + if (*P->p == '+') { + P->p++; + return phdisp_parse_unary (P); + } + return phdisp_parse_power (P); +} + +struct phdisp_val +phdisp_parse_term (struct phdisp_parser* P) { + struct phdisp_val a = phdisp_parse_unary (P); + for (;;) { + phdisp_skipws (P); + if ((*P->p == '*' && P->p[1] != '*') || *P->p == '/') { + char op = *P->p++; + a = phdisp_binop (P, a, phdisp_parse_unary (P), op); + } else + return a; + } +} + +struct phdisp_val +phdisp_parse_expr (struct phdisp_parser* P) { + struct phdisp_val a = phdisp_parse_term (P); + for (;;) { + phdisp_skipws (P); + if (*P->p == '+' || *P->p == '-') { + char op = *P->p++; + a = phdisp_binop (P, a, phdisp_parse_term (P), op); + } else + return a; + } +} + +// comma separated list of scalars -> vector (up to 3 components) +struct phdisp_val +phdisp_parse_list (struct phdisp_parser* P) { + struct phdisp_val r, a; + r.n = 0; + r.v[0] = r.v[1] = r.v[2] = 0; + for (;;) { + phdisp_skipws (P); + if (*P->p == '\0' || *P->p == '\n' || *P->p == ')' || *P->p == ']') + break; + a = phdisp_parse_expr (P); + if (a.n == 1 && r.n < 3) + r.v[r.n++] = a.v[0]; + else if (r.n == 0) + r = a; + else + P->error = 1; + phdisp_skipws (P); + if (*P->p == ',') + P->p++; + else + break; + } + return r; +} + +void +phdisp_read_qgridinfo (const char* dir, struct phdisp_struct* s, const char* comp) { + FILE* fp = phdisp_open (dir, "Qgridinfo", "r", comp); + struct phdisp_parser* P = calloc (1, sizeof (struct phdisp_parser)); + char line[4096]; + int lineno = 0, i; + while (fgets (line, sizeof (line), fp)) { + char *c, *targets[16]; + int ntargets = 0; + lineno++; + if ((c = strchr (line, '#'))) + *c = '\0'; + c = line; + while (*c == ' ' || *c == '\t') + c++; + if (*c == '\0' || *c == '\n' || *c == '\r' || !strncmp (c, "import ", 7) || !strncmp (c, "from ", 5)) + continue; + // split "a = b = expr" on '=' + char* start = c; + while ((c = strchr (start, '=')) && ntargets < 16) { + *c = '\0'; + targets[ntargets++] = start; + start = c + 1; + } + if (ntargets == 0) { + fprintf (stderr, "%s: Qgridinfo line %d ignored: %s", comp, lineno, line); + continue; + } + P->p = start; + P->error = 0; + struct phdisp_val v = phdisp_parse_list (P); + if (P->error) { + fprintf (stderr, "%s: ERROR: could not parse Qgridinfo line %d\n", comp, lineno); + exit (-1); + } + for (i = 0; i < ntargets; i++) { + char name[64]; + if (sscanf (targets[i], " %63[A-Za-z0-9_]", name) != 1) + continue; + if (P->nvars < PHDISP_MAXVARS) { + strcpy (P->names[P->nvars], name); + P->vals[P->nvars++] = v; + } + } + } + fclose (fp); + const char* bn[3] = { "b1", "b2", "b3" }; + const char* nn[3] = { "n1", "n2", "n3" }; + int d, found; + for (d = 0; d < 3; d++) { + found = 0; + for (i = P->nvars - 1; i >= 0 && !found; i--) + if (!strcmp (P->names[i], bn[d]) && P->vals[i].n == 3) { + s->b[d][0] = P->vals[i].v[0]; + s->b[d][1] = P->vals[i].v[1]; + s->b[d][2] = P->vals[i].v[2]; + found = 1; + } + if (!found) { + fprintf (stderr, "%s: ERROR: %s (3-vector) not defined in Qgridinfo\n", comp, bn[d]); + exit (-1); + } + found = 0; + for (i = P->nvars - 1; i >= 0 && !found; i--) + if (!strcmp (P->names[i], nn[d]) && P->vals[i].n == 1) { + s->n[d] = (int)floor (P->vals[i].v[0] + 0.5); + found = 1; + } + if (!found || s->n[d] < 2) { + fprintf (stderr, "%s: ERROR: %s (>=2) not defined in Qgridinfo\n", comp, nn[d]); + exit (-1); + } + } + free (P); +} + +// x_star = (y x z)/vol etc. (mcni::get_inversions) +double +phdisp_inversions (double a[3][3], double inv[3][3]) { + int i, j; + double c[3][3]; + for (i = 0; i < 3; i++) { + int i1 = (i + 1) % 3, i2 = (i + 2) % 3; + c[i][0] = a[i1][1] * a[i2][2] - a[i1][2] * a[i2][1]; + c[i][1] = a[i1][2] * a[i2][0] - a[i1][0] * a[i2][2]; + c[i][2] = a[i1][0] * a[i2][1] - a[i1][1] * a[i2][0]; + } + double vol = a[0][0] * c[0][0] + a[0][1] * c[0][1] + a[0][2] * c[0][2]; + for (i = 0; i < 3; i++) + for (j = 0; j < 3; j++) + inv[i][j] = c[i][j] / vol; + return vol; +} + +void +phdisp_read_dispersion (const char* dir, struct phdisp_struct* s, const char* comp) { + FILE* fp; + int D, Nb, Nq, br, q; + long nE, nP, i; + + phdisp_read_qgridinfo (dir, s, comp); + if (fabs (phdisp_inversions (s->b, s->binv)) < 1e-12) { + fprintf (stderr, "%s: ERROR: b1, b2, b3 in Qgridinfo are not linearly independent\n", comp); + exit (-1); + } + + // Omega2: D, N_b (atoms), N_q, then N_q x (N_b*D) doubles + fp = phdisp_open (dir, "Omega2", "rb", comp); + phdisp_idf_header (fp, "Omega2", comp); + if (fread (&D, sizeof (int), 1, fp) != 1 || fread (&Nb, sizeof (int), 1, fp) != 1 || fread (&Nq, sizeof (int), 1, fp) != 1) { + fprintf (stderr, "%s: ERROR: corrupt Omega2 file\n", comp); + exit (-1); + } + if (D != 3 || Nq != s->n[0] * s->n[1] * s->n[2] || Nb < 1 || Nb > PHDISP_MAXATOMS) { + fprintf (stderr, "%s: ERROR: Omega2 has D=%d, N_atoms=%d, N_q=%d, but Qgridinfo gives %d x %d x %d points\n", comp, D, + Nb, Nq, s->n[0], s->n[1], s->n[2]); + exit (-1); + } + s->n_atoms = Nb; + s->n_branches = 3 * Nb; + nE = (long)Nq * s->n_branches; + s->E = malloc (nE * sizeof (double)); + if (fread (s->E, sizeof (double), nE, fp) != (size_t)nE) { + fprintf (stderr, "%s: ERROR: Omega2 file too short\n", comp); + exit (-1); + } + fclose (fp); + for (i = 0; i < nE; i++) + s->E[i] = (s->E[i] < 0 ? 0 : sqrt (s->E[i])) * PHDISP_HZ2MEV; + + s->Emin = malloc (s->n_branches * sizeof (double)); + s->Emax = malloc (s->n_branches * sizeof (double)); + for (br = 0; br < s->n_branches; br++) { + s->Emin[br] = 1e300; + s->Emax[br] = -1e300; + for (q = 0; q < Nq; q++) { + double e = s->E[(long)q * s->n_branches + br]; + if (e < s->Emin[br]) + s->Emin[br] = e; + if (e > s->Emax[br]) + s->Emax[br] = e; + } + if (br == 0 || s->Emax[br] > s->Emax_all) + s->Emax_all = s->Emax[br]; + } + + // Polarizations: D, N_b, N_q, then N_q x (N_b*D) x N_b x D x 2 doubles + fp = phdisp_open (dir, "Polarizations", "rb", comp); + phdisp_idf_header (fp, "Polarizations", comp); + int D2, Nb2, Nq2; + if (fread (&D2, sizeof (int), 1, fp) != 1 || fread (&Nb2, sizeof (int), 1, fp) != 1 || fread (&Nq2, sizeof (int), 1, fp) != 1 || D2 != D || Nb2 != Nb + || Nq2 != Nq) { + fprintf (stderr, "%s: ERROR: Polarizations file inconsistent with Omega2 file\n", comp); + exit (-1); + } + nP = (long)Nq * s->n_branches * Nb * 3 * 2; + s->pol = malloc (nP * sizeof (double)); + if (fread (s->pol, sizeof (double), nP, fp) != (size_t)nP) { + fprintf (stderr, "%s: ERROR: Polarizations file too short\n", comp); + exit (-1); + } + fclose (fp); +} + +void +phdisp_read_xyz (const char* file, double b_default, double m_default, struct phdisp_struct* s, const char* comp) { + FILE* fp = phdisp_open (NULL, file, "r", comp); + char line[4096]; + double lat[3][3], inv[3][3]; + int natoms = -1, stage = 0, ia = 0, d; + while (fgets (line, sizeof (line), fp)) { + char* c = line; + while (*c == ' ' || *c == '\t') + c++; + if (*c == '\0' || *c == '\n' || *c == '\r' || *c == '#') + continue; + if (stage == 0) { + natoms = atoi (c); + if (natoms != s->n_atoms) { + fprintf (stderr, "%s: ERROR: %s has %d atoms, but the dispersion has %d\n", comp, file, natoms, s->n_atoms); + exit (-1); + } + s->pos = calloc (3 * natoms, sizeof (double)); + s->bc = calloc (natoms, sizeof (double)); + s->m = calloc (natoms, sizeof (double)); + stage = 1; + } else if (stage == 1) { + if (sscanf (c, "%lf %lf %lf %lf %lf %lf %lf %lf %lf", &lat[0][0], &lat[0][1], &lat[0][2], &lat[1][0], &lat[1][1], &lat[1][2], &lat[2][0], &lat[2][1], + &lat[2][2]) + != 9) { + fprintf (stderr, "%s: ERROR: line 2 of %s must contain the 9 components of the lattice vectors\n", comp, file); + exit (-1); + } + stage = 2; + } else if (ia < natoms) { + char sym[64]; + double f[3], extra[3]; + int nread = sscanf (c, "%63s %lf %lf %lf %lf %lf %lf", sym, &f[0], &f[1], &f[2], &extra[0], &extra[1], &extra[2]); + if (nread < 4) { + fprintf (stderr, "%s: ERROR: could not parse atom line in %s: %s", comp, file, line); + exit (-1); + } + for (d = 0; d < 3; d++) + s->pos[3 * ia + d] = f[0] * lat[0][d] + f[1] * lat[1][d] + f[2] * lat[2][d]; + if (nread >= 6) { // Symbol x y z b_coh mass + s->bc[ia] = extra[0]; + s->m[ia] = extra[1]; + } else { + s->bc[ia] = b_default; + s->m[ia] = m_default; + } + if (s->bc[ia] < 0 || s->m[ia] <= 0) { + fprintf (stderr, + "%s: ERROR: atom %d (%s) has no scattering length / mass. Give them as columns 5 and 6 in %s or set b_coh and " + "mass\n", + comp, ia, sym, file); + exit (-1); + } + ia++; + } + } + fclose (fp); + if (ia != natoms) { + fprintf (stderr, "%s: ERROR: expected %d atoms in %s, read %d\n", comp, natoms, file, ia); + exit (-1); + } + s->uc_vol = fabs (phdisp_inversions (lat, inv)); +} + +// DOS (IDF): N_bins, dE [THz], then N_bins doubles. Returns 2W/Q^2 [AA^2] (DWFromDOS.icc) +double +phdisp_dw_core_from_dos (const char* dir, const char* file, double avg_mass, double T, const char* comp) { + FILE* fp = phdisp_open (dir, file, "rb", comp); + int nb, i, nSample = 100, first = -1; + double dE, *Z, area = 0; + phdisp_idf_header (fp, "DOS", comp); + if (fread (&nb, sizeof (int), 1, fp) != 1 || fread (&dE, sizeof (double), 1, fp) != 1 || nb < 2) { + fprintf (stderr, "%s: ERROR: corrupt DOS file\n", comp); + exit (-1); + } + Z = malloc (nb * sizeof (double)); + if (fread (Z, sizeof (double), nb, fp) != (size_t)nb) { + fprintf (stderr, "%s: ERROR: DOS file too short\n", comp); + exit (-1); + } + fclose (fp); + dE *= 2 * PI * 1e12 * PHDISP_HZ2MEV; // THz (not angular) -> meV + for (i = 0; i < nb; i++) + area += Z[i]; + area *= dE; + for (i = 0; i < nb; i++) + Z[i] /= area; + double emin = 0, emax = dE * (nb - 1); + double dw = (emax - emin) / (nSample - 1 + .00000001), core = 0; + double* f = calloc (nSample, sizeof (double)); + for (i = 0; i < nSample; i++) { + double w = dw * i + emin, z = 0; + // linearly interpolated DOS, zero outside [emin, emax) + if (w >= emin && w < emax) { + double x = (w - emin) / dE; + int j = (int)floor (x); + if (j >= nb - 1) { + j = nb - 2; + } + z = Z[j] + (x - j) * (Z[j + 1] - Z[j]); + } + if (w < emax / nSample / 100.) + continue; + if (first == -1) + first = i; + f[i] = (2 / (exp (w / (T * PHDISP_T2E)) - 1) + 1) / w * z; + core += f[i] * (i == nSample - 1 ? 0.5 : 1); + } + if (first == -1 || first + 1 >= nSample) { + fprintf (stderr, "%s: ERROR: invalid DOS (no data for E>0)\n", comp); + exit (-1); + } + double f0 = f[first] - (dw * first + emin) * (f[first + 1] - f[first]) / dw; + core += f0 / 2; + core /= PHDISP_EV * 1e-3; + core *= dw; + core *= PHDISP_HBAR * PHDISP_HBAR / 2 / PHDISP_AMU / avg_mass; + core *= 1e20; + free (f); + free (Z); + return core; +} + +// one-phonon coherent structure factor |sum_d b_d/sqrt(M_d) exp(iQ.d) (Q.e_d)|^2, normalised by the +// total coherent cross section as in MCViNE (scattering_length.icc + kernels): +// returns |..|^2 [fm^2 AA^-2 amu^-1] / 1e30 / (sigma_coh*1e-28) +double +phdisp_structure_factor (struct phdisp_struct* s, int branch, double* Q) { + double sre = 0, sim = 0; + int a; + for (a = 0; a < s->n_atoms; a++) { + double ere[3], eim[3]; + phdisp_polarization (s, branch, a, Q, ere, eim); + double epslen = sqrt (ere[0] * ere[0] + ere[1] * ere[1] + ere[2] * ere[2] + eim[0] * eim[0] + eim[1] * eim[1] + eim[2] * eim[2]); + if (epslen <= 0) + continue; + double qe_re = (Q[0] * ere[0] + Q[1] * ere[1] + Q[2] * ere[2]) / epslen; + double qe_im = (Q[0] * eim[0] + Q[1] * eim[1] + Q[2] * eim[2]) / epslen; + double qd = Q[0] * s->pos[3 * a] + Q[1] * s->pos[3 * a + 1] + Q[2] * s->pos[3 * a + 2]; + double pre = s->bc[a] / sqrt (s->m[a]); + double cr = cos (qd), ci = sin (qd); + sre += pre * (cr * qe_re - ci * qe_im); + sim += pre * (cr * qe_im + ci * qe_re); + } + return (sre * sre + sim * sim) / 1e30 / (s->sigma_coh * 1e-28); +} + +// sum of coherent cross sections and average mass of the unit cell +void +phdisp_cell_sums (struct phdisp_struct* s) { + int ia; + s->sigma_coh = 0; + s->avg_mass = 0; + for (ia = 0; ia < s->n_atoms; ia++) { + s->sigma_coh += 4 * PI * s->bc[ia] * s->bc[ia] / 100.0; // fm^2 -> barn + s->avg_mass += s->m[ia]; + } + s->avg_mass /= s->n_atoms; +} + +void +phdisp_free (struct phdisp_struct* s) { + free (s->E); + free (s->pol); + free (s->Emin); + free (s->Emax); + free (s->pos); + free (s->bc); + free (s->m); +} + +#endif // PHONON_DISPERSION_LIB_C diff --git a/mcstas-comps/share/phonon-dispersion-lib.h b/mcstas-comps/share/phonon-dispersion-lib.h new file mode 100644 index 0000000000..944c70f067 --- /dev/null +++ b/mcstas-comps/share/phonon-dispersion-lib.h @@ -0,0 +1,101 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* Library: share/phonon-dispersion-lib.h +* +* Written by: Peter Willendrup, port of MCViNE (J.Y.Y. Lin et al.) phonon kernel code +* Date: 23.09.2026 +* Origin: DTU Physics +* +* Shared code for the Union processes CoherentPhononSingleXtal_process and +* CoherentPhononPowder_process: +* - readers for phonon dispersions in the MCViNE IDF format (Qgridinfo, Omega2, +* Polarizations, DOS) and MCViNE/diffpy xyz crystal structure files +* - periodic, trilinearly interpolated phonon energies and polarizations +* (PeriodicDispersion_3D, ChangeCoordinateSystem_forDispersion_3D and +* LinearlyInterpolatedDispersionOnGrid_3D in MCViNE) +* - one-phonon coherent structure factor (scattering_length.icc) +* - Debye-Waller factor from the DOS (DWFromDOS.icc) +* - root finding for omega(Q) = |Ei-Ef| (Omega_minus_deltaE, FindRootsEvenly, ZRidd) +* +* Usage: %include "read_table-lib" and then %include "phonon-dispersion-lib" in SHARE. Error messages are prefixed +* with the string "comp" passed by the caller. +* +*******************************************************************************/ + +#ifndef PHONON_DISPERSION_LIB_H +#define PHONON_DISPERSION_LIB_H + +#include + +#define PHDISP_T2E (1.0 / 11.605) // Kelvin to meV, as in MCViNE +#define PHDISP_HBAR 1.05457148e-34 // J s +#define PHDISP_EV 1.60217653e-19 // J +#define PHDISP_AMU 1.66053886e-27 // kg +#define PHDISP_HZ2MEV (PHDISP_HBAR / (1e-3 * PHDISP_EV)) // angular frequency [rad/s] to meV +#define PHDISP_MAXATOMS 1024 + +// Phonon dispersion on a grid + unit cell content +struct phdisp_struct { + // grid dispersion + int n_atoms; + int n_branches; + int n[3]; // grid points along b1, b2, b3 + double b[3][3]; // reciprocal cell of the grid [AA^-1], b[i] is vector b_i + double binv[3][3]; // binv[i] . b[j] = delta_ij, gives fractional coordinates + double* E; // [n1][n2][n3][branch] phonon energy [meV] + double* pol; // [n1][n2][n3][branch][atom][xyz][re/im] + double* Emin; // per branch + double* Emax; // per branch + double Emax_all; // maximum phonon energy + // unit cell + double uc_vol; // [AA^3] + double* pos; // [atom][xyz] cartesian [AA] + double* bc; // coherent scattering length [fm] + double* m; // mass [amu] + double sigma_coh; // sum of coherent cross sections in unit cell [barn] + double avg_mass; // [amu] +}; + +// Context for the function omega(Q(vf)) - |Ei - Ef| +struct phdisp_omega_ctx { + struct phdisp_struct* s; + int branch; + double vf_dir[3]; + double vi[3]; + double vi_l; + int powder; // 1: Q = |ki-kf| * qhat (orientational average of a powder) + double qhat[3]; // unit vector in the crystal frame, used when powder=1 +}; + +// Qgridinfo expression evaluator +#define PHDISP_MAXVARS 64 +struct phdisp_val { + int n; + double v[3]; +}; + +struct phdisp_parser { + const char* p; + int nvars; + char names[PHDISP_MAXVARS][64]; + struct phdisp_val vals[PHDISP_MAXVARS]; + int error; +}; + +// public functions +void phdisp_read_dispersion (const char* dir, struct phdisp_struct* s, const char* comp); +void phdisp_read_xyz (const char* file, double b_default, double m_default, struct phdisp_struct* s, const char* comp); +void phdisp_cell_sums (struct phdisp_struct* s); +double phdisp_dw_core_from_dos (const char* dir, const char* file, double avg_mass, double T, const char* comp); +double phdisp_energy (struct phdisp_struct* s, int branch, double* Q); +void phdisp_polarization (struct phdisp_struct* s, int branch, int atom, double* Q, double* re, double* im); +double phdisp_structure_factor (struct phdisp_struct* s, int branch, double* Q); +double phdisp_omega_minus_dE (double vf, struct phdisp_omega_ctx* c); +int phdisp_find_roots (struct phdisp_omega_ctx* c, double x1, double x2, int nsteps, double xacc, double* roots); +double phdisp_bose (double omega, double T); +void phdisp_free (struct phdisp_struct* s); + +#endif // PHONON_DISPERSION_LIB_H diff --git a/mcstas-comps/share/tinyexpr.c b/mcstas-comps/share/tinyexpr.c index 7f2e91694c..673a041d5c 100644 --- a/mcstas-comps/share/tinyexpr.c +++ b/mcstas-comps/share/tinyexpr.c @@ -26,6 +26,9 @@ // OF THE MCSTAS SOFTWARE PACKAGE +#ifndef TINYEXPR_C +#define TINYEXPR_C + /* COMPILE TIME OPTIONS */ /* Exponentiation associativity: @@ -750,3 +753,5 @@ static void pn (const te_expr *n, int depth) { void te_print(const te_expr *n) { pn(n, 0); } + +#endif /* TINYEXPR_C */ diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 8e6e6176df..41d1225ad8 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -33,6 +33,15 @@ enum process { PhononSimple, Texture, IncoherentPhonon, + Resolution, + ElasticSQ, + IsotropicSqw, + DispersionPowder, + DispersionSingleXtal, + IncoherentElastic, + CoherentPhononPowder, + IncoherentOnePhonon, + CoherentPhononSingleXtal, NCrystal, Non, Template @@ -488,6 +497,15 @@ union data_transfer_union{ struct Single_crystal_physics_storage_struct *pointer_to_a_Single_crystal_physics_storage_struct; struct AF_HB_1D_physics_storage_struct *pointer_to_a_AF_HB_1D_physics_storage_struct; struct IncoherentPhonon_physics_storage_struct *pointer_to_a_IncoherentPhonon_physics_storage_struct; + struct Resolution_physics_storage_struct *pointer_to_a_Resolution_physics_storage_struct; + struct ElasticSQ_physics_storage_struct *pointer_to_a_ElasticSQ_physics_storage_struct; + struct IsotropicSqw_physics_storage_struct *pointer_to_a_IsotropicSqw_physics_storage_struct; + struct DispersionPowder_physics_storage_struct *pointer_to_a_DispersionPowder_physics_storage_struct; + struct DispersionSingleXtal_physics_storage_struct *pointer_to_a_DispersionSingleXtal_physics_storage_struct; + struct IncoherentElastic_physics_storage_struct *pointer_to_a_IncoherentElastic_physics_storage_struct; + struct CoherentPhononPowder_physics_storage_struct *pointer_to_a_CoherentPhononPowder_physics_storage_struct; + struct IncoherentOnePhonon_physics_storage_struct *pointer_to_a_IncoherentOnePhonon_physics_storage_struct; + struct CoherentPhononSingleXtal_physics_storage_struct *pointer_to_a_CoherentPhononSingleXtal_physics_storage_struct; struct PhononSimpleNumeric_physics_storage_struct *pointer_to_a_PhononSimpleNumeric_storage_struct; struct PhononSimple_physics_storage_struct *pointer_to_a_PhononSimple_storage_struct; struct MagnonSimple_physics_storage_struct *pointer_to_a_MagnonSimple_storage_struct; diff --git a/mcstas-comps/share/union-suffix.c b/mcstas-comps/share/union-suffix.c index f83bd9988a..3c0d1bd358 100644 --- a/mcstas-comps/share/union-suffix.c +++ b/mcstas-comps/share/union-suffix.c @@ -51,6 +51,51 @@ int physics_my(enum process choice, double *my,double *k_initial, union data_tra output = IncoherentPhonon_physics_my(my, k_initial, data_transfer, focus_data, _particle); break; #endif + #ifdef PROCESS_RESOLUTION_DETECTOR + case Resolution: + output = Resolution_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_ELASTICSQ_DETECTOR + case ElasticSQ: + output = ElasticSQ_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_ISOTROPICSQW_DETECTOR + case IsotropicSqw: + output = IsotropicSqw_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_DISPERSIONPOWDER_DETECTOR + case DispersionPowder: + output = DispersionPowder_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_DISPERSIONSINGLEXTAL_DETECTOR + case DispersionSingleXtal: + output = DispersionSingleXtal_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_INCOHERENTELASTIC_DETECTOR + case IncoherentElastic: + output = IncoherentElastic_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_COHERENTPHONONPOWDER_DETECTOR + case CoherentPhononPowder: + output = CoherentPhononPowder_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_INCOHERENTONEPHONON_DETECTOR + case IncoherentOnePhonon: + output = IncoherentOnePhonon_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_COHERENTPHONONSINGLEXTAL_DETECTOR + case CoherentPhononSingleXtal: + output = CoherentPhononSingleXtal_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif #ifdef PROCESS_NCRYSTAL_DETECTOR case NCrystal: output = NCrystal_physics_my(my, k_initial, data_transfer, focus_data, _particle); @@ -120,6 +165,51 @@ int physics_scattering(enum process choice, double *k_final, double *k_initial, output = IncoherentPhonon_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); break; #endif + #ifdef PROCESS_RESOLUTION_DETECTOR + case Resolution: + output = Resolution_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_ELASTICSQ_DETECTOR + case ElasticSQ: + output = ElasticSQ_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_ISOTROPICSQW_DETECTOR + case IsotropicSqw: + output = IsotropicSqw_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_DISPERSIONPOWDER_DETECTOR + case DispersionPowder: + output = DispersionPowder_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_DISPERSIONSINGLEXTAL_DETECTOR + case DispersionSingleXtal: + output = DispersionSingleXtal_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_INCOHERENTELASTIC_DETECTOR + case IncoherentElastic: + output = IncoherentElastic_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_COHERENTPHONONPOWDER_DETECTOR + case CoherentPhononPowder: + output = CoherentPhononPowder_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_INCOHERENTONEPHONON_DETECTOR + case IncoherentOnePhonon: + output = IncoherentOnePhonon_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif + #ifdef PROCESS_COHERENTPHONONSINGLEXTAL_DETECTOR + case CoherentPhononSingleXtal: + output = CoherentPhononSingleXtal_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif #ifdef PROCESS_NCRYSTAL_DETECTOR case NCrystal: output = NCrystal_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); diff --git a/mcstas-comps/union/CoherentPhononPowder_process.comp b/mcstas-comps/union/CoherentPhononPowder_process.comp new file mode 100644 index 0000000000..e01dc2ca0e --- /dev/null +++ b/mcstas-comps/union/CoherentPhononPowder_process.comp @@ -0,0 +1,408 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernel by Jiao Lin +* Date: 24.09.2026 +* Origin: DTU Physics +* +* Coherent inelastic one-phonon scattering in a powder (polycrystal), using a phonon +* dispersion (energies and polarization vectors) given on a grid in reciprocal space. +* +* %D +* Port of the MCViNE kernel +* mccomponents/lib/kernels/sample/phonon/CoherentInelastic_PolyXtal +* (J.Y.Y. Lin et al., Caltech/ORNL) to a Union process. +* +* The scattering is the orientational average of the single crystal one-phonon +* cross section of CoherentPhononSingleXtal_process, +* 1/sigma_coh d2sigma/dOmega dEf = kf/ki sum_s < A_s(RQ) delta(|Ei-Ef| - omega_s(RQ)) >_R +* where A_s contains the normalised one-phonon structure factor, 1/|omega|, +* the Debye-Waller and the Bose factor. Dispersion data (MCViNE IDF directory) +* and crystal structure (xyz file) are given exactly as for +* CoherentPhononSingleXtal_process, see that component for the file formats. +* +* Two sampling methods are available: +* +* method=0 (as MCViNE): a phonon branch and a phonon wave vector q (uniformly in a +* sphere of radius Qmax = ki + k(Ei+Emax)) are picked. The phonon energy fixes Ef +* (energy loss or gain), |q| fixes the scattering angle and the azimuthal angle is +* random. Union focusing is not used. Samples that are kinematically impossible +* are absorbed (weight 0) and the weight contains the exact sampled volume. +* MCViNE instead retries until a valid q is found and uses an approximate +* "accessible reciprocal volume", which only approximately normalises the result. +* method=1 (focusing): the final direction is picked using the Union focusing of the +* geometry, a branch and an orientation (direction of q in the crystal) are picked +* at random and the final speed is found by root finding as in +* CoherentPhononSingleXtal_process. Use this with focusing onto detectors. +* +* Both methods give the same result (without focusing), method=0 is faster for 4pi +* scattering, method=1 is much more efficient when focusing on a small detector. +* +* Note on combining processes: as in MCViNE, the inverse penetration depth of this +* process is the full cross section sigma/V_uc, and the kernel physics is carried by +* the weight. When several processes describing parts of the same cross section are +* combined in one material, e.g. +* Powder_process + CoherentPhononPowder_process, +* their attenuations add up, so the attenuation is overestimated for thick samples. +* Single scattering intensities are not affected. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* Example: CoherentPhononPowder_process(dispersion_dir="fccNi-phonons", +* xyz_file="Ni.xyz", T=300, method=1) +* +* %P +* INPUT PARAMETERS: +* dispersion_dir: [string] Directory with the IDF dispersion files (Qgridinfo, Omega2, Polarizations, DOS) +* xyz_file: [string] xyz file with lattice vectors and atoms (fractional coordinates) of the unit cell +* b_coh: [fm] Coherent scattering length used for atoms without explicit value in xyz_file +* mass: [amu] Atomic mass used for atoms without explicit value in xyz_file +* T: [K] Temperature +* DW_core: [AA^2] Debye-Waller core, 2W = DW_core*Q^2. Negative: calculate from the DOS file in dispersion_dir +* DOS_file: [string] Optional DOS file (IDF format) overriding dispersion_dir/DOS +* method: [1] 0: MCViNE sampling of q (no focusing), 1: focusing + root finding +* min_omega: [meV] Phonons with energy below min_omega are ignored (singular 1/omega near Gamma) +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* branch_Emin_factor: [1] If positive, only branches with minimum energy below branch_Emin_factor*Ei are sampled +* root_nsteps: [1] method=1: Number of intervals in which roots of omega(Q)-dE are searched for +* root_xacc: [m/s] method=1: Accuracy of the final neutron speed from the root finder +* deltaV_Jacobi: [1] method=1: Relative velocity step used for the numerical Jacobian +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* CPP_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* %L +* J.Y.Y. Lin et al., Nucl. Instrum. Methods A 810 (2016) 86 +* %L +* G.L. Squires, Introduction to the Theory of Thermal Neutron Scattering, ch. 3 +* +* %E +******************************************************************************/ + +DEFINE COMPONENT CoherentPhononPowder_process + +SETTING PARAMETERS(string dispersion_dir="", string xyz_file="", double b_coh=-1, double mass=-1, + double T=300, double DW_core=-1, string DOS_file="", int method=1, double min_omega=0.01, + double packing_factor=1, double branch_Emin_factor=-1, int root_nsteps=100, double root_xacc=0.1, + double deltaV_Jacobi=0.001, double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this CoherentPhononPowder_process component" + #endif + + %include "read_table-lib" + %include "phonon-dispersion-lib" + + // Very important to add a pointer to this struct in the union-lib.c file + struct CoherentPhononPowder_physics_storage_struct { + struct phdisp_struct disp; // dispersion and unit cell + // physics + double Temp; + double dw_core; + double my_scattering; + double w_min; + // numerics + int meth; + double emin_factor; + int nsteps; + double xacc; + double deltaV; + int* good_branches; // work array + double* roots; // work array + }; + + // Inverse penetration depth: constant sum(sigma_coh)/V_uc as in MCViNE + int + CoherentPhononPowder_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_CoherentPhononPowder_physics_storage_struct->my_scattering; + return 1; + } + + // A_s(Q) = ksquare2E(normalised |F|^2)/|omega| * DW * Bose, Q in the crystal frame + double + cpp_A (struct CoherentPhononPowder_physics_storage_struct* s, int branch, double* Q, double omega) { + double Q2 = Q[0] * Q[0] + Q[1] * Q[1] + Q[2] * Q[2]; + double norm_of_slsum = phdisp_structure_factor (&s->disp, branch, Q); + return VS2E * K2V * K2V * norm_of_slsum / fabs (omega) * exp (-s->dw_core * Q2) * phdisp_bose (omega, s->Temp); + } + + int + CoherentPhononPowder_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + struct CoherentPhononPowder_physics_storage_struct* s = data_transfer.pointer_to_a_CoherentPhononPowder_physics_storage_struct; + int d, br, ngood = 0, index; + double vi[3], vi_l, Ei, ki; + + for (d = 0; d < 3; d++) + vi[d] = K2V * k_initial[d]; + vi_l = sqrt (vi[0] * vi[0] + vi[1] * vi[1] + vi[2] * vi[2]); + Ei = VS2E * vi_l * vi_l; + ki = V2K * vi_l; + + // branches to sample + for (br = 0; br < s->disp.n_branches; br++) + if (s->emin_factor <= 0 || s->disp.Emin[br] < Ei * s->emin_factor) + s->good_branches[ngood++] = br; + if (ngood == 0) + return 0; + index = (int)floor (rand01 () * ngood); + if (index >= ngood) + index = ngood - 1; + br = s->good_branches[index]; + + if (s->meth == 0) { + // --- MCViNE CoherentInelastic_PolyXtal: sample q uniformly in a sphere --- + double Qmax = ki + V2K * sqrt ((Ei + s->disp.Emax_all) / VS2E); + double q[3], ql, omega, Ef, vf_l, vQ, nsign = 1; + do { + q[0] = randpm1 (); + q[1] = randpm1 (); + q[2] = randpm1 (); + } while (q[0] * q[0] + q[1] * q[1] + q[2] * q[2] > 1); + for (d = 0; d < 3; d++) + q[d] *= Qmax; + ql = sqrt (q[0] * q[0] + q[1] * q[1] + q[2] * q[2]); + double pw = phdisp_energy (&s->disp, br, q); + if (pw < s->w_min) + return 0; + // energy loss or gain (pick_Ef) + if (pw < Ei) { + nsign = 2; + Ef = rand01 () >= 0.5 ? Ei + pw : Ei - pw; + } else + Ef = Ei + pw; + omega = Ei - Ef; + vf_l = sqrt (Ef / VS2E); + vQ = K2V * ql; + if (vQ < fabs (vi_l - vf_l) || vQ > vi_l + vf_l) + return 0; + // final direction: angle to ki from |q|, random azimuth (pick_v_f) + double cos_theta = (vi_l * vi_l + vf_l * vf_l - vQ * vQ) / (2 * vi_l * vf_l); + if (cos_theta > 1) + cos_theta = 1; + if (cos_theta < -1) + cos_theta = -1; + double sin_theta = sqrt (1 - cos_theta * cos_theta), phi = 2 * PI * rand01 (); + double e1[3] = { vi[0] / vi_l, vi[1] / vi_l, vi[2] / vi_l }, e2[3], e3[3], n; + if (fabs (e1[0]) > 1e-10 || fabs (e1[1]) > 1e-10) { // e2 = z x e1 + e2[0] = -e1[1]; + e2[1] = e1[0]; + e2[2] = 0; + } else { + e2[0] = 1; + e2[1] = 0; + e2[2] = 0; + } + n = sqrt (e2[0] * e2[0] + e2[1] * e2[1] + e2[2] * e2[2]); + for (d = 0; d < 3; d++) + e2[d] /= n; + e3[0] = e1[1] * e2[2] - e1[2] * e2[1]; + e3[1] = e1[2] * e2[0] - e1[0] * e2[2]; + e3[2] = e1[0] * e2[1] - e1[1] * e2[0]; + for (d = 0; d < 3; d++) + k_final[d] = V2K * vf_l * (sin_theta * cos (phi) * e2[d] + sin_theta * sin (phi) * e3[d] + cos_theta * e1[d]); + // weight: V_sphere * n_branches * n_sign * A / (2 ki^2 |q|) + double Vs = 4.0 / 3.0 * PI * Qmax * Qmax * Qmax; + *weight *= Vs * ngood * nsign * cpp_A (s, br, q, omega) / (2 * ki * ki * ql); + return 1; + } + + // --- method 1: focused final direction, random crystal orientation, root finding --- + struct phdisp_omega_ctx ctx; + double solid_angle, u, cth, sth, ph; + Coords dir; + int nf; + + ctx.s = &s->disp; + ctx.branch = br; + ctx.powder = 1; + ctx.vi_l = vi_l; + for (d = 0; d < 3; d++) + ctx.vi[d] = vi[d]; + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + ctx.vf_dir[0] = dir.x; + ctx.vf_dir[1] = dir.y; + ctx.vf_dir[2] = dir.z; + // direction of Q in the crystal, uniform on the unit sphere (orientational average) + u = randpm1 (); + cth = u; + sth = sqrt (1 - u * u); + ph = 2 * PI * rand01 (); + ctx.qhat[0] = sth * cos (ph); + ctx.qhat[1] = sth * sin (ph); + ctx.qhat[2] = cth; + + nf = phdisp_find_roots (&ctx, 0, 2 * vi_l, s->nsteps, s->xacc, s->roots); + if (nf < 1) + return 0; + index = (int)floor (rand01 () * nf); + if (index >= nf) + index = nf - 1; + double vf_l = s->roots[index]; + if (vf_l <= 0) + return 0; + double delta_v = s->deltaV * vi_l; + double f1 = phdisp_omega_minus_dE (vf_l - delta_v, &ctx); + double f2 = phdisp_omega_minus_dE (vf_l + delta_v, &ctx); + double Jacobi = fabs (f2 - f1) / (2 * delta_v); + if (Jacobi <= 0) + return 0; + double omega = Ei - VS2E * vf_l * vf_l; + if (fabs (omega) < s->w_min) + return 0; + double Q[3], ql = 0; + for (d = 0; d < 3; d++) { + k_final[d] = V2K * vf_l * ctx.vf_dir[d]; + ql += (k_initial[d] - k_final[d]) * (k_initial[d] - k_final[d]); + } + ql = sqrt (ql); + for (d = 0; d < 3; d++) + Q[d] = ql * ctx.qhat[d]; + *weight *= solid_angle * ngood * nf * (2 * VS2E * vf_l / Jacobi) * cpp_A (s, br, Q, omega) * vf_l / vi_l; + return 1; + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_COHERENTPHONONPOWDER_DETECTOR + #define PROCESS_COHERENTPHONONPOWDER_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct CoherentPhononPowder_physics_storage_struct CPP_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + char who[512]; + snprintf (who, 512, "CoherentPhononPowder_process:%s", NAME_CURRENT_COMP); + + memset (&CPP_storage, 0, sizeof (CPP_storage)); + + if (!strlen (dispersion_dir) || !strlen (xyz_file)) { + fprintf (stderr, "%s: ERROR: dispersion_dir and xyz_file must be given\n", who); + exit (-1); + } + if (T <= 0) { + fprintf (stderr, "%s: ERROR: T must be positive\n", who); + exit (-1); + } + if (method != 0 && method != 1) { + fprintf (stderr, "%s: ERROR: method must be 0 or 1\n", who); + exit (-1); + } + if (root_nsteps < 1 || root_xacc <= 0 || deltaV_Jacobi <= 0) { + fprintf (stderr, "%s: ERROR: root_nsteps, root_xacc and deltaV_Jacobi must be positive\n", who); + exit (-1); + } + + phdisp_read_dispersion (dispersion_dir, &CPP_storage.disp, who); + phdisp_read_xyz (xyz_file, b_coh, mass, &CPP_storage.disp, who); + phdisp_cell_sums (&CPP_storage.disp); + if (CPP_storage.disp.sigma_coh <= 0) { + fprintf (stderr, "%s: ERROR: total coherent cross section is zero\n", who); + exit (-1); + } + + CPP_storage.Temp = T; + if (DW_core >= 0) + CPP_storage.dw_core = DW_core; + else if (strlen (DOS_file)) + CPP_storage.dw_core = phdisp_dw_core_from_dos (NULL, DOS_file, CPP_storage.disp.avg_mass, T, who); + else + CPP_storage.dw_core = phdisp_dw_core_from_dos (dispersion_dir, "DOS", CPP_storage.disp.avg_mass, T, who); + + CPP_storage.my_scattering = packing_factor * CPP_storage.disp.sigma_coh / CPP_storage.disp.uc_vol * 100; // barn/AA^3 -> 1/m + CPP_storage.w_min = min_omega; + CPP_storage.meth = method; + CPP_storage.emin_factor = branch_Emin_factor; + CPP_storage.nsteps = root_nsteps; + CPP_storage.xacc = root_xacc; + CPP_storage.deltaV = deltaV_Jacobi; + CPP_storage.good_branches = malloc (CPP_storage.disp.n_branches * sizeof (int)); + CPP_storage.roots = malloc ((root_nsteps + 1) * sizeof (double)); + + MPI_MASTER (printf ("%s: %d atoms, %d branches, %dx%dx%d grid, Emax=%g meV, V_uc=%g AA^3, sigma_coh=%g barn, DW_core=%g AA^2, my_scattering=%g 1/m, method=%d\n", + who, CPP_storage.disp.n_atoms, CPP_storage.disp.n_branches, CPP_storage.disp.n[0], CPP_storage.disp.n[1], CPP_storage.disp.n[2], + CPP_storage.disp.Emax_all, CPP_storage.disp.uc_vol, CPP_storage.disp.sigma_coh, CPP_storage.dw_core, CPP_storage.my_scattering, + method);); + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Isotropic process (powder) + This_process.non_isotropic_rot_index = -1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = CoherentPhononPowder; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_CoherentPhononPowder_physics_storage_struct = &CPP_storage; + This_process.probability_for_scattering_function = &CoherentPhononPowder_physics_my; + This_process.scattering_function = &CoherentPhononPowder_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "%s: Error identifying Union_init component, %s is not a known component name.\n", who, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +FINALLY +%{ + phdisp_free (&CPP_storage.disp); + free (CPP_storage.good_branches); + free (CPP_storage.roots); +%} + +END diff --git a/mcstas-comps/union/CoherentPhononSingleXtal_process.comp b/mcstas-comps/union/CoherentPhononSingleXtal_process.comp new file mode 100644 index 0000000000..ba0d0c8a8e --- /dev/null +++ b/mcstas-comps/union/CoherentPhononSingleXtal_process.comp @@ -0,0 +1,356 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernel by Jiao Lin +* Date: 23.09.2026 +* Origin: DTU Physics +* +* Coherent inelastic one-phonon scattering in a single crystal, using a phonon +* dispersion (energies and polarization vectors) given on a grid in reciprocal space. +* +* %D +* Port of the MCViNE kernel +* mccomponents/lib/kernels/sample/phonon/CoherentInelastic_SingleXtal +* (J.Y.Y. Lin et al., Caltech/ORNL) to a Union process. +* +* The phonon dispersion is read from a directory in the MCViNE "IDF" format, +* containing the files +* Qgridinfo Text file defining b1, b2, b3 [AA^-1] (the reciprocal cell spanned +* by the grid) and n1, n2, n3 (grid points along each of them). +* Simple Python-like expressions are allowed, e.g. +* import math +* twopi = 2*math.pi +* b = 0.284*twopi +* b1 = (b, -b, b) +* n1 = n2 = n3 = 21 +* Omega2 Binary file with squared phonon angular frequencies [rad^2/s^2] +* Polarizations Binary file with (complex) phonon polarization vectors +* DOS Binary file with the phonon DOS (only needed when the Debye-Waller +* factor is calculated from the DOS, i.e. when DW_core < 0) +* The grid covers the fractional coordinates [0,1] (both ends included) along +* b1, b2 and b3, and the dispersion is assumed periodic with that cell. +* +* The crystal structure is read from an MCViNE/diffpy style xyz file: +* line 1: number of atoms in the unit cell +* line 2: the 9 components of the lattice vectors a1 a2 a3 [AA] +* line 3+: Symbol x y z [b_coh mass] +* with x y z being fractional coordinates. The optional columns give the coherent +* scattering length [fm] and the mass [amu] of that atom. When they are left out, +* the component parameters b_coh and mass are used. The order of the atoms must +* match the atom index of the Polarizations file. +* +* Lattice vectors, reciprocal grid vectors and polarization vectors are all given +* in the local coordinate system of this process component, which can be rotated +* using ROTATED to orient the crystal. +* +* Algorithm (as in MCViNE): +* The inverse penetration depth is the constant sum(sigma_coh)/V_uc. At each +* scattering event a final direction is picked (respecting Union focusing), a +* phonon branch is picked among those with a minimum energy below +* branch_Emin_factor*Ei, and the final neutron speeds fulfilling both energy and +* momentum conservation are found by root finding. One of them is picked at +* random and the ray weight is multiplied by the one-phonon structure factor, +* Jacobian, Debye-Waller factor and Bose factor. If no solution is found the +* ray is absorbed (weight 0) which keeps the estimate unbiased; MCViNE instead +* retries up to 100 times, which slightly overestimates the intensity. +* +* Note on combining processes: as in MCViNE, the inverse penetration depth of this +* process is the full cross section sigma/V_uc, and the kernel physics is carried by +* the weight. When several processes describing parts of the same cross section are +* combined in one material, e.g. +* Single_crystal_process + CoherentPhononSingleXtal_process, +* their attenuations add up, so the attenuation is overestimated for thick samples. +* Single scattering intensities are not affected. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* Example: CoherentPhononSingleXtal_process(dispersion_dir="fccNi-phonons", +* xyz_file="Ni.xyz", b_coh=10.3, mass=58.6934, T=300) +* +* %P +* INPUT PARAMETERS: +* dispersion_dir: [string] Directory with the IDF dispersion files (Qgridinfo, Omega2, Polarizations, DOS) +* xyz_file: [string] xyz file with lattice vectors and atoms (fractional coordinates) of the unit cell +* b_coh: [fm] Coherent scattering length used for atoms without explicit value in xyz_file +* mass: [amu] Atomic mass used for atoms without explicit value in xyz_file +* T: [K] Temperature +* DW_core: [AA^2] Debye-Waller core, 2W = DW_core*Q^2. Negative: calculate from the DOS file in dispersion_dir +* DOS_file: [string] Optional DOS file (IDF format) overriding dispersion_dir/DOS +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* branch_Emin_factor: [1] Only branches with minimum energy below branch_Emin_factor*Ei are sampled +* root_nsteps: [1] Number of intervals in which roots of omega(Q)-dE are searched for +* root_xacc: [m/s] Accuracy of the final neutron speed from the root finder +* deltaV_Jacobi: [1] Relative velocity step used for the numerical Jacobian +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* CPSX_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* %L +* J.Y.Y. Lin et al., Nucl. Instrum. Methods A 810 (2016) 86 +* %L +* G.L. Squires, Introduction to the Theory of Thermal Neutron Scattering, ch. 3 +* +* %E +******************************************************************************/ + +DEFINE COMPONENT CoherentPhononSingleXtal_process + +SETTING PARAMETERS(string dispersion_dir="", string xyz_file="", double b_coh=-1, double mass=-1, + double T=300, double DW_core=-1, string DOS_file="", double packing_factor=1, + double branch_Emin_factor=1.5, int root_nsteps=100, double root_xacc=0.1, + double deltaV_Jacobi=0.001, double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this CoherentPhononSingleXtal_process component" + #endif + + %include "read_table-lib" + %include "phonon-dispersion-lib" + + // Very important to add a pointer to this struct in the union-lib.c file + struct CoherentPhononSingleXtal_physics_storage_struct { + struct phdisp_struct disp; // dispersion and unit cell + // physics + double Temp; + double dw_core; + double my_scattering; + // numerics + double emin_factor; + int nsteps; + double xacc; + double deltaV; + int* good_branches; // work array + double* roots; // work array + }; + + // Inverse penetration depth: constant sum(sigma_coh)/V_uc as in MCViNE + int + CoherentPhononSingleXtal_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_CoherentPhononSingleXtal_physics_storage_struct->my_scattering; + return 1; + } + + // CoherentInelastic_SingleXtal::pick_a_final_state + int + CoherentPhononSingleXtal_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + struct CoherentPhononSingleXtal_physics_storage_struct* s = data_transfer.pointer_to_a_CoherentPhononSingleXtal_physics_storage_struct; + struct phdisp_omega_ctx ctx; + int d, br, ngood = 0, nf, index; + double vi_l, Ei, solid_angle; + Coords dir; + + ctx.s = &s->disp; + ctx.powder = 0; + for (d = 0; d < 3; d++) + ctx.vi[d] = K2V * k_initial[d]; + vi_l = sqrt (ctx.vi[0] * ctx.vi[0] + ctx.vi[1] * ctx.vi[1] + ctx.vi[2] * ctx.vi[2]); + ctx.vi_l = vi_l; + Ei = VS2E * vi_l * vi_l; + + // branches not too high compared to Ei + for (br = 0; br < s->disp.n_branches; br++) + if (s->disp.Emin[br] < Ei * s->emin_factor) + s->good_branches[ngood++] = br; + if (ngood == 0) + return 0; + + // final direction + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + ctx.vf_dir[0] = dir.x; + ctx.vf_dir[1] = dir.y; + ctx.vf_dir[2] = dir.z; + + // branch + index = (int)floor (rand01 () * ngood); + if (index >= ngood) + index = ngood - 1; + ctx.branch = s->good_branches[index]; + + // final speed from energy and momentum conservation + nf = phdisp_find_roots (&ctx, 0, 2 * vi_l, s->nsteps, s->xacc, s->roots); + if (nf < 1) + return 0; + index = (int)floor (rand01 () * nf); + if (index >= nf) + index = nf - 1; + double vf_l = s->roots[index]; + if (vf_l <= 0) + return 0; + + // Jacobian of the energy delta function, d(omega - dE)/dvf + double delta_v = s->deltaV * vi_l; + double f1 = phdisp_omega_minus_dE (vf_l - delta_v, &ctx); + double f2 = phdisp_omega_minus_dE (vf_l + delta_v, &ctx); + double Jacobi = fabs (f2 - f1) / (2 * delta_v); + if (Jacobi <= 0) + return 0; + + double Ef = VS2E * vf_l * vf_l; + double omega = Ei - Ef; + if (omega == 0) + return 0; + + double Q[3], Q2 = 0; + for (d = 0; d < 3; d++) { + Q[d] = V2K * (ctx.vi[d] - vf_l * ctx.vf_dir[d]); + Q2 += Q[d] * Q[d]; + } + + // |sum_d b_d/sqrt(M_d) exp(iQ.d) (Q.e_d)|^2 / sigma_coh (scattering_length.icc) + double norm_of_slsum = phdisp_structure_factor (&s->disp, ctx.branch, Q); + // ksquare2E of the normalised |..|^2 over |omega| + double sf = VS2E * K2V * K2V * norm_of_slsum / fabs (omega); + + double DW = exp (-s->dw_core * Q2); + double kf_over_ki = vf_l / vi_l; + + *weight *= solid_angle * ngood * nf * (2 * VS2E * vf_l / Jacobi) * sf * kf_over_ki * DW * phdisp_bose (omega, s->Temp); + + k_final[0] = V2K * vf_l * ctx.vf_dir[0]; + k_final[1] = V2K * vf_l * ctx.vf_dir[1]; + k_final[2] = V2K * vf_l * ctx.vf_dir[2]; + return 1; + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_COHERENTPHONONSINGLEXTAL_DETECTOR + #define PROCESS_COHERENTPHONONSINGLEXTAL_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct CoherentPhononSingleXtal_physics_storage_struct CPSX_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + memset (&CPSX_storage, 0, sizeof (CPSX_storage)); + + if (!strlen (dispersion_dir)) { + fprintf (stderr, "CoherentPhononSingleXtal_process:%s: ERROR: dispersion_dir must be given\n", NAME_CURRENT_COMP); + exit (-1); + } + if (!strlen (xyz_file)) { + fprintf (stderr, "CoherentPhononSingleXtal_process:%s: ERROR: xyz_file must be given\n", NAME_CURRENT_COMP); + exit (-1); + } + if (T <= 0) { + fprintf (stderr, "CoherentPhononSingleXtal_process:%s: ERROR: T must be positive\n", NAME_CURRENT_COMP); + exit (-1); + } + if (root_nsteps < 1 || root_xacc <= 0 || deltaV_Jacobi <= 0) { + fprintf (stderr, "CoherentPhononSingleXtal_process:%s: ERROR: root_nsteps, root_xacc and deltaV_Jacobi must be positive\n", NAME_CURRENT_COMP); + exit (-1); + } + + char who[512]; + snprintf (who, 512, "CoherentPhononSingleXtal_process:%s", NAME_CURRENT_COMP); + phdisp_read_dispersion (dispersion_dir, &CPSX_storage.disp, who); + phdisp_read_xyz (xyz_file, b_coh, mass, &CPSX_storage.disp, who); + phdisp_cell_sums (&CPSX_storage.disp); + if (CPSX_storage.disp.sigma_coh <= 0) { + fprintf (stderr, "%s: ERROR: total coherent cross section is zero\n", who); + exit (-1); + } + + CPSX_storage.Temp = T; + if (DW_core >= 0) + CPSX_storage.dw_core = DW_core; + else if (strlen (DOS_file)) + CPSX_storage.dw_core = phdisp_dw_core_from_dos (NULL, DOS_file, CPSX_storage.disp.avg_mass, T, who); + else + CPSX_storage.dw_core = phdisp_dw_core_from_dos (dispersion_dir, "DOS", CPSX_storage.disp.avg_mass, T, who); + + CPSX_storage.my_scattering = packing_factor * CPSX_storage.disp.sigma_coh / CPSX_storage.disp.uc_vol * 100; // barn/AA^3 -> 1/m + CPSX_storage.emin_factor = branch_Emin_factor; + CPSX_storage.nsteps = root_nsteps; + CPSX_storage.xacc = root_xacc; + CPSX_storage.deltaV = deltaV_Jacobi; + CPSX_storage.good_branches = malloc (CPSX_storage.disp.n_branches * sizeof (int)); + CPSX_storage.roots = malloc ((root_nsteps + 1) * sizeof (double)); + + MPI_MASTER (printf ("CoherentPhononSingleXtal_process:%s: %d atoms, %d branches, %dx%dx%d grid, V_uc=%g AA^3, sigma_coh=%g barn, " + "DW_core=%g AA^2, my_scattering=%g 1/m\n", + NAME_CURRENT_COMP, CPSX_storage.disp.n_atoms, CPSX_storage.disp.n_branches, CPSX_storage.disp.n[0], CPSX_storage.disp.n[1], CPSX_storage.disp.n[2], + CPSX_storage.disp.uc_vol, CPSX_storage.disp.sigma_coh, CPSX_storage.dw_core, CPSX_storage.my_scattering);); + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Single crystal: the process has its own orientation + This_process.non_isotropic_rot_index = 1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = CoherentPhononSingleXtal; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_CoherentPhononSingleXtal_physics_storage_struct = &CPSX_storage; + This_process.probability_for_scattering_function = &CoherentPhononSingleXtal_physics_my; + This_process.scattering_function = &CoherentPhononSingleXtal_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "CoherentPhononSingleXtal_process:%s: Error identifying Union_init component, %s is not a known component name.\n", NAME_CURRENT_COMP, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +FINALLY +%{ + phdisp_free (&CPSX_storage.disp); + free (CPSX_storage.good_branches); + free (CPSX_storage.roots); +%} + +END diff --git a/mcstas-comps/union/DispersionPowder_process.comp b/mcstas-comps/union/DispersionPowder_process.comp new file mode 100644 index 0000000000..2468ad44d7 --- /dev/null +++ b/mcstas-comps/union/DispersionPowder_process.comp @@ -0,0 +1,496 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernels by Jiao Lin +* Date: 24.09.2026 +* Origin: DTU Physics +* +* Powder excitation with an analytical dispersion E(|Q|), optionally broadened +* +* %D +* Port of the MCViNE kernels mccomponents/lib/kernels/sample/E_Q_Kernel, +* Broadened_E_Q_Kernel and LorentzianBroadened_E_Q_Kernel (J.Y.Y. Lin et al., +* Caltech/ORNL) to a Union process. +* +* Scattering from an excitation in a powder or isotropic system, +* 1/sigma d2sigma/dOmega dEf = kf/ki S(Q) R(E - E(Q)) / (4 pi) +* with E = Ei-Ef the energy transfer, Q = |ki-kf| and R either a delta function +* (broadening=0), a Gaussian with standard deviation W(Q) (broadening=1) or a +* Lorentzian with half width at half maximum W(Q) (broadening=2). +* E(Q) [meV], S(Q) [1] and W(Q) [meV] are expressions (tinyexpr syntax, see +* https://github.com/codeplea/tinyexpr) +* in the variables Q [AA^-1], E [meV] (energy transfer, only in S_Q) and the free +* parameters p1..p4, e.g. E_Q="30+5*sin(Q)", S_Q="1", W_Q="0.5+0.1*Q". +* If T > 0 the excitation is also created with energy gain, E = -E(Q), and the +* intensities are multiplied by the Bose factors n(E)+1 and n(E), respectively. For a +* broadened excitation with T > 0 the line shape is taken as that of a damped mode, +* S(Q) [n(E)+1] (R(E-E(Q)) - R(E+E(Q))), +* which fulfils detailed balance and stays finite at E=0. +* +* Two sampling methods are available: +* method=0 (as MCViNE): Q is sampled uniformly in [Qmin, Qmax], E = E(Q) (+ a random +* broadening), the scattering angle follows from Q and the azimuthal angle is random. +* Union focusing is not used. Kinematically impossible samples are absorbed +* (MCViNE retries up to 100 times, which biases the normalisation). +* method=1 (focusing): the final direction is sampled with the Union focusing of the +* geometry. Without broadening, kf is found by root finding of E(Q(kf)) - (Ei-Ef); +* with broadening, Ef is sampled uniformly within Ei-Emax..Ei+Emax (Ef > 0). +* Only Q values within [Qmin, Qmax] contribute. +* +* The inverse penetration depth is sigma/unit_cell_volume. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* %P +* INPUT PARAMETERS: +* E_Q: [string] Expression for the dispersion E(Q) [meV] +* S_Q: [string] Expression for the intensity S(Q,E) [1] +* W_Q: [string] Expression for the broadening W(Q) [meV] (Gaussian sigma or Lorentzian HWHM) +* broadening: [1] 0: none, 1: Gaussian, 2: Lorentzian +* sigma: [barns] Scattering cross section of the unit cell, defines the inverse penetration depth +* unit_cell_volume: [AA^3] Unit cell volume +* Qmin: [AA^-1] Minimum Q of the excitation +* Qmax: [AA^-1] Maximum Q of the excitation +* Emax: [meV] method=1: maximum energy transfer considered, |E| < Emax +* T: [K] Temperature. 0: E(Q) only (no Bose factor). >0: loss and gain with Bose factors +* method: [1] 0: MCViNE sampling in Q (no focusing), 1: focusing +* p1: [1] Free parameter available in the expressions +* p2: [1] Free parameter available in the expressions +* p3: [1] Free parameter available in the expressions +* p4: [1] Free parameter available in the expressions +* root_nsteps: [1] method=1: Number of intervals in kf in which roots are searched for +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* DPW_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* +* %E +******************************************************************************/ + +DEFINE COMPONENT DispersionPowder_process + +SETTING PARAMETERS(string E_Q="", string S_Q="1", string W_Q="0", int broadening=0, double sigma=1, double unit_cell_volume=1, + double Qmin=0, double Qmax=10, double Emax=100, double T=0, int method=1, + double p1=0, double p2=0, double p3=0, double p4=0, + int root_nsteps=200, double packing_factor=1, double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this DispersionPowder_process component" + #endif + + %include "tinyexpr.h" + %include "tinyexpr.c" + + #ifndef DISPERSION_EXPR_TOOLS + #define DISPERSION_EXPR_TOOLS + #define DISP_T2E (1.0 / 11.605) // Kelvin to meV + + // variables available in the expressions + struct disp_expr_vars { + double Qx, Qy, Qz, Q, h, k, l, E, p[4]; + }; + + // compile an expression with the standard variable set + te_expr* + disp_compile (const char* expr, struct disp_expr_vars* v, const char* what, const char* comp) { + int err = 0; + te_variable vars[] = { { "Qx", &v->Qx }, { "Qy", &v->Qy }, { "Qz", &v->Qz }, { "Q", &v->Q }, { "h", &v->h }, + { "k", &v->k }, { "l", &v->l }, { "E", &v->E }, { "p1", &v->p[0] }, { "p2", &v->p[1] }, + { "p3", &v->p[2] }, { "p4", &v->p[3] } }; + te_expr* e = te_compile (expr, vars, 12, &err); + if (!e) { + fprintf (stderr, "%s: ERROR: could not parse %s expression \"%s\" near character %d\n", comp, what, expr, err); + exit (-1); + } + return e; + } + + // Bose factor: n+1 for energy loss (E>0), n for energy gain + double + disp_bose (double E, double T) { + double n = 1.0 / (exp (fabs (E) / (T * DISP_T2E)) - 1.0); + return E > 0 ? 1.0 + n : n; + } + #endif + + // Very important to add a pointer to this struct in the union-lib.c file + struct DispersionPowder_physics_storage_struct { + struct disp_expr_vars v; + te_expr* E_expr; + te_expr* S_expr; + te_expr* W_expr; + int broad, meth, nsteps; + double qmin, qmax, emax, Temp; + double my_scattering; + double* roots; + }; + + double + dpw_eval (struct DispersionPowder_physics_storage_struct* s, te_expr* e, double Q, double E) { + s->v.Q = Q; + s->v.E = E; + return te_eval (e); + } + + // context for f(kf) = sign*E(|ki-kf*dir|) - (Ei - Ef) + struct dpw_ctx { + struct DispersionPowder_physics_storage_struct* s; + double ki[3], Ei, dir[3]; + int sign; + }; + + double + dpw_f (double kf, struct dpw_ctx* c) { + double Q2 = 0, q; + int d; + for (d = 0; d < 3; d++) { + q = c->ki[d] - kf * c->dir[d]; + Q2 += q * q; + } + return c->sign * dpw_eval (c->s, c->s->E_expr, sqrt (Q2), 0) - (c->Ei - VS2E * K2V * K2V * kf * kf); + } + + // Ridder's method with bracketing values fl, fh. Returns 1 on success + int + dpw_ridder (struct dpw_ctx* c, double xl, double xh, double fl, double fh, double xacc, double* root) { + double ans = xl, xm, fm, s, xnew, fnew; + int j; + if (fl == 0) { + *root = xl; + return 1; + } + if (fh == 0) { + *root = xh; + return 1; + } + if (fl * fh > 0) + return 0; + for (j = 0; j < 60; j++) { + xm = 0.5 * (xl + xh); + fm = dpw_f (xm, c); + s = sqrt (fm * fm - fl * fh); + if (s == 0.0) { + *root = xm; + return 1; + } + xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); + if (j && fabs (xnew - ans) <= xacc) { + *root = xnew; + return 1; + } + ans = xnew; + fnew = dpw_f (ans, c); + if (fnew == 0.0) { + *root = ans; + return 1; + } + if ((fnew >= 0 ? fabs (fm) : -fabs (fm)) != fm) { + xl = xm; + fl = fm; + xh = ans; + fh = fnew; + } else if ((fnew >= 0 ? fabs (fl) : -fabs (fl)) != fl) { + xh = ans; + fh = fnew; + } else if ((fnew >= 0 ? fabs (fh) : -fabs (fh)) != fh) { + xl = ans; + fl = fnew; + } else + return 0; + if (fabs (xh - xl) <= xacc) { + *root = ans; + return 1; + } + } + return 0; + } + + // normalised line shape R(x) with width w + double + dpw_lineshape (int broad, double x, double w) { + if (w <= 0) + return 0; + if (broad == 1) + return exp (-0.5 * x * x / (w * w)) / (sqrt (2 * PI) * w); + return w / (PI * (x * x + w * w)); + } + + int + DispersionPowder_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_DispersionPowder_physics_storage_struct->my_scattering; + return 1; + } + + int + DispersionPowder_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + struct DispersionPowder_physics_storage_struct* s = data_transfer.pointer_to_a_DispersionPowder_physics_storage_struct; + double ki, Ei, nsign = 1, C = VS2E * K2V * K2V; + int d, sign = 1; + + ki = sqrt (k_initial[0] * k_initial[0] + k_initial[1] * k_initial[1] + k_initial[2] * k_initial[2]); + Ei = C * ki * ki; + if (s->Temp > 0) { + nsign = 2; + if (rand01 () < 0.5) + sign = -1; + } + + if (s->meth == 0) { + // --- MCViNE E_Q_Kernel / (Lorentzian)Broadened_E_Q_Kernel --- + double Q = s->qmin + rand01 () * (s->qmax - s->qmin); + double E = sign * dpw_eval (s, s->E_expr, Q, 0); + if (s->broad) { + double w = dpw_eval (s, s->W_expr, Q, E); + if (w > 0) { + if (s->broad == 1) + E += w * randnorm (); + else + E += w * tan (PI * (rand01 () - 0.5)); + } + } + double Ef = Ei - E; + if (Ef <= 0) + return 0; + double kf = sqrt (Ef / C); + double cost = (ki * ki + kf * kf - Q * Q) / (2 * ki * kf); + if (cost * cost > 1) + return 0; + double sint = sqrt (1 - cost * cost), phi = 2 * PI * rand01 (); + double e1[3] = { k_initial[0] / ki, k_initial[1] / ki, k_initial[2] / ki }, e2[3], e3[3], n; + if (fabs (e1[0]) > 1e-10 || fabs (e1[1]) > 1e-10) { + e2[0] = -e1[1]; + e2[1] = e1[0]; + e2[2] = 0; + } else { + e2[0] = 1; + e2[1] = 0; + e2[2] = 0; + } + n = sqrt (e2[0] * e2[0] + e2[1] * e2[1] + e2[2] * e2[2]); + for (d = 0; d < 3; d++) + e2[d] /= n; + e3[0] = e1[1] * e2[2] - e1[2] * e2[1]; + e3[1] = e1[2] * e2[0] - e1[0] * e2[2]; + e3[2] = e1[0] * e2[1] - e1[1] * e2[0]; + for (d = 0; d < 3; d++) + k_final[d] = kf * (sint * cos (phi) * e2[d] + sint * sin (phi) * e3[d] + cost * e1[d]); + double S = dpw_eval (s, s->S_expr, Q, E); + if (s->Temp > 0) { + if (s->broad) { + // target [n(E)+1] (R(E-w) - R(E+w)), sampled from (R(E-w) + R(E+w))/2 + double w0 = dpw_eval (s, s->E_expr, Q, E), wd = dpw_eval (s, s->W_expr, Q, E); + double Rm = dpw_lineshape (s->broad, E - w0, wd), Rp = dpw_lineshape (s->broad, E + w0, wd); + if (Rm + Rp <= 0) + return 0; + S *= disp_bose (E, s->Temp) * fabs (Rm - Rp) / (Rm + Rp); + } else + S *= disp_bose (E, s->Temp); + } + *weight *= nsign * S * (kf / ki) * Q * (s->qmax - s->qmin) / (2 * ki * kf); + return 1; + } + + // --- method 1: focusing --- + double solid_angle; + Coords dir; + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + double Efmin = Ei - s->emax > 0 ? Ei - s->emax : 0, Efmax = Ei + s->emax, kf, Ef, pdf_factor; + if (s->broad == 0) { + // sharp excitation: root finding in kf + struct dpw_ctx c; + int i, nf = 0; + c.s = s; + c.Ei = Ei; + c.sign = sign; + for (d = 0; d < 3; d++) + c.ki[d] = k_initial[d]; + c.dir[0] = dir.x; + c.dir[1] = dir.y; + c.dir[2] = dir.z; + if (sign > 0) + Efmax = Ei; + else + Efmin = Ei; + double kmin = sqrt (Efmin / C), kmax = sqrt (Efmax / C); + double step = (kmax - kmin) / s->nsteps, fl, fh, root, xacc = 1e-9 * ki; + fl = dpw_f (kmin, &c); + for (i = 0; i < s->nsteps; i++) { + fh = dpw_f (kmin + (i + 1) * step, &c); + if (dpw_ridder (&c, kmin + i * step, kmin + (i + 1) * step, fl, fh, xacc, &root)) + if (nf == 0 || fabs (root - s->roots[nf - 1]) > 10 * xacc) + s->roots[nf++] = root; + fl = fh; + } + if (nf < 1) + return 0; + int index = (int)floor (rand01 () * nf); + if (index >= nf) + index = nf - 1; + kf = s->roots[index]; + if (kf <= 0) + return 0; + double dk = ki / 200; + double dfdk = fabs (-dpw_f (kf + 2 * dk, &c) + 8 * dpw_f (kf + dk, &c) - 8 * dpw_f (kf - dk, &c) + dpw_f (kf - 2 * dk, &c)) / (12 * dk); + if (!(dfdk > 0)) + return 0; + pdf_factor = nf * 2 * C * kf / dfdk; + Ef = C * kf * kf; + } else { + // broadened excitation: Ef uniform in the allowed window + Ef = Efmin + rand01 () * (Efmax - Efmin); + if (Ef <= 0) + return 0; + kf = sqrt (Ef / C); + pdf_factor = Efmax - Efmin; + } + double Q2 = 0, E = Ei - Ef; + k_final[0] = kf * dir.x; + k_final[1] = kf * dir.y; + k_final[2] = kf * dir.z; + for (d = 0; d < 3; d++) + Q2 += (k_initial[d] - k_final[d]) * (k_initial[d] - k_final[d]); + double Q = sqrt (Q2); + if (Q < s->qmin || Q > s->qmax) + return 0; + double S = dpw_eval (s, s->S_expr, Q, E); + if (s->broad) { + double w0 = dpw_eval (s, s->E_expr, Q, E), wd = dpw_eval (s, s->W_expr, Q, E); + if (s->Temp > 0) { + // detailed balance: [n(E)+1] (R(E-w) - R(E+w)), finite at E=0 + S *= disp_bose (E, s->Temp) * fabs (dpw_lineshape (s->broad, E - w0, wd) - dpw_lineshape (s->broad, E + w0, wd)); + } else + S *= dpw_lineshape (s->broad, E - w0, wd); + nsign = 1; // both signs are contained in the line shape + } else if (s->Temp > 0) + S *= disp_bose (E, s->Temp); + *weight *= solid_angle / (4 * PI) * nsign * pdf_factor * S * kf / ki; + return 1; + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_DISPERSIONPOWDER_DETECTOR + #define PROCESS_DISPERSIONPOWDER_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct DispersionPowder_physics_storage_struct DPW_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + char who[512]; + snprintf (who, 512, "DispersionPowder_process:%s", NAME_CURRENT_COMP); + memset (&DPW_storage, 0, sizeof (DPW_storage)); + + if (!strlen (E_Q)) { + fprintf (stderr, "%s: ERROR: E_Q must be given\n", who); + exit (-1); + } + if (sigma <= 0 || unit_cell_volume <= 0 || Emax <= 0 || root_nsteps < 1 || T < 0 || Qmin < 0 || Qmax <= Qmin) { + fprintf (stderr, "%s: ERROR: need sigma, unit_cell_volume, Emax, root_nsteps > 0, T >= 0 and 0 <= Qmin < Qmax\n", who); + exit (-1); + } + if (broadening < 0 || broadening > 2 || method < 0 || method > 1) { + fprintf (stderr, "%s: ERROR: broadening must be 0, 1 or 2 and method 0 or 1\n", who); + exit (-1); + } + DPW_storage.v.p[0] = p1; + DPW_storage.v.p[1] = p2; + DPW_storage.v.p[2] = p3; + DPW_storage.v.p[3] = p4; + DPW_storage.E_expr = disp_compile (E_Q, &DPW_storage.v, "E_Q", who); + DPW_storage.S_expr = disp_compile (strlen (S_Q) ? S_Q : "1", &DPW_storage.v, "S_Q", who); + DPW_storage.W_expr = disp_compile (strlen (W_Q) ? W_Q : "0", &DPW_storage.v, "W_Q", who); + DPW_storage.broad = broadening; + DPW_storage.meth = method; + DPW_storage.nsteps = root_nsteps; + DPW_storage.qmin = Qmin; + DPW_storage.qmax = Qmax; + DPW_storage.emax = Emax; + DPW_storage.Temp = T; + DPW_storage.my_scattering = packing_factor * sigma / unit_cell_volume * 100; // barn/AA^3 -> 1/m + DPW_storage.roots = malloc ((root_nsteps + 1) * sizeof (double)); + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Isotropic process + This_process.non_isotropic_rot_index = -1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = DispersionPowder; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_DispersionPowder_physics_storage_struct = &DPW_storage; + This_process.probability_for_scattering_function = &DispersionPowder_physics_my; + This_process.scattering_function = &DispersionPowder_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "%s: Error identifying Union_init component, %s is not a known component name.\n", who, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +FINALLY +%{ + te_free (DPW_storage.E_expr); + te_free (DPW_storage.S_expr); + te_free (DPW_storage.W_expr); + free (DPW_storage.roots); +%} + +END diff --git a/mcstas-comps/union/DispersionSingleXtal_process.comp b/mcstas-comps/union/DispersionSingleXtal_process.comp new file mode 100644 index 0000000000..b97251c757 --- /dev/null +++ b/mcstas-comps/union/DispersionSingleXtal_process.comp @@ -0,0 +1,426 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernel by Jiao Lin +* Date: 24.09.2026 +* Origin: DTU Physics +* +* Single crystal excitation with an analytical dispersion E(Q) and intensity S(Q) +* +* %D +* Port of the MCViNE kernel mccomponents/lib/kernels/sample/E_vQ_Kernel +* (J.Y.Y. Lin et al., Caltech/ORNL) to a Union process. +* +* Scattering from a single, sharp excitation branch in a single crystal, +* 1/sigma d2sigma/dOmega dEf = kf/ki S(Q) delta(E - E(Q)) / (4 pi) +* with E = Ei-Ef the energy transfer and Q = ki-kf the scattering vector in the local +* coordinate system of the process (rotate the process component to orient the +* crystal). The dispersion E(Q) [meV] and the intensity S(Q) [1] are given as +* expressions (tinyexpr syntax, see +* https://github.com/codeplea/tinyexpr) +* in the variables +* Qx, Qy, Qz [AA^-1] scattering vector in the local frame +* Q [AA^-1] |Q| +* h, k, l [rlu] Miller indices, h = a.Q/(2 pi) with the lattice vectors a, b, c +* E [meV] energy transfer (only in S_Q) +* p1..p4 free parameters +* e.g. E_Q="p1*sqrt(3-cos(2*pi*h)-cos(2*pi*k)-cos(2*pi*l))" or +* E_Q="10+5*sin(Qx+Qy+Qz)". Available functions include sin, cos, tan, asin, acos, +* atan, atan2, exp, ln, log10, sqrt, pow, abs, floor, ceil, gauss(x,mu,sigma) and +* hvs(x,x0,x1,..) (see Inhomogenous_incoherent_process). +* +* As in MCViNE, E(Q) is taken as an energy loss (0 < E < Emax). If T > 0 the +* excitation is also created with energy gain, E = -E(Q), and the intensities are +* multiplied by the Bose factors n(E)+1 and n(E), respectively. +* +* The final direction is sampled using the Union focusing settings of the geometry, +* and the final wave vector is found by root finding of E(Q(kf)) - (Ei-Ef) in +* root_nsteps intervals. Rays without a solution are absorbed (weight 0). +* The inverse penetration depth is sigma/unit_cell_volume. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* %P +* INPUT PARAMETERS: +* E_Q: [string] Expression for the dispersion E(Q) [meV] +* S_Q: [string] Expression for the intensity S(Q,E) [1] +* sigma: [barns] Scattering cross section of the unit cell, defines the inverse penetration depth +* unit_cell_volume: [AA^3] Unit cell volume +* Emax: [meV] Maximum energy transfer considered (|E| < Emax) +* T: [K] Temperature. 0: energy loss only and no Bose factor. >0: loss and gain with Bose factors +* ax: [AA] First lattice vector, x coordinate (used for h) +* ay: [AA] First lattice vector, y coordinate +* az: [AA] First lattice vector, z coordinate +* bx: [AA] Second lattice vector, x coordinate (used for k) +* by: [AA] Second lattice vector, y coordinate +* bz: [AA] Second lattice vector, z coordinate +* cx: [AA] Third lattice vector, x coordinate (used for l) +* cy: [AA] Third lattice vector, y coordinate +* cz: [AA] Third lattice vector, z coordinate +* p1: [1] Free parameter available in the expressions +* p2: [1] Free parameter available in the expressions +* p3: [1] Free parameter available in the expressions +* p4: [1] Free parameter available in the expressions +* root_nsteps: [1] Number of intervals in kf in which roots are searched for +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* DSX_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* +* %E +******************************************************************************/ + +DEFINE COMPONENT DispersionSingleXtal_process + +SETTING PARAMETERS(string E_Q="", string S_Q="1", double sigma=1, double unit_cell_volume=1, double Emax=100, double T=0, + double ax=6.283185307179586, double ay=0, double az=0, double bx=0, double by=6.283185307179586, double bz=0, + double cx=0, double cy=0, double cz=6.283185307179586, double p1=0, double p2=0, double p3=0, double p4=0, + int root_nsteps=200, double packing_factor=1, double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this DispersionSingleXtal_process component" + #endif + + %include "tinyexpr.h" + %include "tinyexpr.c" + + #ifndef DISPERSION_EXPR_TOOLS + #define DISPERSION_EXPR_TOOLS + #define DISP_T2E (1.0 / 11.605) // Kelvin to meV + + // variables available in the expressions + struct disp_expr_vars { + double Qx, Qy, Qz, Q, h, k, l, E, p[4]; + }; + + // compile an expression with the standard variable set + te_expr* + disp_compile (const char* expr, struct disp_expr_vars* v, const char* what, const char* comp) { + int err = 0; + te_variable vars[] = { { "Qx", &v->Qx }, { "Qy", &v->Qy }, { "Qz", &v->Qz }, { "Q", &v->Q }, { "h", &v->h }, + { "k", &v->k }, { "l", &v->l }, { "E", &v->E }, { "p1", &v->p[0] }, { "p2", &v->p[1] }, + { "p3", &v->p[2] }, { "p4", &v->p[3] } }; + te_expr* e = te_compile (expr, vars, 12, &err); + if (!e) { + fprintf (stderr, "%s: ERROR: could not parse %s expression \"%s\" near character %d\n", comp, what, expr, err); + exit (-1); + } + return e; + } + + // Bose factor: n+1 for energy loss (E>0), n for energy gain + double + disp_bose (double E, double T) { + double n = 1.0 / (exp (fabs (E) / (T * DISP_T2E)) - 1.0); + return E > 0 ? 1.0 + n : n; + } + #endif + + // Very important to add a pointer to this struct in the union-lib.c file + struct DispersionSingleXtal_physics_storage_struct { + struct disp_expr_vars v; + te_expr* E_expr; + te_expr* S_expr; + double a[3][3]; // lattice vectors (rows), h = a[0].Q/(2 pi) + double my_scattering; + double emax, Temp; + int nsteps; + double* roots; + }; + + // context for f(kf) = sign*E(Q(kf)) - (Ei - Ef) + struct dsx_ctx { + struct DispersionSingleXtal_physics_storage_struct* s; + double ki[3], Ei, dir[3]; + int sign; + }; + + double + dsx_E_of_Q (struct DispersionSingleXtal_physics_storage_struct* s, double* Q) { + struct disp_expr_vars* v = &s->v; + v->Qx = Q[0]; + v->Qy = Q[1]; + v->Qz = Q[2]; + v->Q = sqrt (Q[0] * Q[0] + Q[1] * Q[1] + Q[2] * Q[2]); + v->h = (s->a[0][0] * Q[0] + s->a[0][1] * Q[1] + s->a[0][2] * Q[2]) / (2 * PI); + v->k = (s->a[1][0] * Q[0] + s->a[1][1] * Q[1] + s->a[1][2] * Q[2]) / (2 * PI); + v->l = (s->a[2][0] * Q[0] + s->a[2][1] * Q[1] + s->a[2][2] * Q[2]) / (2 * PI); + return te_eval (s->E_expr); + } + + double + dsx_f (double kf, struct dsx_ctx* c) { + double Q[3]; + int d; + for (d = 0; d < 3; d++) + Q[d] = c->ki[d] - kf * c->dir[d]; + return c->sign * dsx_E_of_Q (c->s, Q) - (c->Ei - VS2E * K2V * K2V * kf * kf); + } + + // Ridder's method with bracketing values fl, fh. Returns 1 on success + int + dsx_ridder (struct dsx_ctx* c, double xl, double xh, double fl, double fh, double xacc, double* root) { + double ans = xl, xm, fm, s, xnew, fnew; + int j; + if (fl == 0) { + *root = xl; + return 1; + } + if (fh == 0) { + *root = xh; + return 1; + } + if (fl * fh > 0) + return 0; + for (j = 0; j < 60; j++) { + xm = 0.5 * (xl + xh); + fm = dsx_f (xm, c); + s = sqrt (fm * fm - fl * fh); + if (s == 0.0) { + *root = xm; + return 1; + } + xnew = xm + (xm - xl) * ((fl >= fh ? 1.0 : -1.0) * fm / s); + if (j && fabs (xnew - ans) <= xacc) { + *root = xnew; + return 1; + } + ans = xnew; + fnew = dsx_f (ans, c); + if (fnew == 0.0) { + *root = ans; + return 1; + } + if ((fnew >= 0 ? fabs (fm) : -fabs (fm)) != fm) { + xl = xm; + fl = fm; + xh = ans; + fh = fnew; + } else if ((fnew >= 0 ? fabs (fl) : -fabs (fl)) != fl) { + xh = ans; + fh = fnew; + } else if ((fnew >= 0 ? fabs (fh) : -fabs (fh)) != fh) { + xl = ans; + fl = fnew; + } else + return 0; + if (fabs (xh - xl) <= xacc) { + *root = ans; + return 1; + } + } + return 0; + } + + int + DispersionSingleXtal_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_DispersionSingleXtal_physics_storage_struct->my_scattering; + return 1; + } + + // E_vQ_Kernel::S + int + DispersionSingleXtal_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + struct DispersionSingleXtal_physics_storage_struct* s = data_transfer.pointer_to_a_DispersionSingleXtal_physics_storage_struct; + struct dsx_ctx c; + double ki, solid_angle, kmin, kmax, nsign = 1, Efmin, Efmax; + Coords dir; + int d, i, nf = 0; + + c.s = s; + for (d = 0; d < 3; d++) + c.ki[d] = k_initial[d]; + ki = sqrt (k_initial[0] * k_initial[0] + k_initial[1] * k_initial[1] + k_initial[2] * k_initial[2]); + c.Ei = VS2E * K2V * K2V * ki * ki; + + // energy loss (as MCViNE), or with T>0 also energy gain + c.sign = 1; + if (s->Temp > 0) { + nsign = 2; + if (rand01 () < 0.5) + c.sign = -1; + } + if (c.sign > 0) { + Efmin = c.Ei - s->emax > 0 ? c.Ei - s->emax : 0; + Efmax = c.Ei; + } else { + Efmin = c.Ei; + Efmax = c.Ei + s->emax; + } + kmin = sqrt (Efmin / VS2E) * V2K; + kmax = sqrt (Efmax / VS2E) * V2K; + + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + c.dir[0] = dir.x; + c.dir[1] = dir.y; + c.dir[2] = dir.z; + + // roots of f(kf) in root_nsteps intervals + double step = (kmax - kmin) / s->nsteps, fl, fh, root, xacc = 1e-9 * ki; + fl = dsx_f (kmin, &c); + for (i = 0; i < s->nsteps; i++) { + double xl = kmin + i * step, xh = kmin + (i + 1) * step; + fh = dsx_f (xh, &c); + if (dsx_ridder (&c, xl, xh, fl, fh, xacc, &root)) { + if (nf == 0 || fabs (root - s->roots[nf - 1]) > 10 * xacc) + s->roots[nf++] = root; + } + fl = fh; + } + if (nf < 1) + return 0; + int index = (int)floor (rand01 () * nf); + if (index >= nf) + index = nf - 1; + double kf = s->roots[index]; + if (kf <= 0) + return 0; + + // Jacobian |df/dkf|, 5 point stencil as in MCViNE + double dk = ki / 200; + double dfdk = fabs (-dsx_f (kf + 2 * dk, &c) + 8 * dsx_f (kf + dk, &c) - 8 * dsx_f (kf - dk, &c) + dsx_f (kf - 2 * dk, &c)) / (12 * dk); + if (!(dfdk > 0)) + return 0; + + double Q[3], Ef = VS2E * K2V * K2V * kf * kf, E = c.Ei - Ef; + for (d = 0; d < 3; d++) { + Q[d] = c.ki[d] - kf * c.dir[d]; + k_final[d] = kf * c.dir[d]; + } + if (fabs (E - c.sign * dsx_E_of_Q (s, Q)) > s->emax * 1e-5) + return 0; + s->v.E = E; // Q variables are set by dsx_E_of_Q above + double S = te_eval (s->S_expr); + if (s->Temp > 0) + S *= disp_bose (E, s->Temp); + + *weight *= solid_angle / (4 * PI) * nsign * nf * (2 * VS2E * K2V * K2V * kf / dfdk) * S * kf / ki; + return 1; + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_DISPERSIONSINGLEXTAL_DETECTOR + #define PROCESS_DISPERSIONSINGLEXTAL_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct DispersionSingleXtal_physics_storage_struct DSX_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + char who[512]; + snprintf (who, 512, "DispersionSingleXtal_process:%s", NAME_CURRENT_COMP); + memset (&DSX_storage, 0, sizeof (DSX_storage)); + + if (!strlen (E_Q)) { + fprintf (stderr, "%s: ERROR: E_Q must be given\n", who); + exit (-1); + } + if (sigma <= 0 || unit_cell_volume <= 0 || Emax <= 0 || root_nsteps < 1 || T < 0) { + fprintf (stderr, "%s: ERROR: sigma, unit_cell_volume, Emax and root_nsteps must be positive, T >= 0\n", who); + exit (-1); + } + DSX_storage.v.p[0] = p1; + DSX_storage.v.p[1] = p2; + DSX_storage.v.p[2] = p3; + DSX_storage.v.p[3] = p4; + DSX_storage.E_expr = disp_compile (E_Q, &DSX_storage.v, "E_Q", who); + DSX_storage.S_expr = disp_compile (strlen (S_Q) ? S_Q : "1", &DSX_storage.v, "S_Q", who); + DSX_storage.a[0][0] = ax; + DSX_storage.a[0][1] = ay; + DSX_storage.a[0][2] = az; + DSX_storage.a[1][0] = bx; + DSX_storage.a[1][1] = by; + DSX_storage.a[1][2] = bz; + DSX_storage.a[2][0] = cx; + DSX_storage.a[2][1] = cy; + DSX_storage.a[2][2] = cz; + DSX_storage.my_scattering = packing_factor * sigma / unit_cell_volume * 100; // barn/AA^3 -> 1/m + DSX_storage.emax = Emax; + DSX_storage.Temp = T; + DSX_storage.nsteps = root_nsteps; + DSX_storage.roots = malloc ((root_nsteps + 1) * sizeof (double)); + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Single crystal: the process has its own orientation + This_process.non_isotropic_rot_index = 1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = DispersionSingleXtal; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_DispersionSingleXtal_physics_storage_struct = &DSX_storage; + This_process.probability_for_scattering_function = &DispersionSingleXtal_physics_my; + This_process.scattering_function = &DispersionSingleXtal_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "%s: Error identifying Union_init component, %s is not a known component name.\n", who, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +FINALLY +%{ + te_free (DSX_storage.E_expr); + te_free (DSX_storage.S_expr); + free (DSX_storage.roots); +%} + +END diff --git a/mcstas-comps/union/ElasticSQ_process.comp b/mcstas-comps/union/ElasticSQ_process.comp new file mode 100644 index 0000000000..def2500915 --- /dev/null +++ b/mcstas-comps/union/ElasticSQ_process.comp @@ -0,0 +1,459 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernels by Jiao Lin +* Date: 24.09.2026 +* Origin: DTU Physics +* +* Elastic scattering from a structure factor S(|Q|) (isotropic) or S(Q) (single crystal diffuse) +* +* %D +* Port of the MCViNE kernels mccomponents/lib/kernels/sample/SQkernel (with GridSQ and +* SQ_fromexpression) and SvQkernel (with GridSvQ and SvQ_fromexpression) +* (J.Y.Y. Lin et al., Caltech/ORNL) to a Union process. +* +* Elastic scattering with +* 1/sigma dsigma/dOmega = S(Q) / (4 pi) +* where Q = ki-kf. With crystal=0, S depends on |Q| only (liquids, glasses, powders, +* e.g. a measured structure factor). With crystal=1, S depends on the Q vector in the +* local coordinate system of the process, which can be rotated to orient the crystal +* (e.g. diffuse scattering of a single crystal). +* +* S is given either as an expression (tinyexpr syntax) in the variables +* Q [AA^-1], Qx, Qy, Qz [AA^-1], h, k, l [rlu, h = a.Q/(2 pi)] and p1..p4 +* or as a table in SQ_file: +* crystal=0: two columns "Q S" (Q ascending), linearly interpolated, +* crystal=1: four columns "x y z S" on a regular 3D grid (any order), trilinearly +* interpolated, with x y z = Qx Qy Qz [AA^-1] or, with grid_hkl=1, h k l. +* S is zero outside the tabulated range. +* +* The final direction is sampled with the Union focusing of the geometry. For crystal=0, +* method=0 samples Q uniformly in [Qmin, min(Qmax, 2ki)] instead (as MCViNE, no focusing), +* which is more efficient for a structure factor with sharp features. +* The inverse penetration depth is sigma/unit_cell_volume. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* %P +* INPUT PARAMETERS: +* SQ_file: [string] Table of S: "Q S" (crystal=0) or "x y z S" (crystal=1) +* SQ_expr: [string] Expression for S, used when no SQ_file is given +* crystal: [1] 0: S(|Q|), isotropic; 1: S(Q vector), single crystal +* grid_hkl: [1] crystal=1: 1 if the grid in SQ_file is given in h k l instead of Qx Qy Qz +* sigma: [barns] Scattering cross section of the unit cell, defines the inverse penetration depth +* unit_cell_volume: [AA^3] Unit cell volume +* Qmin: [AA^-1] Minimum |Q| +* Qmax: [AA^-1] Maximum |Q| +* method: [1] crystal=0: 0 MCViNE sampling in Q (no focusing), 1 focusing +* ax: [AA] First lattice vector, x coordinate (used for h) +* ay: [AA] First lattice vector, y coordinate +* az: [AA] First lattice vector, z coordinate +* bx: [AA] Second lattice vector, x coordinate (used for k) +* by: [AA] Second lattice vector, y coordinate +* bz: [AA] Second lattice vector, z coordinate +* cx: [AA] Third lattice vector, x coordinate (used for l) +* cy: [AA] Third lattice vector, y coordinate +* cz: [AA] Third lattice vector, z coordinate +* p1: [1] Free parameter available in SQ_expr +* p2: [1] Free parameter available in SQ_expr +* p3: [1] Free parameter available in SQ_expr +* p4: [1] Free parameter available in SQ_expr +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* ESQ_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* +* %E +******************************************************************************/ + +DEFINE COMPONENT ElasticSQ_process + +SETTING PARAMETERS(string SQ_file="", string SQ_expr="", int crystal=0, int grid_hkl=0, double sigma=1, double unit_cell_volume=1, + double Qmin=0, double Qmax=100, int method=1, + double ax=6.283185307179586, double ay=0, double az=0, double bx=0, double by=6.283185307179586, double bz=0, + double cx=0, double cy=0, double cz=6.283185307179586, double p1=0, double p2=0, double p3=0, double p4=0, + double packing_factor=1, double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this ElasticSQ_process component" + #endif + + %include "read_table-lib" + %include "tinyexpr.h" + %include "tinyexpr.c" + + #ifndef ELASTICSQ_SHARE + #define ELASTICSQ_SHARE + struct esq_vars { + double Qx, Qy, Qz, Q, h, k, l, p[4]; + }; + #endif + + // Very important to add a pointer to this struct in the union-lib.c file + struct ElasticSQ_physics_storage_struct { + struct esq_vars v; + te_expr* expr; + int use_expr, cryst, hkl_grid, meth; + double a[3][3]; + double qmin, qmax, my_scattering; + // 1D table + int n1; + double *tq, *tS; + // 3D grid + int n[3]; + double *gax[3], *gS; // axis values and S[(i*n1+j)*n2+k] + }; + + // index i with a[i] <= x <= a[i+1], -1 if outside + int + esq_locate (double* a, int n, double x) { + int lo = 0, hi = n - 1; + if (n < 2 || x < a[0] || x > a[n - 1]) + return -1; + while (hi - lo > 1) { + int mid = (lo + hi) / 2; + if (a[mid] <= x) + lo = mid; + else + hi = mid; + } + return lo; + } + + double + esq_S (struct ElasticSQ_physics_storage_struct* s, double* Q) { + struct esq_vars* v = &s->v; + v->Qx = Q[0]; + v->Qy = Q[1]; + v->Qz = Q[2]; + v->Q = sqrt (Q[0] * Q[0] + Q[1] * Q[1] + Q[2] * Q[2]); + if (v->Q < s->qmin || v->Q > s->qmax) + return 0; + v->h = (s->a[0][0] * Q[0] + s->a[0][1] * Q[1] + s->a[0][2] * Q[2]) / (2 * PI); + v->k = (s->a[1][0] * Q[0] + s->a[1][1] * Q[1] + s->a[1][2] * Q[2]) / (2 * PI); + v->l = (s->a[2][0] * Q[0] + s->a[2][1] * Q[1] + s->a[2][2] * Q[2]) / (2 * PI); + if (s->use_expr) + return te_eval (s->expr); + if (!s->cryst) { + int i = esq_locate (s->tq, s->n1, v->Q); + if (i < 0) + return 0; + double t = (v->Q - s->tq[i]) / (s->tq[i + 1] - s->tq[i]); + return (1 - t) * s->tS[i] + t * s->tS[i + 1]; + } + double x[3] = { v->Qx, v->Qy, v->Qz }, r[3]; + int i0[3], d, dx, dy, dz; + if (s->hkl_grid) { + x[0] = v->h; + x[1] = v->k; + x[2] = v->l; + } + for (d = 0; d < 3; d++) { + i0[d] = esq_locate (s->gax[d], s->n[d], x[d]); + if (i0[d] < 0) + return 0; + r[d] = (x[d] - s->gax[d][i0[d]]) / (s->gax[d][i0[d] + 1] - s->gax[d][i0[d]]); + } + double res = 0; + for (dx = 0; dx < 2; dx++) + for (dy = 0; dy < 2; dy++) + for (dz = 0; dz < 2; dz++) { + double w = (dx ? r[0] : 1 - r[0]) * (dy ? r[1] : 1 - r[1]) * (dz ? r[2] : 1 - r[2]); + if (w > 0) + res += w * s->gS[((long)(i0[0] + dx) * s->n[1] + (i0[1] + dy)) * s->n[2] + (i0[2] + dz)]; + } + return res; + } + + int + esq_cmp (const void* a, const void* b) { + double x = *(const double*)a, y = *(const double*)b; + return x < y ? -1 : x > y; + } + + // unique sorted values of column c (stride 4) of data + int + esq_axis (double* data, long nrows, int c, double** out) { + double* v = malloc (nrows * sizeof (double)); + long i; + int n = 0; + for (i = 0; i < nrows; i++) + v[i] = data[4 * i + c]; + qsort (v, nrows, sizeof (double), esq_cmp); + for (i = 0; i < nrows; i++) + if (n == 0 || fabs (v[i] - v[n - 1]) > 1e-9 * (fabs (v[i]) + 1)) + v[n++] = v[i]; + *out = v; + return n; + } + + void + esq_read_file (struct ElasticSQ_physics_storage_struct* s, const char* file, const char* comp) { + FILE* fp = Open_File ((char*)file, "r", NULL); + char line[4096]; + int ncol = s->cryst ? 4 : 2; + long nrows = 0, cap = 0, i; + double* data = NULL; + if (!fp) { + fprintf (stderr, "%s: ERROR: cannot open SQ_file %s\n", comp, file); + exit (-1); + } + while (fgets (line, sizeof (line), fp)) { + double x[4]; + char* c = line; + while (*c == ' ' || *c == '\t') + c++; + if (*c == '#' || *c == '%') + continue; + if (sscanf (c, "%lf %lf %lf %lf", &x[0], &x[1], &x[2], &x[3]) < ncol) + continue; + if (nrows == cap) { + cap = cap ? 2 * cap : 1024; + data = realloc (data, cap * 4 * sizeof (double)); + } + for (i = 0; i < 4; i++) + data[4 * nrows + i] = x[i]; + nrows++; + } + fclose (fp); + if (!s->cryst) { + s->n1 = nrows; + s->tq = malloc (nrows * sizeof (double)); + s->tS = malloc (nrows * sizeof (double)); + for (i = 0; i < nrows; i++) { + s->tq[i] = data[4 * i]; + s->tS[i] = data[4 * i + 1]; + if (i && s->tq[i] <= s->tq[i - 1]) { + fprintf (stderr, "%s: ERROR: Q values in %s must be ascending\n", comp, file); + exit (-1); + } + } + if (nrows < 2) { + fprintf (stderr, "%s: ERROR: %s must contain at least 2 rows\n", comp, file); + exit (-1); + } + } else { + int d; + for (d = 0; d < 3; d++) + s->n[d] = esq_axis (data, nrows, d, &s->gax[d]); + if ((long)s->n[0] * s->n[1] * s->n[2] != nrows || s->n[0] < 2 || s->n[1] < 2 || s->n[2] < 2) { + fprintf (stderr, "%s: ERROR: %s is not a complete regular grid (%ld rows, %d x %d x %d axis values)\n", comp, file, nrows, s->n[0], s->n[1], + s->n[2]); + exit (-1); + } + s->gS = calloc (nrows, sizeof (double)); + for (i = 0; i < nrows; i++) { + int idx[3]; + for (d = 0; d < 3; d++) { + // index of the axis value closest to this grid coordinate + double x = data[4 * i + d]; + int j = esq_locate (s->gax[d], s->n[d], x); + if (j < 0) + j = 0; + idx[d] = (j + 1 < s->n[d] && fabs (s->gax[d][j + 1] - x) < fabs (s->gax[d][j] - x)) ? j + 1 : j; + } + s->gS[((long)idx[0] * s->n[1] + idx[1]) * s->n[2] + idx[2]] = data[4 * i + 3]; + } + } + free (data); + } + + int + ElasticSQ_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_ElasticSQ_physics_storage_struct->my_scattering; + return 1; + } + + int + ElasticSQ_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + struct ElasticSQ_physics_storage_struct* s = data_transfer.pointer_to_a_ElasticSQ_physics_storage_struct; + double ki = sqrt (k_initial[0] * k_initial[0] + k_initial[1] * k_initial[1] + k_initial[2] * k_initial[2]), Q[3]; + int d; + + if (!s->cryst && s->meth == 0) { + // --- MCViNE SQkernel: Q uniform in [Qmin, min(Qmax, 2ki)], random azimuth --- + double Q1 = s->qmin, Q2 = s->qmax < 2 * ki ? s->qmax : 2 * ki; + if (Q2 <= Q1) + return 0; + double q = Q1 + rand01 () * (Q2 - Q1); + double cost = 1 - q * q / (2 * ki * ki), sint = sqrt (1 - cost * cost), phi = 2 * PI * rand01 (); + double e1[3] = { k_initial[0] / ki, k_initial[1] / ki, k_initial[2] / ki }, e2[3], e3[3], n; + if (fabs (e1[0]) > 1e-10 || fabs (e1[1]) > 1e-10) { + e2[0] = -e1[1]; + e2[1] = e1[0]; + e2[2] = 0; + } else { + e2[0] = 1; + e2[1] = 0; + e2[2] = 0; + } + n = sqrt (e2[0] * e2[0] + e2[1] * e2[1] + e2[2] * e2[2]); + for (d = 0; d < 3; d++) + e2[d] /= n; + e3[0] = e1[1] * e2[2] - e1[2] * e2[1]; + e3[1] = e1[2] * e2[0] - e1[0] * e2[2]; + e3[2] = e1[0] * e2[1] - e1[1] * e2[0]; + for (d = 0; d < 3; d++) { + k_final[d] = ki * (sint * cos (phi) * e2[d] + sint * sin (phi) * e3[d] + cost * e1[d]); + Q[d] = k_initial[d] - k_final[d]; + } + *weight *= esq_S (s, Q) * q * (Q2 - Q1) / (2 * ki * ki); + return 1; + } + + // --- focusing (SvQkernel samples 4pi uniformly) --- + double solid_angle; + Coords dir; + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + k_final[0] = ki * dir.x; + k_final[1] = ki * dir.y; + k_final[2] = ki * dir.z; + for (d = 0; d < 3; d++) + Q[d] = k_initial[d] - k_final[d]; + *weight *= solid_angle / (4 * PI) * esq_S (s, Q); + return 1; + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_ELASTICSQ_DETECTOR + #define PROCESS_ELASTICSQ_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct ElasticSQ_physics_storage_struct ESQ_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + char who[512]; + snprintf (who, 512, "ElasticSQ_process:%s", NAME_CURRENT_COMP); + memset (&ESQ_storage, 0, sizeof (ESQ_storage)); + + if (sigma <= 0 || unit_cell_volume <= 0 || Qmin < 0 || Qmax <= Qmin) { + fprintf (stderr, "%s: ERROR: need sigma > 0, unit_cell_volume > 0 and 0 <= Qmin < Qmax\n", who); + exit (-1); + } + ESQ_storage.cryst = crystal ? 1 : 0; + ESQ_storage.hkl_grid = grid_hkl ? 1 : 0; + ESQ_storage.meth = method; + ESQ_storage.qmin = Qmin; + ESQ_storage.qmax = Qmax; + ESQ_storage.a[0][0] = ax; + ESQ_storage.a[0][1] = ay; + ESQ_storage.a[0][2] = az; + ESQ_storage.a[1][0] = bx; + ESQ_storage.a[1][1] = by; + ESQ_storage.a[1][2] = bz; + ESQ_storage.a[2][0] = cx; + ESQ_storage.a[2][1] = cy; + ESQ_storage.a[2][2] = cz; + ESQ_storage.v.p[0] = p1; + ESQ_storage.v.p[1] = p2; + ESQ_storage.v.p[2] = p3; + ESQ_storage.v.p[3] = p4; + if (strlen (SQ_file) && strcmp (SQ_file, "NULL") && strcmp (SQ_file, "0")) + esq_read_file (&ESQ_storage, SQ_file, who); + else if (strlen (SQ_expr)) { + int err = 0; + struct esq_vars* v = &ESQ_storage.v; + te_variable vars[] = { { "Qx", &v->Qx }, { "Qy", &v->Qy }, { "Qz", &v->Qz }, { "Q", &v->Q }, { "h", &v->h }, { "k", &v->k }, + { "l", &v->l }, { "p1", &v->p[0] }, { "p2", &v->p[1] }, { "p3", &v->p[2] }, { "p4", &v->p[3] } }; + ESQ_storage.expr = te_compile (SQ_expr, vars, 11, &err); + if (!ESQ_storage.expr) { + fprintf (stderr, "%s: ERROR: could not parse SQ_expr \"%s\" near character %d\n", who, SQ_expr, err); + exit (-1); + } + ESQ_storage.use_expr = 1; + } else { + fprintf (stderr, "%s: ERROR: SQ_file or SQ_expr must be given\n", who); + exit (-1); + } + ESQ_storage.my_scattering = packing_factor * sigma / unit_cell_volume * 100; // barn/AA^3 -> 1/m + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Single crystal has its own orientation, otherwise isotropic + This_process.non_isotropic_rot_index = crystal ? 1 : -1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = ElasticSQ; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_ElasticSQ_physics_storage_struct = &ESQ_storage; + This_process.probability_for_scattering_function = &ElasticSQ_physics_my; + This_process.scattering_function = &ElasticSQ_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "%s: Error identifying Union_init component, %s is not a known component name.\n", who, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +FINALLY +%{ + int d; + if (ESQ_storage.use_expr) + te_free (ESQ_storage.expr); + free (ESQ_storage.tq); + free (ESQ_storage.tS); + free (ESQ_storage.gS); + for (d = 0; d < 3; d++) + free (ESQ_storage.gax[d]); +%} + +END diff --git a/mcstas-comps/union/IncoherentElastic_process.comp b/mcstas-comps/union/IncoherentElastic_process.comp new file mode 100644 index 0000000000..3e86031ec3 --- /dev/null +++ b/mcstas-comps/union/IncoherentElastic_process.comp @@ -0,0 +1,183 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernel by Jiao Lin +* Date: 24.09.2026 +* Origin: DTU Physics +* +* Incoherent elastic scattering with a Debye-Waller factor +* +* %D +* Port of the MCViNE kernel mccomponents/lib/kernels/sample/phonon/IncoherentElastic +* (J.Y.Y. Lin et al., Caltech/ORNL) to a Union process. +* +* Isotropic elastic incoherent scattering reduced by the Debye-Waller factor, +* 1/sigma_inc dsigma/dOmega = exp(-DW_core*Q^2)/(4 pi) +* with DW_core = /3 (2W = DW_core*Q^2). The inverse penetration depth is +* sigma_inc/V, i.e. the Debye-Waller reduction is applied as a weight factor as in +* MCViNE, so the attenuation corresponds to the full incoherent cross section. +* With DW_core=0 the process is identical to Incoherent_process. +* +* The final direction is sampled using the Union focusing settings of the geometry. +* (MCViNE samples the polar angle uniformly and weights with sin(theta), which is +* equivalent without focusing.) +* +* DW_core can be obtained from a phonon DOS, e.g. with the MCViNE script +* mcvine-debye-waller-core-from-phonon-dos, or as printed by +* IncoherentOnePhonon_process for the same DOS, mass and temperature. +* +* Note on combining processes: as in MCViNE, the inverse penetration depth of this +* process is the full cross section sigma/V_uc, and the kernel physics is carried by +* the weight. When several processes describing parts of the same cross section are +* combined in one material, e.g. +* IncoherentElastic_process + IncoherentOnePhonon_process, +* their attenuations add up, so the attenuation is overestimated for thick samples. +* Single scattering intensities are not affected. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* Example: IncoherentElastic_process(sigma_inc=5.08, unit_cell_volume=13.827, DW_core=0.0067) +* +* %P +* INPUT PARAMETERS: +* sigma_inc: [barns] Incoherent scattering cross section of the unit cell +* unit_cell_volume: [AA^3] Unit cell volume +* DW_core: [AA^2] Debye-Waller core, 2W = DW_core*Q^2 (= /3) +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* IEL_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* +* %E +******************************************************************************/ + +DEFINE COMPONENT IncoherentElastic_process + +SETTING PARAMETERS(double sigma_inc=5.08, double unit_cell_volume=13.827, double DW_core=0, double packing_factor=1, + double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this IncoherentElastic_process component" + #endif + + // Very important to add a pointer to this struct in the union-lib.c file + struct IncoherentElastic_physics_storage_struct { + double my_scattering; + double dw_core; + }; + + int + IncoherentElastic_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_IncoherentElastic_physics_storage_struct->my_scattering; + return 1; + } + + int + IncoherentElastic_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + double k_length = sqrt (k_initial[0] * k_initial[0] + k_initial[1] * k_initial[1] + k_initial[2] * k_initial[2]); + double solid_angle, Q2; + Coords k_out; + + focus_data->focusing_function (&k_out, &solid_angle, focus_data); + NORM (k_out.x, k_out.y, k_out.z); + k_final[0] = k_out.x * k_length; + k_final[1] = k_out.y * k_length; + k_final[2] = k_out.z * k_length; + Q2 = (k_initial[0] - k_final[0]) * (k_initial[0] - k_final[0]) + (k_initial[1] - k_final[1]) * (k_initial[1] - k_final[1]) + + (k_initial[2] - k_final[2]) * (k_initial[2] - k_final[2]); + *weight *= solid_angle * 0.25 / PI * exp (-data_transfer.pointer_to_a_IncoherentElastic_physics_storage_struct->dw_core * Q2); + return 1; + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_INCOHERENTELASTIC_DETECTOR + #define PROCESS_INCOHERENTELASTIC_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct IncoherentElastic_physics_storage_struct IEL_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + if (sigma_inc < 0 || unit_cell_volume <= 0 || DW_core < 0) { + fprintf (stderr, "IncoherentElastic_process:%s: ERROR: sigma_inc and DW_core must be >= 0 and unit_cell_volume > 0\n", NAME_CURRENT_COMP); + exit (-1); + } + IEL_storage.my_scattering = packing_factor * sigma_inc / unit_cell_volume * 100; // barn/AA^3 -> 1/m + IEL_storage.dw_core = DW_core; + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Isotropic process + This_process.non_isotropic_rot_index = -1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = IncoherentElastic; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_IncoherentElastic_physics_storage_struct = &IEL_storage; + This_process.probability_for_scattering_function = &IncoherentElastic_physics_my; + This_process.scattering_function = &IncoherentElastic_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "IncoherentElastic_process:%s: Error identifying Union_init component, %s is not a known component name.\n", NAME_CURRENT_COMP, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +END diff --git a/mcstas-comps/union/IncoherentOnePhonon_process.comp b/mcstas-comps/union/IncoherentOnePhonon_process.comp new file mode 100644 index 0000000000..31c36c04cc --- /dev/null +++ b/mcstas-comps/union/IncoherentOnePhonon_process.comp @@ -0,0 +1,589 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernel by Jiao Lin +* Date: 23.09.2026 +* Origin: DTU Physics +* +* Incoherent inelastic one-phonon scattering from a phonon density of states (DOS) +* +* %D +* Port of the MCViNE kernels +* mccomponents/lib/kernels/sample/phonon/IncoherentInelastic and +* IncoherentInelastic_EnergyFocusing (J.Y.Y. Lin et al., Caltech/ORNL) to a Union process. +* +* One-phonon incoherent scattering in the incoherent (and cubic) approximation, +* d2sigma/dOmega dEf = sigma_inc/(4 pi) kf/ki exp(-2W) hbar^2Q^2/(2M) Z(|w|)/|w| [n(w)+1 or n(w)] +* with Z the normalised phonon DOS and n the Bose factor. The Debye-Waller factor +* 2W = DW_core*Q^2 is by default calculated from the DOS. +* +* The DOS file may be either +* - a 2-column ASCII file "E Z" (additional columns ignored). The energy unit is +* taken from DOS_unit, or else from a comment line whose first word contains +* "meV" or "TeraHz"/"THz" (e.g. "# E(meV) Z" or "# Frequency(TeraHz) DOS", as in +* MCViNE). THz means frequency, not angular frequency. Default is meV. +* - a binary MCViNE IDF "DOS" file (frequency in THz) +* The DOS is prepared as in MCViNE (utils.nice_dos): curves with fewer than 500 points +* are resampled on a 500 point grid starting at E=0, the low energy part is replaced +* by a fitted parabola (smoothing the curve first if the fit is bad), and the curve +* is normalised to unit area. +* +* The final energy is sampled uniformly in [Ei-Emax, Ei+Emax] (limited to Ef>0), with +* Emax the maximum energy of the DOS. The final direction is sampled using the Union +* focusing settings of the geometry (default 4pi). With dEf_focus > 0 the final +* energy is instead only sampled within [Ef_focus-dEf_focus/2, Ef_focus+dEf_focus/2] +* (energy focusing, e.g. for indirect geometry spectrometers). +* +* Note on combining processes: as in MCViNE, the inverse penetration depth of this +* process is the full cross section sigma/V_uc, and the kernel physics is carried by +* the weight. When several processes describing parts of the same cross section are +* combined in one material, e.g. +* IncoherentElastic_process + IncoherentOnePhonon_process, +* their attenuations add up, so the attenuation is overestimated for thick samples. +* Single scattering intensities are not affected. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* Example: IncoherentOnePhonon_process(DOS_file="V-dos.dat", sigma_inc=5.08, +* unit_cell_volume=27.66, mass=50.94, T=300) +* +* %P +* INPUT PARAMETERS: +* DOS_file: [string] Phonon DOS, 2-column ASCII (E Z) or MCViNE IDF binary +* DOS_unit: [string] Energy unit of an ASCII DOS_file: "meV" or "THz". Empty: read from file comments, default meV +* sigma_inc: [barns] Incoherent scattering cross section of the unit cell +* unit_cell_volume: [AA^3] Unit cell volume +* mass: [amu] Average atomic mass +* T: [K] Temperature +* DW_core: [AA^2] Debye-Waller core, 2W = DW_core*Q^2. Negative: calculate from the DOS +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* Ef_focus: [meV] Centre of final energy window for energy focusing +* dEf_focus: [meV] Width of final energy window. 0: no energy focusing +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* IOP_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* %L +* J.Y.Y. Lin et al., Nucl. Instrum. Methods A 810 (2016) 86 +* %L +* G.L. Squires, Introduction to the Theory of Thermal Neutron Scattering, ch. 3 +* +* %E +******************************************************************************/ + +DEFINE COMPONENT IncoherentOnePhonon_process + +SETTING PARAMETERS(string DOS_file="", string DOS_unit="", double sigma_inc=-1, double unit_cell_volume=-1, double mass=-1, + double T=300, double DW_core=-1, double packing_factor=1, double Ef_focus=0, double dEf_focus=0, + double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this IncoherentOnePhonon_process component" + #endif + + %include "read_table-lib" + #include + + #ifndef INCOHERENTONEPHONON_SHARE + #define INCOHERENTONEPHONON_SHARE + + #define IOP_T2E (1.0 / 11.605) // Kelvin to meV, as in MCViNE + #define IOP_HBAR 1.05457148e-34 // J s + #define IOP_EV 1.60217653e-19 // J + #define IOP_AMU 1.66053886e-27 // kg + #define IOP_HZ2MEV (IOP_HBAR / (1e-3 * IOP_EV)) // angular frequency [rad/s] to meV + #define IOP_SMALL_OMEGA 1e-2 // |omega| < IOP_SMALL_OMEGA*Emax uses the parabolic DOS limit + + // Very important to add a pointer to this struct in the union-lib.c file + struct IncoherentOnePhonon_physics_storage_struct { + // DOS on a regular grid, normalised to unit area + int nE; + double E0, dE, Emax; + double* Z; + double sod; // Z(E) ~ sod*E^2 for small E + // material + double M; // average mass [amu] + double Temp; // [K] + double dw_core; // [AA^2] + double my_scattering; + // energy focusing + double Ef_c, dEf_c; + }; + + // linearly interpolated, normalised DOS, zero outside [E0, Emax) + double + iop_dos (struct IncoherentOnePhonon_physics_storage_struct* s, double e) { + if (e < s->E0 || e >= s->Emax) + return 0; + double x = (e - s->E0) / s->dE; + int j = (int)floor (x); + if (j >= s->nE - 1) + j = s->nE - 2; + return s->Z[j] + (x - j) * (s->Z[j + 1] - s->Z[j]); + } + + // Bose factor: n+1 for phonon creation (omega>0), n for annihilation + double + iop_bose (double omega, double T) { + if (omega == 0.0) + return 1.0; + double n = 1.0 / (exp (fabs (omega) / (T * IOP_T2E)) - 1.0); + return omega > 0 ? 1.0 + n : n; + } + + // fit y = c x, return c and R^2 (mccomponents/math/regression/linear1.h) + double + iop_linreg (double* x, double* y, int N, double* R2) { + double xys = 0, xxs = 0, ys = 0, sstot = 0, sserr = 0, c; + int i; + for (i = 0; i < N; i++) { + xys += x[i] * y[i]; + xxs += x[i] * x[i]; + ys += y[i]; + } + c = xys / xxs; + for (i = 0; i < N; i++) { + sstot += (y[i] - ys / N) * (y[i] - ys / N); + sserr += (y[i] - c * x[i]) * (y[i] - c * x[i]); + } + *R2 = sstot > 0 ? 1 - sserr / sstot : 1; + return c; + } + + // fit first N (100 down to 20) points to c*E^2 (utils.fitparabolic). Returns 1 if the fit is good + int + iop_fitparabolic (double* E, double* g, int n, int force, const char* comp) { + int N = n < 100 ? n : 100, minN = 20, bad = 1, i; + double c = 0, R2; + double* x = malloc (N * sizeof (double)); + while (N > minN) { + for (i = 0; i < N; i++) + x[i] = E[i] * E[i]; + c = iop_linreg (x, g, N, &R2); + if (R2 < 0.9) + N--; + else { + bad = 0; + break; + } + } + free (x); + if (bad && !force) + return 0; + if (bad) + fprintf (stderr, "IncoherentOnePhonon_process:%s: WARNING: unable to fit the low energy part of the DOS to a parabola\n", comp); + for (i = 0; i < N; i++) + g[i] = c * E[i] * E[i]; + return 1; + } + + // Hanning window smoothing with reflected ends, as utils.smooth(window_len=21) + void + iop_smooth (double* g, int n) { + int L = 21, i, k, ns = n + 2 * (L - 1); + double *s, *w, *y, wsum = 0; + if (n < L) + return; + s = malloc (ns * sizeof (double)); + w = malloc (L * sizeof (double)); + y = malloc ((ns - L + 1) * sizeof (double)); + for (i = 0; i < L - 1; i++) + s[i] = g[L - 1 - i]; // x[window_len-1:0:-1] + for (i = 0; i < n; i++) + s[L - 1 + i] = g[i]; + for (i = 0; i < L - 1; i++) + s[L - 1 + n + i] = g[n - 1 - i]; // x[-1:-window_len:-1] + for (i = 0; i < L; i++) { + w[i] = 0.5 - 0.5 * cos (2 * PI * i / (L - 1)); + wsum += w[i]; + } + for (i = 0; i < ns - L + 1; i++) { + y[i] = 0; + for (k = 0; k < L; k++) + y[i] += w[k] / wsum * s[i + k]; + } + for (i = 0; i < n; i++) // y[(L/2-1):-(L/2)-1] + g[i] = y[L / 2 - 1 + i]; + free (s); + free (w); + free (y); + } + + // Read DOS (ASCII or IDF), convert to meV, prepare as MCViNE, store in s + void + iop_load_dos (const char* file, const char* unit, struct IncoherentOnePhonon_physics_storage_struct* s, const char* comp) { + FILE* fp = Open_File ((char*)file, "rb", NULL); + double *E = NULL, *g = NULL, scale = 1; + int n = 0, i, cap = 0; + char head[64]; + if (!fp) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: cannot open DOS file %s\n", comp, file); + exit (-1); + } + if (fread (head, 1, 64, fp) == 64 && !strncmp (head, "DOS", 4)) { + // MCViNE IDF: 64 char type, int version, 1024 char comment, int N, double dE [THz], N doubles + int version, nb; + double dnu; + char comment[1024]; + if (fread (&version, sizeof (int), 1, fp) != 1 || fread (comment, 1, 1024, fp) != 1024 || fread (&nb, sizeof (int), 1, fp) != 1 + || fread (&dnu, sizeof (double), 1, fp) != 1 || nb < 2) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: corrupt IDF DOS file %s\n", comp, file); + exit (-1); + } + n = nb; + E = malloc (n * sizeof (double)); + g = malloc (n * sizeof (double)); + if (fread (g, sizeof (double), n, fp) != (size_t)n) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: IDF DOS file %s too short\n", comp, file); + exit (-1); + } + for (i = 0; i < n; i++) + E[i] = i * dnu; + scale = 2 * PI * 1e12 * IOP_HZ2MEV; + } else { + // ASCII + char line[4096], found_unit[16] = ""; + rewind (fp); + while (fgets (line, sizeof (line), fp)) { + char* c = line; + double x, y; + while (*c == ' ' || *c == '\t') + c++; + if (*c == '#') { + char tok[256] = ""; + if (!found_unit[0] && sscanf (c + 1, "%255s", tok) == 1) { + if (strstr (tok, "meV")) + strcpy (found_unit, "meV"); + else if (strstr (tok, "TeraHz") || strstr (tok, "THz")) + strcpy (found_unit, "THz"); + } + continue; + } + if (sscanf (c, "%lf %lf", &x, &y) != 2) + continue; + if (n == cap) { + cap = cap ? 2 * cap : 256; + E = realloc (E, cap * sizeof (double)); + g = realloc (g, cap * sizeof (double)); + } + E[n] = x; + g[n] = y; + n++; + } + if (unit && strlen (unit)) { + if (!strcmp (unit, "meV")) + scale = 1; + else if (!strcmp (unit, "THz") || !strcmp (unit, "TeraHz")) + scale = 2 * PI * 1e12 * IOP_HZ2MEV; + else { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: DOS_unit must be \"meV\" or \"THz\"\n", comp); + exit (-1); + } + } else if (!strcmp (found_unit, "THz")) + scale = 2 * PI * 1e12 * IOP_HZ2MEV; + else { + if (!found_unit[0]) + MPI_MASTER (printf ("IncoherentOnePhonon_process:%s: no energy unit found in %s, assuming meV\n", comp, file);); + scale = 1; + } + } + fclose (fp); + if (n < 2) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: DOS file %s has fewer than 2 points\n", comp, file); + exit (-1); + } + for (i = 0; i < n; i++) { + E[i] *= scale; + if (i && E[i] <= E[i - 1]) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: DOS energies must be strictly ascending\n", comp); + exit (-1); + } + } + if (E[0] < 0) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: DOS energies must be positive\n", comp); + exit (-1); + } + + // nice_dos: resample to 500 points from E=0 if the curve is coarse; otherwise ensure an equidistant grid + int uniform = 1; + for (i = 2; i < n; i++) + if (fabs ((E[i] - E[i - 1]) - (E[1] - E[0])) > 1e-6 * (E[n - 1] - E[0])) + uniform = 0; + if (n < 500 || !uniform) { + // numpy: E1 = arange(0, E[-1], E[-1]/500.), which may give 501 points due to rounding + double e0 = n < 500 ? 0 : E[0]; + double de = n < 500 ? E[n - 1] / 500. : (E[n - 1] - E[0]) / (n - 1); + int n1 = n < 500 ? (int)ceil (E[n - 1] / de) : n; + double *E1 = malloc (n1 * sizeof (double)), *g1 = malloc (n1 * sizeof (double)); + int j = 0; + for (i = 0; i < n1; i++) { + double x = e0 + i * de; + E1[i] = x; + if (x <= E[0]) + g1[i] = g[0]; + else if (x >= E[n - 1]) + g1[i] = g[n - 1]; + else { + while (E[j + 1] < x) + j++; + g1[i] = g[j] + (x - E[j]) / (E[j + 1] - E[j]) * (g[j + 1] - g[j]); + } + } + free (E); + free (g); + E = E1; + g = g1; + n = n1; + } + if (!iop_fitparabolic (E, g, n, 0, comp)) { + iop_smooth (g, n); + g[0] = 0; + iop_fitparabolic (E, g, n, 1, comp); + } + double area = 0; + for (i = 0; i < n; i++) + area += g[i]; + area *= E[1] - E[0]; + if (area <= 0) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: DOS has zero area\n", comp); + exit (-1); + } + for (i = 0; i < n; i++) + g[i] /= area; + + s->nE = n; + s->E0 = E[0]; + s->dE = E[1] - E[0]; + s->Emax = E[0] + s->dE * (n - 1); + s->Z = g; + + // second order coefficient at E=0 from the first 20 points (LinearlyInterpolatedDOS::_compute_sod) + int Nfit = n < 20 ? n : 20; + double x[20], R2; + for (i = 0; i < Nfit; i++) + x[i] = (E[0] + s->dE * i) * (E[0] + s->dE * i); + s->sod = iop_linreg (x, g, Nfit, &R2); + if (R2 < .9) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: failed to fit the first %d points of the DOS to a parabola\n", comp, Nfit); + exit (-1); + } + free (E); + } + + // 2W/Q^2 [AA^2] from the DOS (DWFromDOS.icc, 100 sampling points) + double + iop_dw_core (struct IncoherentOnePhonon_physics_storage_struct* s, double avg_mass, double T, const char* comp) { + int nSample = 100, first = -1, i; + double emin = s->E0, emax = s->Emax; + double dw = (emax - emin) / (nSample - 1 + .00000001), core = 0; + double* f = calloc (nSample, sizeof (double)); + for (i = 0; i < nSample; i++) { + double w = dw * i + emin; + double z = iop_dos (s, w); + if (w < emax / nSample / 100.) + continue; + if (first == -1) + first = i; + f[i] = (2 / (exp (w / (T * IOP_T2E)) - 1) + 1) / w * z; + core += f[i] * (i == nSample - 1 ? 0.5 : 1); + } + if (first == -1 || first + 1 >= nSample) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: invalid DOS (no data for E>0)\n", comp); + exit (-1); + } + double f0 = f[first] - (dw * first + emin) * (f[first + 1] - f[first]) / dw; + core += f0 / 2; + core /= IOP_EV * 1e-3; + core *= dw; + core *= IOP_HBAR * IOP_HBAR / 2 / IOP_AMU / avg_mass; + core *= 1e20; + free (f); + return core; + } + + #endif // INCOHERENTONEPHONON_SHARE + + // Inverse penetration depth: constant sigma_inc/V_uc as in MCViNE + int + IncoherentOnePhonon_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_IncoherentOnePhonon_physics_storage_struct->my_scattering; + return 1; + } + + // IncoherentInelastic::S / IncoherentInelastic_EnergyFocusing::S + int + IncoherentOnePhonon_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + struct IncoherentOnePhonon_physics_storage_struct* s = data_transfer.pointer_to_a_IncoherentOnePhonon_physics_storage_struct; + double vi[3], vi_l, Ei, Ef, e_range, solid_angle; + int d; + Coords dir; + + for (d = 0; d < 3; d++) + vi[d] = K2V * k_initial[d]; + vi_l = sqrt (vi[0] * vi[0] + vi[1] * vi[1] + vi[2] * vi[2]); + Ei = VS2E * vi_l * vi_l; + + // final energy + if (s->dEf_c > 0) { + double lo = s->Ef_c - s->dEf_c / 2, hi = s->Ef_c + s->dEf_c / 2; + if (lo < Ei - s->Emax) + lo = Ei - s->Emax; + if (lo < 0) + lo = 0; + if (hi > Ei + s->Emax) + hi = Ei + s->Emax; + if (hi <= lo) + return 0; + e_range = hi - lo; + Ef = lo + rand01 () * e_range; + } else if (Ei > s->Emax) { + e_range = 2 * s->Emax; + Ef = Ei - s->Emax + rand01 () * e_range; + } else { + e_range = Ei + s->Emax; + Ef = rand01 () * e_range; + } + if (Ef <= 0) + return 0; + double omega = Ei - Ef; + double vf_l = sqrt (Ef / VS2E); + + // final direction + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + + double Q2 = 0; + k_final[0] = V2K * vf_l * dir.x; + k_final[1] = V2K * vf_l * dir.y; + k_final[2] = V2K * vf_l * dir.z; + for (d = 0; d < 3; d++) + Q2 += (k_initial[d] - k_final[d]) * (k_initial[d] - k_final[d]); + double EQ = VS2E * K2V * K2V * Q2; // hbar^2 Q^2 / 2m_n [meV] + + double beta = 1. / (s->Temp * IOP_T2E); + double p = e_range / s->M * vf_l / vi_l * exp (-s->dw_core * Q2) * solid_angle / (4 * PI); + if (fabs (omega) < IOP_SMALL_OMEGA * s->Emax) + // omega -> 0: Bose factor -> 1/(beta omega), DOS -> sod*omega^2 + p *= s->sod / beta * EQ; + else + p *= iop_bose (omega, s->Temp) * iop_dos (s, fabs (omega)) * EQ / fabs (omega); + + *weight *= p; + return 1; + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_INCOHERENTONEPHONON_DETECTOR + #define PROCESS_INCOHERENTONEPHONON_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct IncoherentOnePhonon_physics_storage_struct IOP_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + memset (&IOP_storage, 0, sizeof (IOP_storage)); + + if (!strlen (DOS_file)) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: DOS_file must be given\n", NAME_CURRENT_COMP); + exit (-1); + } + if (sigma_inc <= 0 || unit_cell_volume <= 0 || mass <= 0 || T <= 0) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: sigma_inc, unit_cell_volume, mass and T must be positive\n", NAME_CURRENT_COMP); + exit (-1); + } + if (dEf_focus < 0 || (dEf_focus > 0 && Ef_focus <= 0)) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: ERROR: energy focusing needs Ef_focus > 0 and dEf_focus > 0\n", NAME_CURRENT_COMP); + exit (-1); + } + + iop_load_dos (DOS_file, DOS_unit, &IOP_storage, NAME_CURRENT_COMP); + IOP_storage.M = mass; + IOP_storage.Temp = T; + IOP_storage.dw_core = DW_core >= 0 ? DW_core : iop_dw_core (&IOP_storage, mass, T, NAME_CURRENT_COMP); + IOP_storage.my_scattering = packing_factor * sigma_inc / unit_cell_volume * 100; // barn/AA^3 -> 1/m + IOP_storage.Ef_c = Ef_focus; + IOP_storage.dEf_c = dEf_focus; + + MPI_MASTER (printf ("IncoherentOnePhonon_process:%s: DOS %d points up to %g meV, DW_core=%g AA^2, my_scattering=%g 1/m\n", NAME_CURRENT_COMP, IOP_storage.nE, + IOP_storage.Emax, IOP_storage.dw_core, IOP_storage.my_scattering);); + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Isotropic process (powder / polycrystal) + This_process.non_isotropic_rot_index = -1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = IncoherentOnePhonon; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_IncoherentOnePhonon_physics_storage_struct = &IOP_storage; + This_process.probability_for_scattering_function = &IncoherentOnePhonon_physics_my; + This_process.scattering_function = &IncoherentOnePhonon_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "IncoherentOnePhonon_process:%s: Error identifying Union_init component, %s is not a known component name.\n", NAME_CURRENT_COMP, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +FINALLY +%{ + free (IOP_storage.Z); +%} + +END diff --git a/mcstas-comps/union/IsotropicSqw_process.comp b/mcstas-comps/union/IsotropicSqw_process.comp new file mode 100644 index 0000000000..e1b31a8ebd --- /dev/null +++ b/mcstas-comps/union/IsotropicSqw_process.comp @@ -0,0 +1,543 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernels by Jiao Lin +* Date: 24.09.2026 +* Origin: DTU Physics +* +* Isotropic scattering from a dynamical structure factor S(|Q|,E) given on a grid or as an expression +* +* %D +* Port of the MCViNE kernels mccomponents/lib/kernels/sample/SQEkernel (with GridSQE and +* SQE_fromexpression) and SQE_EnergyFocusing_Kernel (J.Y.Y. Lin et al., Caltech/ORNL) +* to a Union process. +* +* 1/sigma d2sigma/dOmega dEf = kf/ki S(Q,E) / (4 pi) +* with E = Ei-Ef the energy transfer [meV] and Q = |ki-kf| [AA^-1], i.e. the same +* convention as Isotropic_Sqw. S(Q,E) [1/meV] is taken from either +* - Sqw_file: a McStas .sqw file as used by Isotropic_Sqw (a block with the m Q values, +* a block with the n E values and a block with m rows of n S values, blocks separated +* by comment lines, see e.g. data/He4_liq_coh.sqw). S is interpolated bilinearly and +* is zero outside the grid. The file header values "sigma_coh"/"sigma_inc", "V_rho" +* and "Temperature" are used when sigma, unit_cell_volume or T are not given. +* - SQE_expr: an expression (tinyexpr syntax) in Q, E and the free parameters p1..p4, +* e.g. "exp(-(E-5*Q)^2/2)/sqrt(2*pi)". Qmin, Qmax, Emin and Emax must then be given. +* If T > 0 and the E axis of the file only contains E >= 0, the energy gain side is added +* using detailed balance, S(Q,-E) = exp(-E/kT) S(Q,E). +* As in MCViNE, S is used as given (norm=1). With norm=-1, S from a file is renormalised +* as in Isotropic_Sqw (default there), which makes intensities comparable to Isotropic_Sqw. +* +* Two sampling methods are available: +* method=0 (as MCViNE): E and Q are sampled uniformly within the kinematically allowed +* range and the S(Q,E) range, the azimuthal angle is random. Union focusing is not used. +* method=1 (focusing): the final direction is sampled with the Union focusing of the +* geometry and Ef uniformly within the allowed range. +* With dEf_focus > 0 the final energy is restricted to [Ef_focus-dEf_focus/2, Ef_focus+dEf_focus/2] +* (energy focusing, SQE_EnergyFocusing_Kernel). +* +* The inverse penetration depth is sigma/unit_cell_volume, so S(Q,E) should be normalised +* accordingly (as in MCViNE); the weight carries the S(Q,E) dependence. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* %P +* INPUT PARAMETERS: +* Sqw_file: [string] McStas .sqw file with S(Q,E) +* SQE_expr: [string] Expression for S(Q,E) [1/meV], used when no Sqw_file is given +* sigma: [barns] Scattering cross section of the unit cell. <=0: from the Sqw_file header +* unit_cell_volume: [AA^3] Unit cell volume. <=0: from V_rho in the Sqw_file header +* T: [K] Temperature for detailed balance. <0: from the Sqw_file header, 0: none +* Qmin: [AA^-1] Minimum Q. With a file: further restricts the grid range +* Qmax: [AA^-1] Maximum Q (<=Qmin: use the grid range) +* Emin: [meV] Minimum energy transfer +* Emax: [meV] Maximum energy transfer (<=Emin: use the grid range) +* method: [1] 0: MCViNE sampling in (Q,E) (no focusing), 1: focusing +* norm: [1] Multiplier for S(Q,E). -1 (Sqw_file only): normalise as Isotropic_Sqw with the sum rule +* int Q^2 S(Q) dQ = Qmax^3/3 - 2 pi^2 rho (coherent) or Qmax^3/3 (incoherent), rho=1/unit_cell_volume +* incoherent: [1] 1: the Sqw_file is an incoherent S(Q,E) (only used for norm=-1) +* Ef_focus: [meV] Centre of final energy window for energy focusing +* dEf_focus: [meV] Width of final energy window. 0: no energy focusing +* p1: [1] Free parameter available in SQE_expr +* p2: [1] Free parameter available in SQE_expr +* p3: [1] Free parameter available in SQE_expr +* p4: [1] Free parameter available in SQE_expr +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* SQW_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* +* %E +******************************************************************************/ + +DEFINE COMPONENT IsotropicSqw_process + +SETTING PARAMETERS(string Sqw_file="", string SQE_expr="", double sigma=-1, double unit_cell_volume=-1, double T=-1, + double Qmin=0, double Qmax=0, double Emin=0, double Emax=0, int method=1, + double Ef_focus=0, double dEf_focus=0, double norm=1, int incoherent=0, + double p1=0, double p2=0, double p3=0, double p4=0, + double packing_factor=1, double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this IsotropicSqw_process component" + #endif + + %include "read_table-lib" + %include "tinyexpr.h" + %include "tinyexpr.c" + + #ifndef ISOTROPICSQW_SHARE + #define ISOTROPICSQW_SHARE + #define SQW_T2E (1.0 / 11.605) // Kelvin to meV + + struct sqw_grid { + int nq, nw; + double *q, *w, *S; // S[iq*nw + iw] + }; + + // index i with a[i] <= x < a[i+1], -1 if outside + int + sqw_locate (double* a, int n, double x) { + int lo = 0, hi = n - 1; + if (x < a[0] || x > a[n - 1]) + return -1; + while (hi - lo > 1) { + int mid = (lo + hi) / 2; + if (a[mid] <= x) + lo = mid; + else + hi = mid; + } + return lo; + } + + double + sqw_grid_eval (struct sqw_grid* g, double Q, double E) { + int i = sqw_locate (g->q, g->nq, Q), j = sqw_locate (g->w, g->nw, E); + if (i < 0 || j < 0) + return 0; + double tq = (Q - g->q[i]) / (g->q[i + 1] - g->q[i]), tw = (E - g->w[j]) / (g->w[j + 1] - g->w[j]); + double s00 = g->S[i * g->nw + j], s01 = g->S[i * g->nw + j + 1], s10 = g->S[(i + 1) * g->nw + j], s11 = g->S[(i + 1) * g->nw + j + 1]; + return (1 - tq) * ((1 - tw) * s00 + tw * s01) + tq * ((1 - tw) * s10 + tw * s11); + } + + // Read a McStas .sqw file: blocks of numbers separated by comment lines + // (q values, w values, S matrix), and the header values sigma_coh, sigma_inc, V_rho, Temperature + void + sqw_read_file (const char* file, struct sqw_grid* g, double* sig, double* vrho, double* temp, const char* comp) { + FILE* fp = Open_File ((char*)file, "r", NULL); + char line[65536]; + int block = -1, in_data = 0, cap[3] = { 0, 0, 0 }, n[3] = { 0, 0, 0 }, nlines[3] = { 0, 0, 0 }, rowlen = -1; + double* v[3] = { NULL, NULL, NULL }; + double sig_coh = 0, sig_inc = 0; + if (!fp) { + fprintf (stderr, "%s: ERROR: cannot open Sqw_file %s\n", comp, file); + exit (-1); + } + while (fgets (line, sizeof (line), fp)) { + char* c = line; + while (*c == ' ' || *c == '\t') + c++; + if (*c == '#' || *c == '%') { + char key[64]; + double val; + if (sscanf (c + 1, " %63s %lf", key, &val) == 2) { + if (!strcmp (key, "sigma_coh")) + sig_coh = val; + else if (!strcmp (key, "sigma_inc")) + sig_inc = val; + else if (!strcmp (key, "V_rho")) + *vrho = val; + else if (!strcmp (key, "Temperature")) + *temp = val; + } + in_data = 0; + continue; + } + // numeric line? + char* end; + double x = strtod (c, &end); + if (end == c) + continue; + if (!in_data) { + block++; + in_data = 1; + } + int b = block < 2 ? block : 2; + int cnt = 0; + while (end != c) { + if (n[b] == cap[b]) { + cap[b] = cap[b] ? 2 * cap[b] : 1024; + v[b] = realloc (v[b], cap[b] * sizeof (double)); + } + v[b][n[b]++] = x; + cnt++; + c = end; + x = strtod (c, &end); + } + if (b == 2 && rowlen < 0) + rowlen = cnt; + nlines[b]++; + } + fclose (fp); + if (block < 2) { + fprintf (stderr, "%s: ERROR: %s must contain three blocks of numbers (q, w, S) separated by comment lines\n", comp, file); + exit (-1); + } + g->nq = n[0]; + g->nw = n[1]; + g->q = v[0]; + g->w = v[1]; + g->S = v[2]; + if (n[2] != g->nq * g->nw) { + // maybe stored transposed (n rows of m values) + fprintf (stderr, "%s: ERROR: %s: S matrix has %d values, expected %d x %d\n", comp, file, n[2], g->nq, g->nw); + exit (-1); + } + if (g->nq < 2 || g->nw < 2) { + fprintf (stderr, "%s: ERROR: %s: need at least 2 q and 2 w values\n", comp, file); + exit (-1); + } + int i; + for (i = 1; i < g->nq; i++) + if (g->q[i] <= g->q[i - 1]) { + fprintf (stderr, "%s: ERROR: %s: q values must be ascending\n", comp, file); + exit (-1); + } + for (i = 1; i < g->nw; i++) + if (g->w[i] <= g->w[i - 1]) { + fprintf (stderr, "%s: ERROR: %s: w values must be ascending\n", comp, file); + exit (-1); + } + *sig = sig_coh > 0 ? sig_coh : sig_inc; + } + + // add the energy gain side by detailed balance S(Q,-E) = exp(-E/kT) S(Q,E) + void + sqw_detailed_balance (struct sqw_grid* g, double T) { + int skip0 = g->w[0] == 0 ? 1 : 0, nneg = g->nw - skip0, nw2 = g->nw + nneg, i, j; + double* w2 = malloc (nw2 * sizeof (double)); + double* S2 = malloc ((long)g->nq * nw2 * sizeof (double)); + for (j = 0; j < nneg; j++) + w2[j] = -g->w[g->nw - 1 - j]; + for (j = 0; j < g->nw; j++) + w2[nneg + j] = g->w[j]; + for (i = 0; i < g->nq; i++) { + for (j = 0; j < nneg; j++) { + double e = g->w[g->nw - 1 - j]; + S2[(long)i * nw2 + j] = exp (-e / (T * SQW_T2E)) * g->S[(long)i * g->nw + g->nw - 1 - j]; + } + for (j = 0; j < g->nw; j++) + S2[(long)i * nw2 + nneg + j] = g->S[(long)i * g->nw + j]; + } + free (g->w); + free (g->S); + g->w = w2; + g->S = S2; + g->nw = nw2; + } + #endif + + // Very important to add a pointer to this struct in the union-lib.c file + struct IsotropicSqw_physics_storage_struct { + struct sqw_grid grid; + int use_expr; + te_expr* expr; + double vQ, vE, vp[4]; + double qmin, qmax, emin, emax; + int meth; + double Ef_c, dEf_c; + double my_scattering; + double scale; + }; + + double + sqw_S (struct IsotropicSqw_physics_storage_struct* s, double Q, double E) { + if (Q < s->qmin || Q > s->qmax || E < s->emin || E > s->emax) + return 0; + if (s->use_expr) { + s->vQ = Q; + s->vE = E; + return s->scale * te_eval (s->expr); + } + return sqw_grid_eval (&s->grid, Q, E); + } + + int + IsotropicSqw_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_IsotropicSqw_physics_storage_struct->my_scattering; + return 1; + } + + // SQEkernel::S / SQE_EnergyFocusing_Kernel::S + int + IsotropicSqw_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + struct IsotropicSqw_physics_storage_struct* s = data_transfer.pointer_to_a_IsotropicSqw_physics_storage_struct; + double C = VS2E * K2V * K2V, ki, Ei, Efmin, Efmax; + int d; + + ki = sqrt (k_initial[0] * k_initial[0] + k_initial[1] * k_initial[1] + k_initial[2] * k_initial[2]); + Ei = C * ki * ki; + // final energy range from the E range of S and the energy window + Efmin = Ei - s->emax; + Efmax = Ei - s->emin; + if (Efmin < 0) + Efmin = 0; + if (s->dEf_c > 0) { + if (Efmin < s->Ef_c - s->dEf_c / 2) + Efmin = s->Ef_c - s->dEf_c / 2; + if (Efmax > s->Ef_c + s->dEf_c / 2) + Efmax = s->Ef_c + s->dEf_c / 2; + } + if (Efmax <= Efmin) + return 0; + double Ef = Efmin + rand01 () * (Efmax - Efmin), E = Ei - Ef, kf = sqrt (Ef / C); + + if (s->meth == 0) { + // --- MCViNE: Q uniform in the allowed range, random azimuth --- + double Q1 = fabs (ki - kf) > s->qmin ? fabs (ki - kf) : s->qmin; + double Q2 = ki + kf < s->qmax ? ki + kf : s->qmax; + if (Q2 <= Q1) + return 0; + double Q = Q1 + rand01 () * (Q2 - Q1); + double cost = (kf * kf + ki * ki - Q * Q) / (2 * kf * ki); + if (cost > 1) + cost = 1; + if (cost < -1) + cost = -1; + double sint = sqrt (1 - cost * cost), phi = 2 * PI * rand01 (); + double e1[3] = { k_initial[0] / ki, k_initial[1] / ki, k_initial[2] / ki }, e2[3], e3[3], n; + if (fabs (e1[0]) > 1e-10 || fabs (e1[1]) > 1e-10) { + e2[0] = -e1[1]; + e2[1] = e1[0]; + e2[2] = 0; + } else { + e2[0] = 1; + e2[1] = 0; + e2[2] = 0; + } + n = sqrt (e2[0] * e2[0] + e2[1] * e2[1] + e2[2] * e2[2]); + for (d = 0; d < 3; d++) + e2[d] /= n; + e3[0] = e1[1] * e2[2] - e1[2] * e2[1]; + e3[1] = e1[2] * e2[0] - e1[0] * e2[2]; + e3[2] = e1[0] * e2[1] - e1[1] * e2[0]; + for (d = 0; d < 3; d++) + k_final[d] = kf * (sint * cos (phi) * e2[d] + sint * sin (phi) * e3[d] + cost * e1[d]); + *weight *= sqw_S (s, Q, E) * Q * (Q2 - Q1) * (Efmax - Efmin) / (2 * ki * ki); + return 1; + } + + // --- method 1: focused direction --- + double solid_angle, Q2 = 0; + Coords dir; + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + k_final[0] = kf * dir.x; + k_final[1] = kf * dir.y; + k_final[2] = kf * dir.z; + for (d = 0; d < 3; d++) + Q2 += (k_initial[d] - k_final[d]) * (k_initial[d] - k_final[d]); + *weight *= solid_angle / (4 * PI) * (Efmax - Efmin) * kf / ki * sqw_S (s, sqrt (Q2), E); + return 1; + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_ISOTROPICSQW_DETECTOR + #define PROCESS_ISOTROPICSQW_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct IsotropicSqw_physics_storage_struct SQW_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + char who[512]; + double sig = sigma, vuc = unit_cell_volume, temp = T, file_sig = 0, file_vrho = 0, file_T = 0; + snprintf (who, 512, "IsotropicSqw_process:%s", NAME_CURRENT_COMP); + memset (&SQW_storage, 0, sizeof (SQW_storage)); + SQW_storage.scale = 1; + + if (strlen (Sqw_file) && strcmp (Sqw_file, "NULL") && strcmp (Sqw_file, "0")) { + sqw_read_file (Sqw_file, &SQW_storage.grid, &file_sig, &file_vrho, &file_T, who); + if (sig <= 0) + sig = file_sig; + if (vuc <= 0 && file_vrho > 0) + vuc = 1 / file_vrho; + if (temp < 0) + temp = file_T; + if (temp > 0 && SQW_storage.grid.w[0] >= 0) + sqw_detailed_balance (&SQW_storage.grid, temp); + SQW_storage.qmin = SQW_storage.grid.q[0]; + SQW_storage.qmax = SQW_storage.grid.q[SQW_storage.grid.nq - 1]; + SQW_storage.emin = SQW_storage.grid.w[0]; + SQW_storage.emax = SQW_storage.grid.w[SQW_storage.grid.nw - 1]; + if (Qmin > SQW_storage.qmin) + SQW_storage.qmin = Qmin; + if (Qmax > Qmin && Qmax < SQW_storage.qmax) + SQW_storage.qmax = Qmax; + if (Emax > Emin) { + if (Emin > SQW_storage.emin) + SQW_storage.emin = Emin; + if (Emax < SQW_storage.emax) + SQW_storage.emax = Emax; + } + } else if (strlen (SQE_expr)) { + int err = 0; + te_variable vars[] = { { "Q", &SQW_storage.vQ }, { "E", &SQW_storage.vE }, { "p1", &SQW_storage.vp[0] }, + { "p2", &SQW_storage.vp[1] }, { "p3", &SQW_storage.vp[2] }, { "p4", &SQW_storage.vp[3] } }; + SQW_storage.vp[0] = p1; + SQW_storage.vp[1] = p2; + SQW_storage.vp[2] = p3; + SQW_storage.vp[3] = p4; + SQW_storage.expr = te_compile (SQE_expr, vars, 6, &err); + if (!SQW_storage.expr) { + fprintf (stderr, "%s: ERROR: could not parse SQE_expr \"%s\" near character %d\n", who, SQE_expr, err); + exit (-1); + } + SQW_storage.use_expr = 1; + if (Qmax <= Qmin || Emax <= Emin || Qmin < 0) { + fprintf (stderr, "%s: ERROR: with SQE_expr, 0 <= Qmin < Qmax and Emin < Emax must be given\n", who); + exit (-1); + } + SQW_storage.qmin = Qmin; + SQW_storage.qmax = Qmax; + SQW_storage.emin = Emin; + SQW_storage.emax = Emax; + } else { + fprintf (stderr, "%s: ERROR: Sqw_file or SQE_expr must be given\n", who); + exit (-1); + } + if (norm < 0 && !SQW_storage.use_expr && vuc > 0) { + // sum rule normalisation as Isotropic_Sqw (H.E. Fischer et al, Rep. Prog. Phys. 69 (2006) 233, Eq 2.44) + struct sqw_grid* g = &SQW_storage.grid; + double iq2Sq = 0, alpha; + int i, j; + for (i = 0; i < g->nq; i++) { + double sq = 0, dq = (i == 0 ? g->q[1] - g->q[0] : (i == g->nq - 1 ? g->q[i] - g->q[i - 1] : (g->q[i + 1] - g->q[i - 1]) / 2)); + for (j = 0; j < g->nw; j++) { + double dw = (j == 0 ? g->w[1] - g->w[0] : (j == g->nw - 1 ? g->w[j] - g->w[j - 1] : (g->w[j + 1] - g->w[j - 1]) / 2)); + sq += g->S[(long)i * g->nw + j] * dw; + } + iq2Sq += g->q[i] * g->q[i] * sq * dq; + } + double qm = g->q[g->nq - 1]; + alpha = iq2Sq > 0 ? (qm * qm * qm / 3 - (incoherent ? 0 : 2 * PI * PI / vuc)) / iq2Sq : 0; + if (alpha <= 0) { + MPI_MASTER (printf ("%s: WARNING: normalisation factor %g is not positive, keeping S as given\n", who, alpha);); + alpha = 1; + } else + MPI_MASTER (printf ("%s: S(Q,E) normalised by the factor %g (sum rule)\n", who, alpha);); + for (i = 0; i < g->nq * g->nw; i++) + g->S[i] *= alpha; + } else if (norm > 0 && norm != 1) { + if (SQW_storage.use_expr) + SQW_storage.scale = norm; + else { + int i; + for (i = 0; i < SQW_storage.grid.nq * SQW_storage.grid.nw; i++) + SQW_storage.grid.S[i] *= norm; + } + } + if (sig <= 0 || vuc <= 0) { + fprintf (stderr, "%s: ERROR: sigma and unit_cell_volume must be positive (or given in the Sqw_file header)\n", who); + exit (-1); + } + if (method != 0 && method != 1) { + fprintf (stderr, "%s: ERROR: method must be 0 or 1\n", who); + exit (-1); + } + SQW_storage.meth = method; + SQW_storage.Ef_c = Ef_focus; + SQW_storage.dEf_c = dEf_focus; + SQW_storage.my_scattering = packing_factor * sig / vuc * 100; // barn/AA^3 -> 1/m + + MPI_MASTER (printf ("%s: Q=[%g,%g] AA^-1, E=[%g,%g] meV, sigma=%g barn, V=%g AA^3, my_scattering=%g 1/m%s\n", who, SQW_storage.qmin, SQW_storage.qmax, + SQW_storage.emin, SQW_storage.emax, sig, vuc, SQW_storage.my_scattering, + (!SQW_storage.use_expr && temp > 0) ? " (detailed balance applied if needed)" : "");); + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Isotropic process + This_process.non_isotropic_rot_index = -1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = IsotropicSqw; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_IsotropicSqw_physics_storage_struct = &SQW_storage; + This_process.probability_for_scattering_function = &IsotropicSqw_physics_my; + This_process.scattering_function = &IsotropicSqw_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "%s: Error identifying Union_init component, %s is not a known component name.\n", who, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +FINALLY +%{ + if (SQW_storage.use_expr) + te_free (SQW_storage.expr); + else { + free (SQW_storage.grid.q); + free (SQW_storage.grid.w); + free (SQW_storage.grid.S); + } +%} + +END diff --git a/mcstas-comps/union/Resolution_process.comp b/mcstas-comps/union/Resolution_process.comp new file mode 100644 index 0000000000..8a312ee06f --- /dev/null +++ b/mcstas-comps/union/Resolution_process.comp @@ -0,0 +1,260 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Peter Willendrup, port of the MCViNE kernels by Jiao Lin +* Date: 24.09.2026 +* Origin: DTU Physics +* +* Resolution / test scattering kernels: constant energy transfer, constant (|Q|,E), constant (Q vector,E), TOF window +* +* %D +* Port of the MCViNE kernels mccomponents/lib/kernels/sample/ConstantEnergyTransferKernel, +* ConstantQEKernel, ConstantvQEKernel and DGSSXResKernel (J.Y.Y. Lin et al., Caltech/ORNL) +* to a Union process. These kernels are mostly useful for resolution studies and tests, +* e.g. of an instrument with sample environment modelled with Union geometries. +* +* mode=0 (ConstantEnergyTransfer): Ef = Ei - E, direction sampled with the Union focusing +* of the geometry (isotropic in 4pi without focusing). Weight: solid angle/(4 pi). +* mode=1 (ConstantQE): Ef = Ei - E and |Q| = Q, i.e. scattering on a cone around ki with a +* random azimuthal angle (powder). No weight change; rays for which (Q,E) is +* kinematically impossible are absorbed. +* mode=2 (ConstantvQE): kf = ki - Q with the fixed vector Q = (Qx,Qy,Qz) in the local +* frame of the process (rotate the process to orient the "crystal"), weighted by a +* Gaussian exp(-(Ei-Ef-E)^2/(2 dE^2)) of the energy transfer. +* mode=3 (DGSSXRes): the final direction is sampled with the Union focusing of the +* geometry (use a target and e.g. focus_r) and the arrival time at the focusing target +* is sampled uniformly in [tof-dtof/2, tof+dtof/2]; the final speed follows from the +* distance to the target and the time of the scattering event. Weight: +* solid angle/(4 pi) * dtof * dEf/dt * kf/ki (the flight time from the scattering point +* is used in dEf/dt, MCViNE uses the total time of flight). +* +* The inverse penetration depth is sigma/unit_cell_volume. +* +* Part of the Union components, a set of components that work together and thus +* separates geometry and physics within McStas. +* The use of this component requires other components to be used. +* +* 1) One specifies a number of processes using process components like this one +* 2) These are gathered into material definitions using Union_make_material +* 3) Geometries are placed using Union_box / Union_cylinder, assigned a material +* 4) A Union_master component placed after all of the above +* +* Only in step 4 will any simulation happen, and per default all geometries +* defined before the master, but after the previous will be simulated here. +* +* Absorption is not handled by this process, set it in Union_make_material. +* +* %P +* INPUT PARAMETERS: +* mode: [1] 0: constant E, 1: constant |Q| and E, 2: constant Q vector and E, 3: TOF window at target +* E: [meV] Energy transfer (modes 0, 1, 2) +* Q: [AA^-1] Momentum transfer (mode 1) +* Qx: [AA^-1] Momentum transfer vector, x component (mode 2) +* Qy: [AA^-1] Momentum transfer vector, y component (mode 2) +* Qz: [AA^-1] Momentum transfer vector, z component (mode 2) +* dE: [meV] Standard deviation of the energy transfer (mode 2) +* tof: [s] Centre of the time of flight window at the focusing target (mode 3) +* dtof: [s] Width of the time of flight window (mode 3) +* sigma: [barns] Scattering cross section of the unit cell, defines the inverse penetration depth +* unit_cell_volume: [AA^3] Unit cell volume +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* init: [string] Name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* RES_storage: Storage struct transferred to Union_master +* +* %L +* MCViNE: https://github.com/mcvine/mcvine +* +* %E +******************************************************************************/ + +DEFINE COMPONENT Resolution_process + +SETTING PARAMETERS(int mode=0, double E=0, double Q=1, double Qx=0, double Qy=0, double Qz=1, double dE=1, + double tof=0.01, double dtof=1e-4, double sigma=5, double unit_cell_volume=13.8, double packing_factor=1, + double interact_fraction=-1, string init="init") + +SHARE +%{ + #ifndef Union + #error "The Union_init component must be included before this Resolution_process component" + #endif + + // Very important to add a pointer to this struct in the union-lib.c file + struct Resolution_physics_storage_struct { + int md; + double E0, Q0, Qv[3], sE, t0, dt; + double my_scattering; + }; + + int + Resolution_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + *my = data_transfer.pointer_to_a_Resolution_physics_storage_struct->my_scattering; + return 1; + } + + int + Resolution_physics_scattering (double* k_final, double* k_initial, double* weight, union data_transfer_union data_transfer, + struct focus_data_struct* focus_data, _class_particle* _particle) { + struct Resolution_physics_storage_struct* s = data_transfer.pointer_to_a_Resolution_physics_storage_struct; + double C = VS2E * K2V * K2V; + double ki = sqrt (k_initial[0] * k_initial[0] + k_initial[1] * k_initial[1] + k_initial[2] * k_initial[2]); + double Ei = C * ki * ki, Ef, kf, solid_angle; + Coords dir; + int d; + + switch (s->md) { + case 0: // ConstantEnergyTransferKernel + Ef = Ei - s->E0; + if (Ef <= 0) + return 0; + kf = sqrt (Ef / C); + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + k_final[0] = kf * dir.x; + k_final[1] = kf * dir.y; + k_final[2] = kf * dir.z; + *weight *= solid_angle / (4 * PI); + return 1; + case 1: { // ConstantQEKernel + Ef = Ei - s->E0; + if (Ef <= 0) + return 0; + kf = sqrt (Ef / C); + double cost = (ki * ki + kf * kf - s->Q0 * s->Q0) / (2 * ki * kf); + if (cost * cost > 1) + return 0; + double sint = sqrt (1 - cost * cost), phi = 2 * PI * rand01 (); + double e1[3] = { k_initial[0] / ki, k_initial[1] / ki, k_initial[2] / ki }, e2[3], e3[3], n; + if (fabs (e1[0]) > 1e-10 || fabs (e1[1]) > 1e-10) { + e2[0] = -e1[1]; + e2[1] = e1[0]; + e2[2] = 0; + } else { + e2[0] = 1; + e2[1] = 0; + e2[2] = 0; + } + n = sqrt (e2[0] * e2[0] + e2[1] * e2[1] + e2[2] * e2[2]); + for (d = 0; d < 3; d++) + e2[d] /= n; + e3[0] = e1[1] * e2[2] - e1[2] * e2[1]; + e3[1] = e1[2] * e2[0] - e1[0] * e2[2]; + e3[2] = e1[0] * e2[1] - e1[1] * e2[0]; + for (d = 0; d < 3; d++) + k_final[d] = kf * (sint * cos (phi) * e2[d] + sint * sin (phi) * e3[d] + cost * e1[d]); + return 1; + } + case 2: { // ConstantvQEKernel + double kf2 = 0; + for (d = 0; d < 3; d++) { + k_final[d] = k_initial[d] - s->Qv[d]; + kf2 += k_final[d] * k_final[d]; + } + double x = (Ei - C * kf2 - s->E0) / s->sE; + *weight *= exp (-x * x / 2); + return 1; + } + default: { // DGSSXResKernel + focus_data->focusing_function (&dir, &solid_angle, focus_data); + NORM (dir.x, dir.y, dir.z); + double L = sqrt (focus_data->RayAim.x * focus_data->RayAim.x + focus_data->RayAim.y * focus_data->RayAim.y + + focus_data->RayAim.z * focus_data->RayAim.z); + double tflight = s->t0 + (rand01 () - 0.5) * s->dt - _particle->t; + if (tflight <= 0 || L <= 0) + return 0; + double vf = L / tflight; + Ef = VS2E * vf * vf; + kf = V2K * vf; + k_final[0] = kf * dir.x; + k_final[1] = kf * dir.y; + k_final[2] = kf * dir.z; + *weight *= solid_angle / (4 * PI) * s->dt * (2 * Ef / tflight) * kf / ki; + return 1; + } + } + } + + // These lines help with future error correction, and tell other Union components + // that at least one process have been defined. + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_RESOLUTION_DETECTOR + #define PROCESS_RESOLUTION_DETECTOR dummy + #endif +%} + +DECLARE +%{ + struct Resolution_physics_storage_struct RES_storage; + + // Needed for transport to the main component, will be the same for all processes + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; +%} + +INITIALIZE +%{ + if (mode < 0 || mode > 3 || sigma <= 0 || unit_cell_volume <= 0 || (mode == 2 && dE <= 0) || (mode == 3 && (tof <= 0 || dtof <= 0))) { + fprintf (stderr, "Resolution_process:%s: ERROR: invalid parameters (mode 0-3, sigma, unit_cell_volume > 0, dE > 0 for mode 2, tof, dtof > 0 for mode 3)\n", + NAME_CURRENT_COMP); + exit (-1); + } + RES_storage.md = mode; + RES_storage.E0 = E; + RES_storage.Q0 = Q; + RES_storage.Qv[0] = Qx; + RES_storage.Qv[1] = Qy; + RES_storage.Qv[2] = Qz; + RES_storage.sE = dE; + RES_storage.t0 = tof; + RES_storage.dt = dtof; + RES_storage.my_scattering = packing_factor * sigma / unit_cell_volume * 100; // barn/AA^3 -> 1/m + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Only mode 2 depends on the orientation of the process + This_process.non_isotropic_rot_index = mode == 2 ? 1 : -1; + // Cross section does not depend on the focusing + This_process.needs_cross_section_focus = -1; + + // The type of the process must be saved in the global enum process + This_process.eProcess = Resolution; + + // Packing the data into a structure that is transported to the main component + This_process.data_transfer.pointer_to_a_Resolution_physics_storage_struct = &RES_storage; + This_process.probability_for_scattering_function = &Resolution_physics_my; + This_process.scattering_function = &Resolution_physics_scattering; + + // This will be the same for all process's, and can thus be moved to an include. + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + rot_copy (This_process.rotation_matrix, ROT_A_CURRENT_COMP); + sprintf (global_process_element.name, "%s", NAME_CURRENT_COMP); + global_process_element.component_index = INDEX_CURRENT_COMP; + global_process_element.p_scattering_process = &This_process; + + if (_getcomp_index (init) < 0) { + fprintf (stderr, "Resolution_process:%s: Error identifying Union_init component, %s is not a known component name.\n", NAME_CURRENT_COMP, init); + exit (-1); + } + + struct pointer_to_global_process_list* global_process_list = COMP_GETPAR3 (Union_init, init, global_process_list); + add_element_to_process_list (global_process_list, global_process_element); +%} + +TRACE +%{ + // Trace should be empty, the simulation is done in Union_master +%} + +END