diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.instr b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.instr new file mode 100755 index 000000000..d4e5af60e --- /dev/null +++ b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.instr @@ -0,0 +1,601 @@ +/******************************************************************************* +* McStas instrument definition URL=http://www.mcstas.org +* +* Instrument: Test_inhomogenous_process +* +* %Identification +* Written by: Daniel Lomholt Christensen +* Date: 18/03/2026 +* Origin: UCPH@NBI +* %INSTRUMENT_SITE: Tests_union +* +* Test of inhomogenous incoherent union component, +* and comparison to incoherent process, and incoherent component. +* +* %Description +* Test of 9 different setups. +* For all setups, the idea is the same. A source shines neutrons on the box, +* and a detector placed right up on the side of the box, measures the +* incoherently scattered neutrons. +* The first three have a constant scattering cross section throughout their +* depths: +* 1 is the normal incoherent component, +* 1 is a union box with an incoherent process +* 1 is a union box with the inhomogenous incoherent process +* Then the next three setups showcase a linearly rising scattering cross section +* in the first example, this linearly increasing cross section is implemented +* by 10 boxes utilizing the original inocherent process union component. +* +* Then the inhomogenous incoherent component implements it with its simple syntax +* +* Then a union box that uses the inhomogenous process +* and the original incoherent process is implemented. +* +* +* The 7th example examines the if the component works when a phonon component +* is present. It is a lot slower than the others, and is therefore run with +* fewer neutrons +* +* The 8th and 9th examples show an intuitive image of the transmitted beam. +* Here again the scattering cross section scales linearly, but due to the +* nature of the two different processes, the one process curve becomes +* continous (Inhomogenous_incoherent_process), +* and the other discrete (Incoherent_process). +* +* +* %Example: sample=thin Detector: lin_det_I=0.007387 +* %Example: sample=union_inc Detector: lin_det_I=0.007387 +* %Example: sample=union_inho Detector: lin_det_I=0.007387 +* %Example: sample=inc_linear Detector: lin_det_I=0.017457 +* %Example: sample=inho_linear Detector: lin_det_I=0.017457 +* %Example: sample=inc_and_inho Detector: lin_det_I=0.017457 +* %Example: -n 1e5 sample=phonon Detector: lin_det_I=0.0174609 +* %Example: sample=inc_trans Detector: lin_det_I=0.081029 +* %Example: sample=inho_trans Detector: lin_det_I=0.081040 +* +* +* +* +* %Parameters +* thick [m]: Thickness of box samples +* width [m]: Width of the box samples +* sample [string]: Parameter choosing the specific sample. see %Examples for the options +* +* +* %Link +* +* %End +*******************************************************************************/ +DEFINE INSTRUMENT Test_inhomogenous_process(thick=0.01, + width = 0.03, + string sample="thin") + +/* The DECLARE section allows us to declare variables or small */ +/* functions in C syntax. These may be used in the whole instrument. */ + +DECLARE +%{ + double E_i; + double eps; + int activate_inc; + int activate_inho; + int activate_inho_linear; + int activate_inc_linear; + int activate_inc_and_inho_linear; + int activate_phonon; + double rot_samples; + +%} + +/* The INITIALIZE section is executed when the simulation starts */ +/* (C code). You may use them as component parameter values. */ +INITIALIZE +%{ + eps=1e-6; + + E_i = 25.3; + + activate_inc = 0; + activate_inho = 0; + activate_inc_linear = 0; + activate_inho_linear = 0; + activate_inc_and_inho_linear = 0; + activate_phonon = 0; + rot_samples = 0; + + if (!strcmp("union_inc", sample)){ + activate_inc = 1; + } else if (!strcmp("union_inho", sample)){ + activate_inho = 1; + } else if (!strcmp("inho_linear", sample)){ + activate_inho_linear=1; + } else if (!strcmp("inc_linear", sample)){ + activate_inc_linear=1; + } else if (!strcmp("inc_and_inho", sample)){ + activate_inc_and_inho_linear = 1; + } else if (!strcmp("phonon", sample)){ + activate_phonon = 1; + } else if (!strcmp("inc_trans", sample)) { + rot_samples = 90; + activate_inc_linear = 1; + width = 0.01; + } else if (!strcmp("inho_trans", sample)){ + rot_samples = 90; + activate_inho_linear = 1; + width = 0.01; + } + +%} + +/* Here comes the TRACE section, where the actual */ +/* instrument is defined as a sequence of components. */ +TRACE + +COMPONENT origin = Progress_bar() + AT (0,0,0) ABSOLUTE + +COMPONENT source = Source_simple( + yheight = 0.035, + xwidth = 0.035, + dist = 10.5, + focus_xw = .03, + focus_yh = .03, + E0 = E_i, + dE = 0.1, + flux = 1e4) AT (0, 0, 0) RELATIVE origin +ROTATED (0, 0, 0) RELATIVE origin + +COMPONENT pre_samp_arm = Arm( + ) AT (0, 0, 10.475) RELATIVE source + +COMPONENT pre_samp_det = PSD_monitor( + nx = 200, + ny = 200, + filename="pre_det", + xwidth = 0.01, + yheight = 0.01, + restore_neutron = 1 +) AT (0, 0, -thick) RELATIVE pre_samp_arm + + +COMPONENT unrot_sample_arm = Arm( + +) +AT (0,0,0.025) RELATIVE pre_samp_arm + + +COMPONENT arm_sample = Arm( + +) AT (0, 0, 0) RELATIVE unrot_sample_arm +ROTATED (0, rot_samples, 0) RELATIVE unrot_sample_arm + +//============================================================================== +//=========================== START SAMPLE +//============================================================================== + +// =================== Simple Incoherent component +COMPONENT box_sample = Incoherent( + xwidth=width, + yheight=0.03, + zdepth=thick, + sigma_inc = 5.08, + sigma_abs = 5.08 +) WHEN (!strcmp((const char*)"thin", (const char*)sample)) +AT (0, 0, 0) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + + +COMPONENT init = Union_init() +AT (0, 0, 0) RELATIVE arm_sample +ROTATED (0.0, 0.0, 0.0) RELATIVE arm_sample + +// =================== Incoherent process constant + +COMPONENT inc_proc = Incoherent_process( + sigma = 5.08, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box = Union_box( + priority = 10, + xwidth=width, + yheight = 0.03, + zdepth = thick, + material_string = "Inco", + number_of_activations = activate_inc +) AT (0,0, 0) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + + +// =================== Inhomogenous_incoherent_process constant + +COMPONENT inho_proc = Inhomogenous_incoherent_process( + sigma_expr = "5.08", unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inho_inco = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inho_proc" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inho_box = Union_box( + priority = 10, + xwidth=width, + yheight = 0.03, + zdepth = thick, + material_string = "inho_inco", + number_of_activations = activate_inho +) AT (0,0, 0) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +// =================== Incoherent_process Linear +// Linear means a gradient from 5.08 to 50 over the 10 cm of the box +COMPONENT inc_proc_0 = Incoherent_process( + sigma = 5.08, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_0 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_0" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_0 = Union_box( + priority = 10, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string ="Inco_0", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +COMPONENT inc_proc_1 = Incoherent_process( + sigma = 5.08 + 9.58/2, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_1 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_1" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_1 = Union_box( + priority = 11, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_1", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20 + thick/10*1) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + + +COMPONENT inc_proc_2 = Incoherent_process( + sigma = 5.08 + 9.58/2*2, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_2 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_2" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_2 = Union_box( + priority = 12, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_2", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20 + thick/10*2) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +COMPONENT inc_proc_3 = Incoherent_process( + sigma = 5.08 + 9.58/2*3, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_3 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_3" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_3 = Union_box( + priority = 13, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_3", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20 + thick/10*3) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +COMPONENT inc_proc_4 = Incoherent_process( + sigma = 5.08 + 9.58/2*4, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_4 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_4" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_4 = Union_box( + priority = 14, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_4", + number_of_activations = activate_inc_linear +) AT (0,0,-thick/2 + thick/20 + thick/10*4 ) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +COMPONENT inc_proc_5 = Incoherent_process( + sigma = 5.08 + 9.58/2*5, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_5 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_5" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_5 = Union_box( + priority = 15, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_5", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20 + thick/10*5) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + + +COMPONENT inc_proc_6 = Incoherent_process( + sigma = 5.08 + 9.58/2*6, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_6 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_6" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_6 = Union_box( + priority = 16, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_6", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20 + thick/10*6) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +COMPONENT inc_proc_7 = Incoherent_process( + sigma = 5.08 + 9.58/2*7, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_7 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_7" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_7 = Union_box( + priority = 17, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_7", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20 + thick/10*7) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + + +COMPONENT inc_proc_8 = Incoherent_process( + sigma = 5.08 + 9.58/2*8, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_8 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_8" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_8 = Union_box( + priority = 18, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_8", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20 + thick/10*8) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +COMPONENT inc_proc_9 = Incoherent_process( + sigma = 5.08 + 9.58/2*9, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inco_9 = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inc_proc_9" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_box_9 = Union_box( + priority = 19, + xwidth=width, + yheight = 0.03, + zdepth = thick/10 - eps, + material_string = "Inco_9", + number_of_activations = activate_inc_linear +) AT (0,0, -thick/2 + thick/20 + thick/10*9) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + + + +// =================== Inhomogenous_incoherent_process Linear +COMPONENT inho_lin_proc = Inhomogenous_incoherent_process( + sigma_expr = "5.08 + 9.58/2*((z+0.0045)*1000)", unit_cell_volume = 13.827, + // verbose=1, + number_of_sample_points = 10 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Inho_lin = Union_make_material( + my_absorption = 100*5.08/13.827, + process_string="inho_lin_proc" +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inho_lin_box = Union_box( + priority = 10, + xwidth=width, + yheight = 0.03, + zdepth = thick, + material_string = "Inho_lin", + number_of_activations = activate_inho_linear +) AT (0,0, 0) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + + +// ============================ Inhomogenous + homogenous + +COMPONENT only_inc_proc = Incoherent_process( + sigma = 5.08, unit_cell_volume = 13.827 +) +AT (0,0,0) RELATIVE PREVIOUS + +COMPONENT only_inho_proc = Inhomogenous_incoherent_process( + sigma_expr = "9.58/2*((z+0.0045)*1000)", unit_cell_volume = 13.827, + // verbose=1, + number_of_sample_points = 10 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT inc_and_inho = Union_make_material( + my_absorption = 100 * 5.08 / 13.827, + process_string = "only_inc_proc,only_inho_proc" +) +AT (0,0,0) RELATIVE PREVIOUS + +COMPONENT inho_and_inc_box = Union_box( + priority = 10, + xwidth=width, + yheight = 0.03, + zdepth = thick, + material_string = "Inho_lin", + number_of_activations = activate_inc_and_inho_linear +) AT (0,0, 0) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + +// ============================ Inhomogenous + homogenous + phononsimple process + +COMPONENT phonon_proc = PhononSimple_process( + +) +AT (0,0,0) RELATIVE arm_sample + + + +COMPONENT phonon_mat = Union_make_material( + my_absorption = 100 * 5.08 / 13.827, + process_string = "only_inc_proc,only_inho_proc,phonon_proc" +) +AT (0,0,0) RELATIVE PREVIOUS + +COMPONENT phonon_box = Union_box( + priority = 10, + xwidth=width, + yheight = 0.03, + zdepth = thick, + material_string = "phonon_mat", + number_of_activations = activate_phonon +) AT (0,0, 0) RELATIVE arm_sample +ROTATED (0, 0, 0) RELATIVE arm_sample + + + + + + +// =================== Sample environment finalizer +COMPONENT space_logger = Union_logger_2D_space( + D_direction_1 = "x", + D_direction_2 = "z", + D1_min = -0.02, D1_max = 0.02, + D2_min = -0.01, D2_max = 0.01, + n1 = 360, + n2 = 360, + order_total = 1 +) +AT (0,0,0) RELATIVE arm_sample + +COMPONENT Sample_environment = Union_master( + verbal = 0) +AT (0, 0, 0) RELATIVE arm_sample +ROTATED (0.0, 0.0, 0.0) RELATIVE arm_sample + +COMPONENT stop = Union_stop() +AT (0, 0, 0) RELATIVE arm_sample +ROTATED (0.0, 0.0, 0.0) RELATIVE arm_sample + +//============================================================================== +//=========================== END SAMPLE +//============================================================================== + +COMPONENT arm_det = Arm( +) AT (0, 0, 0) RELATIVE unrot_sample_arm +ROTATED (0,90 - rot_samples,0) RELATIVE unrot_sample_arm + + + + + +COMPONENT psd_det = PSD_monitor( + nx = 200, + ny = 200, + filename="post_det", + xwidth = 0.02, + yheight = 0.1, + restore_neutron = 1 +) AT (0, 0, 0.0151) RELATIVE arm_det + +COMPONENT lin_det = PSD_monitor( + nx = 200, + ny = 1, + filename="lin_det", + xwidth = 0.02, + yheight = 0.1, + restore_neutron = 1 +) AT (0, 0, .0151) RELATIVE arm_det + + + +/* This section is executed when the simulation ends (C code). Other */ +/* optional sections are : SAVE */ +FINALLY +%{ +%} +/* The END token marks the instrument definition end */ +END + diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.md b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.md new file mode 100644 index 000000000..f51d61977 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Test_inhomogenous_process.md @@ -0,0 +1,35 @@ +# Test instrument for the Inhomogenous_incoherent component + +### Author: Daniel Lomholt Christensen @NBI + +This instrument is a test for the Inhomogenous incoherent process component, that works with the McStas Union system. + + +Below I have tested the implementation in three important setups. For each setup, I use three different approaches. +First I build the union materials with the Incoherent_process component (and here 10 Union_boxes are used). +Secondly I build the union materials with the Inhomogenous_incoherent_process (here using only 1 component, but setting the number of sampling points to 10). +Thirdly I build a material that uses both the Inhomogenous_incoherent and incoherent processes, to perform the same action. + + +The first two setups, I have a detector placed on the side of a sample, with 10 different boxes that follow each other. +For the first experiment, the boxes have the same scattering cross section. +For the second experiment, the boxes vary according to a linear rise in scattering cross section. +For a visualization of these setups, see the below figure, +![Sketch of instrument](experiment.png) + +And the results of the first experiment is, +![Results experiment 1](sigma_const.png) + +While the results of the second experiment is, +![Results experiment 2](linear_sigma.png) + +For the last setup, the detector is placed in a transmission geometry, with the 10 boxes rotated 90 degrees. +This means that they vary their scattering cross section linearly across x in the mcstas coordinate system. +The result of this is, +![Transmission geometry](trans_linear.png) + +Another experiment is performed, just to check that adding a directional process, such as the PhononSimple_process does not break the system. +This material is created with the linearly rising scattering cross section, implemented with the Inhomogenous process, as well as a PhononSimple_process. + + + diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/comparison.png b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/comparison.png new file mode 100644 index 000000000..1481efed0 Binary files /dev/null and b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/comparison.png differ diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/const_sigma.png b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/const_sigma.png new file mode 100644 index 000000000..bfd9844ef Binary files /dev/null and b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/const_sigma.png differ diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/experiment.png b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/experiment.png new file mode 100644 index 000000000..f1b40b241 Binary files /dev/null and b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/experiment.png differ diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/linear_sigma.png b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/linear_sigma.png new file mode 100644 index 000000000..c61eda1f4 Binary files /dev/null and b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/linear_sigma.png differ diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/trans_linear.png b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/trans_linear.png new file mode 100644 index 000000000..1481efed0 Binary files /dev/null and b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/trans_linear.png differ diff --git a/mcstas-comps/share/tinyexpr.c b/mcstas-comps/share/tinyexpr.c new file mode 100644 index 000000000..2a0ddfc4f --- /dev/null +++ b/mcstas-comps/share/tinyexpr.c @@ -0,0 +1,734 @@ +// SPDX-License-Identifier: Zlib +/* + * TINYEXPR - Tiny recursive descent parser and evaluation engine in C + * + * Copyright (c) 2015-2020 Lewis Van Winkle + * + * http://CodePlea.com + * + * This software is provided 'as-is', without any express or implied + * warranty. In no event will the authors be held liable for any damages + * arising from the use of this software. + * + * Permission is granted to anyone to use this software for any purpose, + * including commercial applications, and to alter it and redistribute it + * freely, subject to the following restrictions: + * + * 1. The origin of this software must not be misrepresented; you must not + * claim that you wrote the original software. If you use this software + * in a product, an acknowledgement in the product documentation would be + * appreciated but is not required. + * 2. Altered source versions must be plainly marked as such, and must not be + * misrepresented as being the original software. + * 3. This notice may not be removed or altered from any source distribution. + */ + +/* COMPILE TIME OPTIONS */ + +/* Exponentiation associativity: +For a^b^c = (a^b)^c and -a^b = (-a)^b do nothing. +For a^b^c = a^(b^c) and -a^b = -(a^b) uncomment the next line.*/ +/* #define TE_POW_FROM_RIGHT */ + +/* Logarithms +For log = base 10 log do nothing +For log = natural log uncomment the next line. */ +/* #define TE_NAT_LOG */ + +// #include "tinyexpr.h" +#include +#include +#include +#include +#include +#include + +#ifndef NAN +#define NAN (0.0/0.0) +#endif + +#ifndef INFINITY +#define INFINITY (1.0/0.0) +#endif + + +typedef double (*te_fun2)(double, double); + +enum { + TOK_NULL = TE_CLOSURE7+1, TOK_ERROR, TOK_END, TOK_SEP, + TOK_OPEN, TOK_CLOSE, TOK_NUMBER, TOK_VARIABLE, TOK_INFIX +}; + + +enum {TE_CONSTANT = 1}; + + +typedef struct state { + const char *start; + const char *next; + int type; + union {double value; const double *bound; const void *function;}; + void *context; + + const te_variable *lookup; + int lookup_len; +} state; + + +#define TYPE_MASK(TYPE) ((TYPE)&0x0000001F) + +#define IS_PURE(TYPE) (((TYPE) & TE_FLAG_PURE) != 0) +#define IS_FUNCTION(TYPE) (((TYPE) & TE_FUNCTION0) != 0) +#define IS_CLOSURE(TYPE) (((TYPE) & TE_CLOSURE0) != 0) +#define ARITY(TYPE) ( ((TYPE) & (TE_FUNCTION0 | TE_CLOSURE0)) ? ((TYPE) & 0x00000007) : 0 ) +#define NEW_EXPR(type, ...) new_expr((type), (const te_expr*[]){__VA_ARGS__}) +#define CHECK_NULL(ptr, ...) if ((ptr) == NULL) { __VA_ARGS__; return NULL; } + +static te_expr *new_expr(const int type, const te_expr *parameters[]) { + const int arity = ARITY(type); + const int psize = sizeof(void*) * arity; + const int size = (sizeof(te_expr) - sizeof(void*)) + psize + (IS_CLOSURE(type) ? sizeof(void*) : 0); + te_expr *ret = malloc(size); + CHECK_NULL(ret); + + memset(ret, 0, size); + if (arity && parameters) { + memcpy(ret->parameters, parameters, psize); + } + ret->type = type; + ret->bound = 0; + return ret; +} + + +void te_free_parameters(te_expr *n) { + if (!n) return; + switch (TYPE_MASK(n->type)) { + case TE_FUNCTION7: case TE_CLOSURE7: te_free(n->parameters[6]); /* Falls through. */ + case TE_FUNCTION6: case TE_CLOSURE6: te_free(n->parameters[5]); /* Falls through. */ + case TE_FUNCTION5: case TE_CLOSURE5: te_free(n->parameters[4]); /* Falls through. */ + case TE_FUNCTION4: case TE_CLOSURE4: te_free(n->parameters[3]); /* Falls through. */ + case TE_FUNCTION3: case TE_CLOSURE3: te_free(n->parameters[2]); /* Falls through. */ + case TE_FUNCTION2: case TE_CLOSURE2: te_free(n->parameters[1]); /* Falls through. */ + case TE_FUNCTION1: case TE_CLOSURE1: te_free(n->parameters[0]); + } +} + + +void te_free(te_expr *n) { + if (!n) return; + te_free_parameters(n); + free(n); +} + + +static double pi(void) {return 3.14159265358979323846;} +static double e(void) {return 2.71828182845904523536;} +static double fac(double a) {/* simplest version of fac */ + if (a < 0.0) + return NAN; + if (a > UINT_MAX) + return INFINITY; + unsigned int ua = (unsigned int)(a); + unsigned long int result = 1, i; + for (i = 1; i <= ua; i++) { + if (i > ULONG_MAX / result) + return INFINITY; + result *= i; + } + return (double)result; +} +static double ncr(double n, double r) { + if (n < 0.0 || r < 0.0 || n < r) return NAN; + if (n > UINT_MAX || r > UINT_MAX) return INFINITY; + unsigned long int un = (unsigned int)(n), ur = (unsigned int)(r), i; + unsigned long int result = 1; + if (ur > un / 2) ur = un - ur; + for (i = 1; i <= ur; i++) { + if (result > ULONG_MAX / (un - ur + i)) + return INFINITY; + result *= un - ur + i; + result /= i; + } + return result; +} +static double npr(double n, double r) {return ncr(n, r) * fac(r);} + +#ifdef _MSC_VER +#pragma function (ceil) +#pragma function (floor) +#endif + +static const te_variable functions[] = { + /* must be in alphabetical order */ + {"abs", fabs, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"acos", acos, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"asin", asin, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"atan", atan, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"atan2", atan2, TE_FUNCTION2 | TE_FLAG_PURE, 0}, + {"ceil", ceil, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"cos", cos, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"cosh", cosh, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"e", e, TE_FUNCTION0 | TE_FLAG_PURE, 0}, + {"exp", exp, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"fac", fac, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"floor", floor, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"ln", log, TE_FUNCTION1 | TE_FLAG_PURE, 0}, +#ifdef TE_NAT_LOG + {"log", log, TE_FUNCTION1 | TE_FLAG_PURE, 0}, +#else + {"log", log10, TE_FUNCTION1 | TE_FLAG_PURE, 0}, +#endif + {"log10", log10, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"ncr", ncr, TE_FUNCTION2 | TE_FLAG_PURE, 0}, + {"npr", npr, TE_FUNCTION2 | TE_FLAG_PURE, 0}, + {"pi", pi, TE_FUNCTION0 | TE_FLAG_PURE, 0}, + {"pow", pow, TE_FUNCTION2 | TE_FLAG_PURE, 0}, + {"sin", sin, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"sinh", sinh, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"sqrt", sqrt, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"tan", tan, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {"tanh", tanh, TE_FUNCTION1 | TE_FLAG_PURE, 0}, + {0, 0, 0, 0} +}; + +static const te_variable *find_builtin(const char *name, int len) { + int imin = 0; + int imax = sizeof(functions) / sizeof(te_variable) - 2; + + /*Binary search.*/ + while (imax >= imin) { + const int i = (imin + ((imax-imin)/2)); + int c = strncmp(name, functions[i].name, len); + if (!c) c = '\0' - functions[i].name[len]; + if (c == 0) { + return functions + i; + } else if (c > 0) { + imin = i + 1; + } else { + imax = i - 1; + } + } + + return 0; +} + +static const te_variable *find_lookup(const state *s, const char *name, int len) { + int iters; + const te_variable *var; + if (!s->lookup) return 0; + + for (var = s->lookup, iters = s->lookup_len; iters; ++var, --iters) { + if (strncmp(name, var->name, len) == 0 && var->name[len] == '\0') { + return var; + } + } + return 0; +} + + + +static double add(double a, double b) {return a + b;} +static double sub(double a, double b) {return a - b;} +static double mul(double a, double b) {return a * b;} +static double divide(double a, double b) {return a / b;} +static double negate(double a) {return -a;} +static double comma(double a, double b) {(void)a; return b;} + + +void next_token(state *s) { + s->type = TOK_NULL; + + do { + + if (!*s->next){ + s->type = TOK_END; + return; + } + + /* Try reading a number. */ + if ((s->next[0] >= '0' && s->next[0] <= '9') || s->next[0] == '.') { + s->value = strtod(s->next, (char**)&s->next); + s->type = TOK_NUMBER; + } else { + /* Look for a variable or builtin function call. */ + if (isalpha(s->next[0])) { + const char *start; + start = s->next; + while (isalpha(s->next[0]) || isdigit(s->next[0]) || (s->next[0] == '_')) s->next++; + + const te_variable *var = find_lookup(s, start, s->next - start); + if (!var) var = find_builtin(start, s->next - start); + + if (!var) { + s->type = TOK_ERROR; + } else { + switch(TYPE_MASK(var->type)) + { + case TE_VARIABLE: + s->type = TOK_VARIABLE; + s->bound = var->address; + break; + + case TE_CLOSURE0: case TE_CLOSURE1: case TE_CLOSURE2: case TE_CLOSURE3: /* Falls through. */ + case TE_CLOSURE4: case TE_CLOSURE5: case TE_CLOSURE6: case TE_CLOSURE7: /* Falls through. */ + s->context = var->context; /* Falls through. */ + + case TE_FUNCTION0: case TE_FUNCTION1: case TE_FUNCTION2: case TE_FUNCTION3: /* Falls through. */ + case TE_FUNCTION4: case TE_FUNCTION5: case TE_FUNCTION6: case TE_FUNCTION7: /* Falls through. */ + s->type = var->type; + s->function = var->address; + break; + } + } + + } else { + /* Look for an operator or special character. */ + switch (s->next++[0]) { + case '+': s->type = TOK_INFIX; s->function = add; break; + case '-': s->type = TOK_INFIX; s->function = sub; break; + case '*': s->type = TOK_INFIX; s->function = mul; break; + case '/': s->type = TOK_INFIX; s->function = divide; break; + case '^': s->type = TOK_INFIX; s->function = pow; break; + case '%': s->type = TOK_INFIX; s->function = fmod; break; + case '(': s->type = TOK_OPEN; break; + case ')': s->type = TOK_CLOSE; break; + case ',': s->type = TOK_SEP; break; + case ' ': case '\t': case '\n': case '\r': break; + default: s->type = TOK_ERROR; break; + } + } + } + } while (s->type == TOK_NULL); +} + + +static te_expr *list(state *s); +static te_expr *expr(state *s); +static te_expr *power(state *s); + +static te_expr *base(state *s) { + /* = | | {"(" ")"} | | "(" {"," } ")" | "(" ")" */ + te_expr *ret; + int arity; + + switch (TYPE_MASK(s->type)) { + case TOK_NUMBER: + ret = new_expr(TE_CONSTANT, 0); + CHECK_NULL(ret); + + ret->value = s->value; + next_token(s); + break; + + case TOK_VARIABLE: + ret = new_expr(TE_VARIABLE, 0); + CHECK_NULL(ret); + + ret->bound = s->bound; + next_token(s); + break; + + case TE_FUNCTION0: + case TE_CLOSURE0: + ret = new_expr(s->type, 0); + CHECK_NULL(ret); + + ret->function = s->function; + if (IS_CLOSURE(s->type)) ret->parameters[0] = s->context; + next_token(s); + if (s->type == TOK_OPEN) { + next_token(s); + if (s->type != TOK_CLOSE) { + s->type = TOK_ERROR; + } else { + next_token(s); + } + } + break; + + case TE_FUNCTION1: + case TE_CLOSURE1: + ret = new_expr(s->type, 0); + CHECK_NULL(ret); + + ret->function = s->function; + if (IS_CLOSURE(s->type)) ret->parameters[1] = s->context; + next_token(s); + ret->parameters[0] = power(s); + CHECK_NULL(ret->parameters[0], te_free(ret)); + break; + + case TE_FUNCTION2: case TE_FUNCTION3: case TE_FUNCTION4: + case TE_FUNCTION5: case TE_FUNCTION6: case TE_FUNCTION7: + case TE_CLOSURE2: case TE_CLOSURE3: case TE_CLOSURE4: + case TE_CLOSURE5: case TE_CLOSURE6: case TE_CLOSURE7: + arity = ARITY(s->type); + + ret = new_expr(s->type, 0); + CHECK_NULL(ret); + + ret->function = s->function; + if (IS_CLOSURE(s->type)) ret->parameters[arity] = s->context; + next_token(s); + + if (s->type != TOK_OPEN) { + s->type = TOK_ERROR; + } else { + int i; + for(i = 0; i < arity; i++) { + next_token(s); + ret->parameters[i] = expr(s); + CHECK_NULL(ret->parameters[i], te_free(ret)); + + if(s->type != TOK_SEP) { + break; + } + } + if(s->type != TOK_CLOSE || i != arity - 1) { + s->type = TOK_ERROR; + } else { + next_token(s); + } + } + + break; + + case TOK_OPEN: + next_token(s); + ret = list(s); + CHECK_NULL(ret); + + if (s->type != TOK_CLOSE) { + s->type = TOK_ERROR; + } else { + next_token(s); + } + break; + + default: + ret = new_expr(0, 0); + CHECK_NULL(ret); + + s->type = TOK_ERROR; + ret->value = NAN; + break; + } + + return ret; +} + + +static te_expr *power(state *s) { + /* = {("-" | "+")} */ + int sign = 1; + while (s->type == TOK_INFIX && (s->function == add || s->function == sub)) { + if (s->function == sub) sign = -sign; + next_token(s); + } + + te_expr *ret; + + if (sign == 1) { + ret = base(s); + } else { + te_expr *b = base(s); + CHECK_NULL(b); + + ret = NEW_EXPR(TE_FUNCTION1 | TE_FLAG_PURE, b); + CHECK_NULL(ret, te_free(b)); + + ret->function = negate; + } + + return ret; +} + +#ifdef TE_POW_FROM_RIGHT +static te_expr *factor(state *s) { + /* = {"^" } */ + te_expr *ret = power(s); + CHECK_NULL(ret); + + int neg = 0; + + if (ret->type == (TE_FUNCTION1 | TE_FLAG_PURE) && ret->function == negate) { + te_expr *se = ret->parameters[0]; + free(ret); + ret = se; + neg = 1; + } + + te_expr *insertion = 0; + + while (s->type == TOK_INFIX && (s->function == pow)) { + te_fun2 t = s->function; + next_token(s); + + if (insertion) { + /* Make exponentiation go right-to-left. */ + te_expr *p = power(s); + CHECK_NULL(p, te_free(ret)); + + te_expr *insert = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, insertion->parameters[1], p); + CHECK_NULL(insert, te_free(p), te_free(ret)); + + insert->function = t; + insertion->parameters[1] = insert; + insertion = insert; + } else { + te_expr *p = power(s); + CHECK_NULL(p, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, p); + CHECK_NULL(ret, te_free(p), te_free(prev)); + + ret->function = t; + insertion = ret; + } + } + + if (neg) { + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION1 | TE_FLAG_PURE, ret); + CHECK_NULL(ret, te_free(prev)); + + ret->function = negate; + } + + return ret; +} +#else +static te_expr *factor(state *s) { + /* = {"^" } */ + te_expr *ret = power(s); + CHECK_NULL(ret); + + while (s->type == TOK_INFIX && (s->function == pow)) { + te_fun2 t = s->function; + next_token(s); + te_expr *p = power(s); + CHECK_NULL(p, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, p); + CHECK_NULL(ret, te_free(p), te_free(prev)); + + ret->function = t; + } + + return ret; +} +#endif + + + +static te_expr *term(state *s) { + /* = {("*" | "/" | "%") } */ + te_expr *ret = factor(s); + CHECK_NULL(ret); + + while (s->type == TOK_INFIX && (s->function == mul || s->function == divide || s->function == fmod)) { + te_fun2 t = s->function; + next_token(s); + te_expr *f = factor(s); + CHECK_NULL(f, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, f); + CHECK_NULL(ret, te_free(f), te_free(prev)); + + ret->function = t; + } + + return ret; +} + + +static te_expr *expr(state *s) { + /* = {("+" | "-") } */ + te_expr *ret = term(s); + CHECK_NULL(ret); + + while (s->type == TOK_INFIX && (s->function == add || s->function == sub)) { + te_fun2 t = s->function; + next_token(s); + te_expr *te = term(s); + CHECK_NULL(te, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, te); + CHECK_NULL(ret, te_free(te), te_free(prev)); + + ret->function = t; + } + + return ret; +} + + +static te_expr *list(state *s) { + /* = {"," } */ + te_expr *ret = expr(s); + CHECK_NULL(ret); + + while (s->type == TOK_SEP) { + next_token(s); + te_expr *e = expr(s); + CHECK_NULL(e, te_free(ret)); + + te_expr *prev = ret; + ret = NEW_EXPR(TE_FUNCTION2 | TE_FLAG_PURE, ret, e); + CHECK_NULL(ret, te_free(e), te_free(prev)); + + ret->function = comma; + } + + return ret; +} + + +#define TE_FUN(...) ((double(*)(__VA_ARGS__))n->function) +#define M(e) te_eval(n->parameters[e]) + + +double te_eval(const te_expr *n) { + if (!n) return NAN; + + switch(TYPE_MASK(n->type)) { + case TE_CONSTANT: return n->value; + case TE_VARIABLE: return *n->bound; + + case TE_FUNCTION0: case TE_FUNCTION1: case TE_FUNCTION2: case TE_FUNCTION3: + case TE_FUNCTION4: case TE_FUNCTION5: case TE_FUNCTION6: case TE_FUNCTION7: + switch(ARITY(n->type)) { + case 0: return TE_FUN(void)(); + case 1: return TE_FUN(double)(M(0)); + case 2: return TE_FUN(double, double)(M(0), M(1)); + case 3: return TE_FUN(double, double, double)(M(0), M(1), M(2)); + case 4: return TE_FUN(double, double, double, double)(M(0), M(1), M(2), M(3)); + case 5: return TE_FUN(double, double, double, double, double)(M(0), M(1), M(2), M(3), M(4)); + case 6: return TE_FUN(double, double, double, double, double, double)(M(0), M(1), M(2), M(3), M(4), M(5)); + case 7: return TE_FUN(double, double, double, double, double, double, double)(M(0), M(1), M(2), M(3), M(4), M(5), M(6)); + default: return NAN; + } + + case TE_CLOSURE0: case TE_CLOSURE1: case TE_CLOSURE2: case TE_CLOSURE3: + case TE_CLOSURE4: case TE_CLOSURE5: case TE_CLOSURE6: case TE_CLOSURE7: + switch(ARITY(n->type)) { + case 0: return TE_FUN(void*)(n->parameters[0]); + case 1: return TE_FUN(void*, double)(n->parameters[1], M(0)); + case 2: return TE_FUN(void*, double, double)(n->parameters[2], M(0), M(1)); + case 3: return TE_FUN(void*, double, double, double)(n->parameters[3], M(0), M(1), M(2)); + case 4: return TE_FUN(void*, double, double, double, double)(n->parameters[4], M(0), M(1), M(2), M(3)); + case 5: return TE_FUN(void*, double, double, double, double, double)(n->parameters[5], M(0), M(1), M(2), M(3), M(4)); + case 6: return TE_FUN(void*, double, double, double, double, double, double)(n->parameters[6], M(0), M(1), M(2), M(3), M(4), M(5)); + case 7: return TE_FUN(void*, double, double, double, double, double, double, double)(n->parameters[7], M(0), M(1), M(2), M(3), M(4), M(5), M(6)); + default: return NAN; + } + + default: return NAN; + } + +} + +#undef TE_FUN +#undef M + +static void optimize(te_expr *n) { + /* Evaluates as much as possible. */ + if (n->type == TE_CONSTANT) return; + if (n->type == TE_VARIABLE) return; + + /* Only optimize out functions flagged as pure. */ + if (IS_PURE(n->type)) { + const int arity = ARITY(n->type); + int known = 1; + int i; + for (i = 0; i < arity; ++i) { + optimize(n->parameters[i]); + if (((te_expr*)(n->parameters[i]))->type != TE_CONSTANT) { + known = 0; + } + } + if (known) { + const double value = te_eval(n); + te_free_parameters(n); + n->type = TE_CONSTANT; + n->value = value; + } + } +} + + +te_expr *te_compile(const char *expression, const te_variable *variables, int var_count, int *error) { + state s; + s.start = s.next = expression; + s.lookup = variables; + s.lookup_len = var_count; + + next_token(&s); + te_expr *root = list(&s); + if (root == NULL) { + if (error) *error = -1; + return NULL; + } + + if (s.type != TOK_END) { + te_free(root); + if (error) { + *error = (s.next - s.start); + if (*error == 0) *error = 1; + } + return 0; + } else { + optimize(root); + if (error) *error = 0; + return root; + } +} + + +double te_interp(const char *expression, int *error) { + te_expr *n = te_compile(expression, 0, 0, error); + + double ret; + if (n) { + ret = te_eval(n); + te_free(n); + } else { + ret = NAN; + } + return ret; +} + +static void pn (const te_expr *n, int depth) { + int i, arity; + printf("%*s", depth, ""); + + switch(TYPE_MASK(n->type)) { + case TE_CONSTANT: printf("%f\n", n->value); break; + case TE_VARIABLE: printf("bound %p\n", n->bound); break; + + case TE_FUNCTION0: case TE_FUNCTION1: case TE_FUNCTION2: case TE_FUNCTION3: + case TE_FUNCTION4: case TE_FUNCTION5: case TE_FUNCTION6: case TE_FUNCTION7: + case TE_CLOSURE0: case TE_CLOSURE1: case TE_CLOSURE2: case TE_CLOSURE3: + case TE_CLOSURE4: case TE_CLOSURE5: case TE_CLOSURE6: case TE_CLOSURE7: + arity = ARITY(n->type); + printf("f%d", arity); + for(i = 0; i < arity; i++) { + printf(" %p", n->parameters[i]); + } + printf("\n"); + for(i = 0; i < arity; i++) { + pn(n->parameters[i], depth + 1); + } + break; + } +} + + +void te_print(const te_expr *n) { + pn(n, 0); +} diff --git a/mcstas-comps/share/tinyexpr.h b/mcstas-comps/share/tinyexpr.h new file mode 100644 index 000000000..c2cbe1a30 --- /dev/null +++ b/mcstas-comps/share/tinyexpr.h @@ -0,0 +1,87 @@ +// SPDX-License-Identifier: Zlib +/* + * TINYEXPR - Tiny recursive descent parser and evaluation engine in C + * + * Copyright (c) 2015-2020 Lewis Van Winkle + * + * http://CodePlea.com + * + * This software is provided 'as-is', without any express or implied + * warranty. In no event will the authors be held liable for any damages + * arising from the use of this software. + * + * Permission is granted to anyone to use this software for any purpose, + * including commercial applications, and to alter it and redistribute it + * freely, subject to the following restrictions: + * + * 1. The origin of this software must not be misrepresented; you must not + * claim that you wrote the original software. If you use this software + * in a product, an acknowledgement in the product documentation would be + * appreciated but is not required. + * 2. Altered source versions must be plainly marked as such, and must not be + * misrepresented as being the original software. + * 3. This notice may not be removed or altered from any source distribution. + */ + +#ifndef TINYEXPR_H +#define TINYEXPR_H + + +#ifdef __cplusplus +extern "C" { +#endif + + + +typedef struct te_expr { + int type; + union {double value; const double *bound; const void *function;}; + void *parameters[1]; +} te_expr; + + +enum { + TE_VARIABLE = 0, + + TE_FUNCTION0 = 8, TE_FUNCTION1, TE_FUNCTION2, TE_FUNCTION3, + TE_FUNCTION4, TE_FUNCTION5, TE_FUNCTION6, TE_FUNCTION7, + + TE_CLOSURE0 = 16, TE_CLOSURE1, TE_CLOSURE2, TE_CLOSURE3, + TE_CLOSURE4, TE_CLOSURE5, TE_CLOSURE6, TE_CLOSURE7, + + TE_FLAG_PURE = 32 +}; + +typedef struct te_variable { + const char *name; + const void *address; + int type; + void *context; +} te_variable; + + + +/* Parses the input expression, evaluates it, and frees it. */ +/* Returns NaN on error. */ +double te_interp(const char *expression, int *error); + +/* Parses the input expression and binds variables. */ +/* Returns NULL on error. */ +te_expr *te_compile(const char *expression, const te_variable *variables, int var_count, int *error); + +/* Evaluates the expression. */ +double te_eval(const te_expr *n); + +/* Prints debugging information on the syntax tree. */ +void te_print(const te_expr *n); + +/* Frees the expression. */ +/* This is safe to call on NULL pointers. */ +void te_free(te_expr *n); + + +#ifdef __cplusplus +} +#endif + +#endif /*TINYEXPR_H*/ diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index eac2f889a..755352aef 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -25,6 +25,7 @@ enum shape { }; enum process { + Inhomogenous_incoherent, Incoherent, Powder, Single_crystal, @@ -458,6 +459,7 @@ struct pointer_to_1d_int_list mask_intersect_list; int number_of_faces; struct surface_stack_struct **surface_stack_for_each_face; struct surface_stack_struct *internal_cut_surface_stack; + }; struct physics_struct @@ -475,10 +477,19 @@ struct scattering_process_struct *p_scattering_array; int has_refraction_info; double refraction_scattering_length_density; // [AA^-2] double refraction_Qc; + +// Numerical integration +int sampling_points; +double *cumul_transmission_prob; +double dist; +double *cumul_dists; +double **mus; +double *total_mus; }; union data_transfer_union{ // List of pointers to storage structs for all supported physical processes + struct Inhomogenous_incoherent_struct *Inhomogenous_incoherent_struct; struct Incoherent_physics_storage_struct *pointer_to_a_Incoherent_physics_storage_struct; struct Powder_physics_storage_struct *pointer_to_a_Powder_physics_storage_struct; struct Single_crystal_physics_storage_struct *pointer_to_a_Single_crystal_physics_storage_struct; @@ -499,36 +510,46 @@ union data_transfer_union{ struct scattering_process_struct { -char name[256]; // User defined process name -enum process eProcess; // enum value corresponding to this process GPU -double process_p_interact; // double between 0 and 1 that describes the fraction of events forced to undergo this process. -1 for disable -int non_isotropic_rot_index; // -1 if process is isotrpic, otherwise is the index of the process rotation matrix in the volume -int needs_cross_section_focus; // 1 if physics_my needs to call focus functions, otherwise -1 -Rotation rotation_matrix; // rotation matrix of process, reported by component in local frame, transformed and moved to volume struct in main - -union data_transfer_union data_transfer; // The way to reach the storage space allocated for this process (see examples in process.comp files) - -// probability_for_scattering_functions calculates this probability given k_i and parameters -int (*probability_for_scattering_function)(double*,double*,union data_transfer_union,struct focus_data_struct*, _class_particle *_particle); -// prop, k_i, ,parameters , focus data / function - -// A scattering_function takes k_i and parameters, returns k_f -int (*scattering_function)(double*,double*,double*,union data_transfer_union,struct focus_data_struct*, _class_particle *_particle); -// k_f, k_i, weight, parameters , focus data / function -}; - -//Utility function for initialising a scattering_process_struct with default -//values: -void scattering_process_struct_init( struct scattering_process_struct * sps ) + char name[256]; // User defined process name + enum process eProcess; // enum value corresponding to this process GPU + double process_p_interact; // double between 0 and 1 that describes the fraction of events forced to undergo this process. -1 for disable + int non_isotropic_rot_index; // -1 if process is isotrpic, otherwise is the index of the process rotation matrix in the volume + int needs_cross_section_focus; // 1 if physics_my needs to call focus functions, otherwise -1 + int needs_numerical_integration; // 1 if the process is inhomogenous and therefore needs numerical integration, otherwise -1. + Rotation rotation_matrix; // rotation matrix of process, reported by component in local frame, transformed and moved to volume struct in main + double *inhomogenous_cumul_prob; // The cumulative probability of a process in case of inhomogenous processes + double *inhomogenous_distances; // The distance of each step in which the cumulative probabilities will be calculated. + double *inhomogenous_cumul_distances; // The cumulative distances + double *inhomogenous_mu; // The different attenuation coefficients that are sampled in the numerical integration + double *inhomogenous_prob; // The probability of the process at the different sampled points. + double *inhomogenous_t; // The different times at which mu must be sampled. + int sampling_points; // Maximum number of samplings performed. If it is -1, no sampling has been done, and the arrays must be malloc'ed. + union data_transfer_union data_transfer; // The way to reach the storage space allocated for this process (see examples in process.comp files) + + // probability_for_scattering_functions calculates this probability given k_i and parameters + int (*probability_for_scattering_function)(double *, double *, union data_transfer_union, struct focus_data_struct *, _class_particle *_particle); + // prop, k_i, ,parameters , focus data / function + + // A scattering_function takes k_i and parameters, returns k_f + int (*scattering_function)(double *, double *, double *, union data_transfer_union, struct focus_data_struct *, _class_particle *_particle); + // k_f, k_i, weight, parameters , focus data / function +}; + +// Utility function for initialising a scattering_process_struct with default +// values: +void scattering_process_struct_init(struct scattering_process_struct *sps) { - memset(sps,0,sizeof(struct scattering_process_struct));//catch all + memset(sps, 0, sizeof(struct scattering_process_struct)); // catch all sps->name[0] = '\0'; sps->probability_for_scattering_function = NULL; sps->scattering_function = NULL; sps->non_isotropic_rot_index = -1; sps->needs_cross_section_focus = -1; + sps->needs_numerical_integration = -1; + sps->sampling_points = -1; } + union surface_data_transfer_union { struct Mirror_surface_storage_struct *pointer_to_a_Mirror_surface_storage_struct; diff --git a/mcstas-comps/share/union-suffix.c b/mcstas-comps/share/union-suffix.c index 1f7d2db17..a1bdb00b1 100644 --- a/mcstas-comps/share/union-suffix.c +++ b/mcstas-comps/share/union-suffix.c @@ -11,6 +11,11 @@ int physics_my(enum process choice, double *my,double *k_initial, union data_tra int output = 0; // Error return value #ifdef PROCESS_DETECTOR switch(choice) { + #ifdef PROCESS_INHOMOGENOUS_INCOHERENT_DETECTOR + case Inhomogenous_incoherent: + output = Inhomogenous_incoherent_physics_my(my, k_initial, data_transfer, focus_data, _particle); + break; + #endif #ifdef PROCESS_INCOHERENT_DETECTOR case Incoherent: output = Incoherent_physics_my(my, k_initial, data_transfer, focus_data, _particle); @@ -75,6 +80,11 @@ int physics_scattering(enum process choice, double *k_final, double *k_initial, int output = 0; // Error return value #ifdef PROCESS_DETECTOR switch(choice) { + #ifdef PROCESS_INHOMOGENOUS_INCOHERENT_DETECTOR + case Inhomogenous_incoherent: + output = Inhomogenous_incoherent_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); + break; + #endif #ifdef PROCESS_INCOHERENT_DETECTOR case Incoherent: output = Incoherent_physics_scattering(k_final, k_initial, weight, data_transfer, focus_data, _particle); diff --git a/mcstas-comps/union/Inhomogenous_incoherent_process.comp b/mcstas-comps/union/Inhomogenous_incoherent_process.comp new file mode 100755 index 000000000..3e581ff68 --- /dev/null +++ b/mcstas-comps/union/Inhomogenous_incoherent_process.comp @@ -0,0 +1,379 @@ +/******************************************************************************* +* +* McStas, neutron ray-tracing package +* Copyright(C) 2007 Risoe National Laboratory. +* +* %I +* Written by: Daniel Lomholt Christensen +* Date: 26/01/2026 +* Version: $Revision: 0.1 $ +* Origin: University of Copenhagen +* +* A sample component to separate geometry and phsysics +* +* %D +* +* This Union_process is based on the Incoherent_process.comp component +* originally written by Mads Bertelsen inspired by Kim Lefmann and +* Kristian Nielsen +* +* Part of the Union components, a set of components that work together and thus +* sperates 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. +* +* There is a dedicated manual available for the Union components +* +* Algorithm: +* The general algorithm for the Union system is described elsewhere. +* +* I here give a brief introduction as to what changes occur when using an +* inhomogenous process in your Union make material. It is expected that you +* understand the basic algorithm of the Union system before reading this. +* +* In Union, the neutron moves through a network of objects in a 3 dimensional world. +* When the neutron hits a material, the probability to scatter is calculated, +* and a Monte Carlo choice is taken, as to whether that neutron should scatter, +* or pass through. For a homogenous material (i.e constant attenuation coefficient $\mu$), +* this probability is the Beer-Lambert law, +* +*
+* $P_s = 1 - e^{-\mu l}$ +*
+* +* Where $P_s$ is the scattering probability, +* and $l$ is length of the neutron path throughout the object. +* +* +* For an inhomogenous material, this Beer-Lambert law must be modified, as $\mu$ is a function of the position. +* Therefore the Beer-Lambert law becomes, +* +*
+* $P_s = \int^l_0 1 - e^{-\mu(l')l'}dl'$ +*
+* +* Calculating this $\mu$ in the inhomogenous case is often trivial, but not feasible, +* from a software development point of view (seeing as many different functions of $\mu$ might be wanted). +* Instead the inhomogenous processes performs an approximate integral, by +* evaluating $\mu$ at a number of points along the neutron path (This number is in fact number_of_sample_points). +* +* For this incoherent process, the linear attenuation coefficient is, +* +*
+* $\mu = pack/V_u * 100 * \sigma$ +*
+* +* Where $pack$ is the packing factor of the material (defaults to 1), $V_u$ is the +* Unit cell volume, and $\sigma$ is the scattering cross section in barns. +* $\mu$ therefore has units of $m^{-1}$. +* +* For this component each factor in the attenuation coefficient can be a "tiny expression". +* This means that it can be a mathematical equation such as $\sigma_{expr} = "5.08 + 1000 * z * 2.35"$. +* When the attenuation coefficient is calculated, then the current value of $z$ is used to get $\sigma$. +* +* The parameters that the tiny expression can rely upon are currently: +* The positions, $x, y, z$ +* The velocities $vx, vy, vz$ +* and the time $t$ +* +* For more information on tiny expressions, see the link below. +* +* +* +* %P +* INPUT PARAMETERS: +* sigma: [barns]  Incoherent scattering cross section +* sigma_expr: [string] Tiny expression to be calculated as replacement for sigma +* packing_factor: [1] How dense is the material compared to optimal 0-1 +* packing_factor_expr: [string] Tiny expression to be calculated as replacement for packing factor +* unit_cell_volume: [AA^3] Unit cell volume +* unit_cell_volume_expr: [string] Tiny expression to be calculated as replacement for the unit cell volume +* gamma: [meV] Lorentzian width of quasielastic broadening (HWHM) [1] +* gamma_expr: [meV] Tiny expression to be calculated as replacement for the gamma value. +* f_QE: [1] Fraction of quasielastic scattering (rest is elastic) [1] +* number_of_sample_points [1] Number of points that are sampled along the neutron path through a material +* interact_fraction: [1] How large a part of the scattering events should use this process 0-1 (sum of all processes in material = 1) +* verbose [1] Flag that prints out the values calculated in the cross section calculation +* init: [string] name of Union_init component (typically "init", default) +* +* CALCULATED PARAMETERS: +* +* %L +* For information on how to write a tiny expression, see their github repository +* +* %E +******************************************************************************/ + +DEFINE COMPONENT Inhomogenous_incoherent_process + +SETTING PARAMETERS( sigma=0, string sigma_expr = "", + packing_factor=1, string packing_factor_expr = "", + unit_cell_volume=0, string unit_cell_volume_expr = "", + gamma=0, string gamma_expr = "", + f_QE=0, + number_of_sample_points = 20, + interact_fraction=-1, + int verbose = 0, + string init="init") + + +/* Neutron parameters: (x,y,z,vx,vy,vz,t,sx,sy,sz,p) */ + +SHARE +%{ + %include "tinyexpr.h" + %include "tinyexpr.c" + #ifndef Union + #error "The Union_init component must be included before this Incoherent_process component" + #endif + + struct Inhomogenous_incoherent_struct { + // Variables that needs to be transfered between any of the following places: + // The initialize in this component + // The function for calculating my + // The function for calculating scattering + double QE_sampling_frequency; + double lorentzian_width; + te_expr* lorentzian_width_expr; + double sig; + te_expr* sig_expr; + double unit_cell_vol; + te_expr* unit_cell_vol_expr; + double pack_fact; + te_expr* pack_fact_expr; + double* particle[7]; + int verbosity; + }; + + // Function for calculating my in Incoherent case + int + Inhomogenous_incoherent_physics_my (double* my, double* k_initial, union data_transfer_union data_transfer, struct focus_data_struct* focus_data, + _class_particle* _particle) { + double sigma; + double unit_cell_volume; + double packing_factor; + *data_transfer.Inhomogenous_incoherent_struct->particle[0] = _particle->x; + *data_transfer.Inhomogenous_incoherent_struct->particle[1] = _particle->y; + *data_transfer.Inhomogenous_incoherent_struct->particle[2] = _particle->z; + *data_transfer.Inhomogenous_incoherent_struct->particle[3] = _particle->vx; + *data_transfer.Inhomogenous_incoherent_struct->particle[4] = _particle->vy; + *data_transfer.Inhomogenous_incoherent_struct->particle[5] = _particle->vz; + *data_transfer.Inhomogenous_incoherent_struct->particle[6] = _particle->t; + + if (data_transfer.Inhomogenous_incoherent_struct->sig) { + sigma = data_transfer.Inhomogenous_incoherent_struct->sig; + } else { + // printf("\n Setting sigma with new values\npointer_z=%g\tpart_z=%g\tpointer=%p\n", + // *data_transfer.Inhomogenous_incoherent_struct->particle[2], _particle->z, data_transfer.Inhomogenous_incoherent_struct->particle[2]); + sigma = te_eval (data_transfer.Inhomogenous_incoherent_struct->sig_expr); + // printf("sigma=%g\n", sigma); + } + if (data_transfer.Inhomogenous_incoherent_struct->unit_cell_vol) { + unit_cell_volume = data_transfer.Inhomogenous_incoherent_struct->unit_cell_vol; + } else { + unit_cell_volume = te_eval (data_transfer.Inhomogenous_incoherent_struct->unit_cell_vol_expr); + } + if (data_transfer.Inhomogenous_incoherent_struct->pack_fact) { + packing_factor = data_transfer.Inhomogenous_incoherent_struct->pack_fact; + } else { + packing_factor = te_eval (data_transfer.Inhomogenous_incoherent_struct->pack_fact_expr); + } + + *my = ((packing_factor / unit_cell_volume) * 100 * sigma); + if (data_transfer.Inhomogenous_incoherent_struct->verbosity) { + printf ("\nDEBUG STATEMENT FOR INHOMOGENOUS INCOHERENT PROCESS:\n" + "mu=%g,\t sigma = %g,\t unit cell = %g,\t packing factor = %g\n" + "Neutron ray x,y,z = %g, %g, %g,\t Neutron speed = %g, %g, %g\t Neutron time = %g\n", + *my, sigma, unit_cell_volume, packing_factor, _particle->x, _particle->y, _particle->z, _particle->vx, _particle->vy, _particle->vz, _particle->t); + } + // printf("\nMu=%g\tsig=%g\tunit=%g\tpack=%g\n", *my, sigma, unit_cell_volume, packing_factor); + return 1; + }; + + // Function for basic incoherent scattering event + int + Inhomogenous_incoherent_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 lorentzian_width; + double QE_sampling_frequency = data_transfer.Inhomogenous_incoherent_struct->QE_sampling_frequency; + if (data_transfer.Inhomogenous_incoherent_struct->lorentzian_width) { + lorentzian_width = data_transfer.Inhomogenous_incoherent_struct->lorentzian_width; + } else { + lorentzian_width = te_eval (data_transfer.Inhomogenous_incoherent_struct->lorentzian_width_expr); + } + + // New version of incoherent scattering + double k_length = sqrt (k_initial[0] * k_initial[0] + k_initial[1] * k_initial[1] + k_initial[2] * k_initial[2]); + + Coords k_out; + // Here is the focusing system in action, get a vector + double solid_angle; + focus_data->focusing_function (&k_out, &solid_angle, focus_data); + NORM (k_out.x, k_out.y, k_out.z); + *weight *= solid_angle * 0.25 / PI; + + double v_i, v_f, E_i, dE, E_f; + + if (rand01 () < QE_sampling_frequency) { + v_i = k_length * K2V; + E_i = VS2E * v_i * v_i; + dE = lorentzian_width * tan (PI / 2 * randpm1 ()); + E_f = E_i + dE; + if (E_f <= 0) + return 0; + v_f = SE2V * sqrt (E_f); + k_length = v_f * V2K; + } + + k_final[0] = k_out.x * k_length; + k_final[1] = k_out.y * k_length; + k_final[2] = k_out.z * k_length; + return 1; + }; + + #ifndef PROCESS_DETECTOR + #define PROCESS_DETECTOR dummy + #endif + + #ifndef PROCESS_INHOMOGENOUS_INCOHERENT_DETECTOR + #define PROCESS_INHOMOGENOUS_INCOHERENT_DETECTOR dummy + #endif +%} + +DECLARE +%{ + // Needed for transport to the main component + struct global_process_element_struct global_process_element; + struct scattering_process_struct This_process; + + // Declare for this component, to do calculations on the input / store in the transported data + struct Inhomogenous_incoherent_struct Inhomogenous_storage; +%} + +INITIALIZE +%{ + // ========================================================================= + // ========================== Input sanitation ============================= + // ========================================================================= + if (!strcmp (sigma_expr, "") && !sigma) { + fprintf (stderr, "\nERROR! No sigma set, either through sigma_expr or sigma\tEXITING!\n\n"); + exit (1); + } + if (!strcmp (unit_cell_volume_expr, "") && !unit_cell_volume) { + fprintf (stderr, "\nERROR! No unit cell volume set, either through unit_cell_volume_expr or unit_cell_volume\tEXITING!\n\n"); + exit (1); + } + + // Store variable names and pointers. + Inhomogenous_storage.particle[0] = malloc (sizeof (double)); + Inhomogenous_storage.particle[1] = malloc (sizeof (double)); + Inhomogenous_storage.particle[2] = malloc (sizeof (double)); + Inhomogenous_storage.particle[3] = malloc (sizeof (double)); + Inhomogenous_storage.particle[4] = malloc (sizeof (double)); + Inhomogenous_storage.particle[5] = malloc (sizeof (double)); + Inhomogenous_storage.particle[6] = malloc (sizeof (double)); + te_variable vars[] = { { "x", Inhomogenous_storage.particle[0] }, { "y", Inhomogenous_storage.particle[1] }, { "z", Inhomogenous_storage.particle[2] }, + { "vx", Inhomogenous_storage.particle[3] }, { "vy", Inhomogenous_storage.particle[4] }, { "vz", Inhomogenous_storage.particle[5] }, + { "t", Inhomogenous_storage.particle[6] } }; + + // =========================================================================== + // ==================== Compile the expression with variables. =============== + // =========================================================================== + int err; + if (strcmp (sigma_expr, "")) { + te_expr* sigma_expression = te_compile (sigma_expr, vars, 7, &err); + if (!sigma_expression) { + printf ("Parse error at %d\n", err); + exit (1); + } + Inhomogenous_storage.sig_expr = sigma_expression; + } else + Inhomogenous_storage.sig = sigma; + if (strcmp (unit_cell_volume_expr, "")) { + te_expr* unit_cell_volume_expression = te_compile (unit_cell_volume_expr, vars, 7, &err); + if (!unit_cell_volume_expression) { + printf ("Parse error at %d\n", err); + exit (1); + } + Inhomogenous_storage.unit_cell_vol_expr = unit_cell_volume_expression; + } else + Inhomogenous_storage.unit_cell_vol = unit_cell_volume; + if (strcmp (packing_factor_expr, "")) { + te_expr* packing_fac_expression = te_compile (packing_factor_expr, vars, 7, &err); + if (!packing_fac_expression) { + printf ("Parse error at %d\n", err); + exit (1); + } + Inhomogenous_storage.pack_fact_expr = packing_fac_expression; + } else + Inhomogenous_storage.pack_fact = packing_factor; + if (strcmp (gamma_expr, "")) { + te_expr* gamma_expression = te_compile (gamma_expr, vars, 7, &err); + if (!gamma_expression) { + printf ("Parse error at %d\n", err); + exit (1); + } + Inhomogenous_storage.lorentzian_width_expr = gamma_expression; + } else + Inhomogenous_storage.lorentzian_width = gamma; + Inhomogenous_storage.QE_sampling_frequency = f_QE; + Inhomogenous_storage.verbosity = verbose; + + // =========================================================================== + // ==================== Perform Union specific dependencies ================= + // =========================================================================== + + // First initialise This_process with default values: + scattering_process_struct_init (&This_process); + + // Need to specify if this process is isotropic + This_process.non_isotropic_rot_index = -1; // Yes (powder) + // This_process.non_isotropic_rot_index = 1; // No (single crystal) + + // Need to specify if this process need to use focusing in calculation of inverse penetration depth (physics_my) + // This_process.needs_cross_section_focus = 1; // Yes + This_process.needs_cross_section_focus = -1; // No + + // The type of the process must be saved in the global enum process + This_process.eProcess = Inhomogenous_incoherent; + + // Packing the data into a structure that is transported to the main component + sprintf (This_process.name, "%s", NAME_CURRENT_COMP); + This_process.process_p_interact = interact_fraction; + This_process.data_transfer.Inhomogenous_incoherent_struct = &Inhomogenous_storage; + This_process.probability_for_scattering_function = &Inhomogenous_incoherent_physics_my; + This_process.scattering_function = &Inhomogenous_incoherent_physics_scattering; + This_process.needs_numerical_integration = 1; + This_process.sampling_points = number_of_sample_points; + + // This will be the same for all process's, and can thus be moved to an include. + 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, "Incoherent_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 +%{ +%} + +FINALLY +%{ + // Since the process and it's storage is a static allocation, there is nothing to deallocate +%} + +END diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 98dd5d6a5..ae0f31345 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -846,6 +846,31 @@ INITIALIZE } } // Initialization for each volume done + // ========================================================================= + // Loop over all physics structs, and assign arrays and values + // relevant to line integral approach + // ========================================================================= + + for (int mat_idx = 0; mat_idx < global_material_list_master->num_elements; mat_idx++) { + struct physics_struct* physics = global_material_list_master->elements[mat_idx].physics; + for (int proc_idx = 0; proc_idx < physics->number_of_processes; proc_idx++) { + struct scattering_process_struct spec_process = physics->p_scattering_array[proc_idx]; + if (spec_process.needs_numerical_integration != 1) + continue; + if (physics->sampling_points < spec_process.sampling_points) + physics->sampling_points = spec_process.sampling_points; + } + if (physics->sampling_points > 0) { + physics->cumul_transmission_prob = malloc (sizeof (double) * physics->sampling_points); + physics->cumul_dists = malloc (sizeof (double) * physics->sampling_points); + physics->mus = malloc (sizeof (double*) * physics->number_of_processes); + physics->total_mus = malloc (sizeof (double) * physics->sampling_points); + for (int i = 0; i < physics->number_of_processes; i++) { + physics->mus[i] = malloc (sizeof (double) * physics->sampling_points); + } + } + } + // ------- Initialization of ray-tracing algorithm ------------------------------------ my_trace = malloc (max_number_of_processes * sizeof (double)); @@ -1097,6 +1122,7 @@ TRACE double weight_limit; weight_limit = p * weight_ratio_limit; + double* sampling_mus; // Initialize logic done = 0; @@ -1537,6 +1563,8 @@ TRACE } } else { // Since there is a non-zero number of processes in this material, all the scattering cross section for these are calculated + struct physics_struct* current_p_physics = Volumes[current_volume]->p_physics; + int selected_sampling = -1; my_sum = 0; k[0] = V2K * vx; k[1] = V2K * vy; @@ -1549,7 +1577,7 @@ TRACE Coords ray_position_geometry; // If any process in this material needs focusing, sample scattering position and update focus_data accordingly - if (Volumes[current_volume]->p_physics->any_process_needs_cross_section_focus == 1) { + if (current_p_physics->any_process_needs_cross_section_focus == 1) { // Sample length_to_scattering in linear manner forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); @@ -1581,9 +1609,9 @@ TRACE // update focus data for this ray (could limit this to only update the necessary focus_data element, but there are typically very few) int f_index; for (f_index=0; f_index < Volumes[current_volume]->geometry.focus_data_array.num_elements; f_index++) { - this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[f_index]; - // Coords ray_position_geometry_rotated = rot_apply(this_focus_data.absolute_rotation, ray_position_geometry); - this_focus_data->RayAim = coords_sub(this_focus_data->Aim, ray_position_geometry); // Aim vector for this ray + this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[f_index]; + // Coords ray_position_geometry_rotated = rot_apply(this_focus_data.absolute_rotation, ray_position_geometry); + this_focus_data->RayAim = coords_sub(this_focus_data->Aim, ray_position_geometry); // Aim vector for this ray } printf("calculated forced_length_to_scattering = %lf, new RayAim \n", forced_length_to_scattering); @@ -1636,11 +1664,51 @@ TRACE process = &Volumes[current_volume]->p_physics->p_scattering_array[p_index]; // GPU Allowed int physics_output; - physics_output = physics_my (process->eProcess, p_my_trace, k_rotated, process->data_transfer, this_focus_data, _particle); + double mu; + current_p_physics->dist = length_to_boundary / current_p_physics->sampling_points - safety_distance; + + if (process->sampling_points != -1 || (process->needs_cross_section_focus == 1 && current_p_physics->sampling_points != 0)) { + // Populate length and probability arrays: + Coords original_position = coords_set (x, y, z); + for (int i = 0; i < current_p_physics->sampling_points; i++) { + current_p_physics->cumul_dists[i] = (i > 0) ? current_p_physics->cumul_dists[i - 1] + current_p_physics->dist : current_p_physics->dist / 2; + // Transport neutron to place inside geometry + ray_velocity = coords_set (vx, vy, vz); + // Find location of scattering point in master coordinate system without changing main position / velocity variables + Coords direction = coords_scalar_mult (ray_velocity, 1.0 / length_of_position_vector (ray_velocity)); + Coords sampling_displacement = coords_scalar_mult (direction, current_p_physics->cumul_dists[i]); + Coords sampling_point = coords_add (ray_position, sampling_displacement); + Coords sampling_point_geometry = coords_sub (sampling_point, Volumes[current_volume]->geometry.center); + // Also focus the ray at this point, if the component needs focusing + if (process->needs_cross_section_focus) { + this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; + if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { + int non_isotropic_rot_index = Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; + sampling_point_geometry + = rot_apply (Volumes[current_volume]->geometry.process_rot_matrix_array[non_isotropic_rot_index], sampling_point_geometry); + } + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, sampling_point_geometry); + } + // Calculate mu and probability + coords_get (sampling_point_geometry, &x, &y, &z); + physics_my (process->eProcess, &mu, k_rotated, process->data_transfer, this_focus_data, _particle); + current_p_physics->mus[p_index][i] = mu; + } + + coords_get (original_position, &x, &y, &z); + *p_my_trace = 0; + for (int i = 0; i < current_p_physics->sampling_points; i++) + *p_my_trace += current_p_physics->mus[p_index][i] / current_p_physics->sampling_points; + } else { + physics_output = physics_my (process->eProcess, p_my_trace, k_rotated, process->data_transfer, this_focus_data, _particle); + if (current_p_physics->sampling_points != 0) { + current_p_physics->mus[p_index][0] = *p_my_trace; + } + } my_sum += *p_my_trace; #ifdef Union_trace_verbal_setting - printf ("my_trace = %f, my_sum = %f\n", *p_my_trace, my_sum); + printf ("my_trace = %f, my_sum = %f\n", *p_my_trace, my_sum); #endif // increment the pointers so that it point to the next element (max number of process in any material is allocated) p_my_trace++; @@ -1663,7 +1731,42 @@ TRACE } else if (length_to_boundary < safety_distance2) { scattering_event = 0; } else { - real_transmission_probability = exp (-length_to_boundary * my_sum_plus_abs); + if (current_p_physics->sampling_points == 0) { + real_transmission_probability = exp (-length_to_boundary * my_sum_plus_abs); + } else if (current_p_physics->sampling_points != 0) { + // Calculate the probabilities and then add them cumulatively + memset (current_p_physics->total_mus, 0, sizeof (double) * current_p_physics->sampling_points); + + for (int i = 0; i < Volumes[current_volume]->p_physics->number_of_processes; i++) { + struct scattering_process_struct* process_i = &Volumes[current_volume]->p_physics->p_scattering_array[i]; + if (process_i->needs_numerical_integration != 1) + for (int j = 0; j < current_p_physics->sampling_points; j++) { + current_p_physics->total_mus[j] += current_p_physics->mus[i][0] * current_p_physics->dist; + } + else + for (int j = 0; j < current_p_physics->sampling_points; j++) { + current_p_physics->total_mus[j] += current_p_physics->mus[i][j] * current_p_physics->dist; + } + } + // for (int i =0;isampling_points;i++){ + // printf("\nTotalmu=%g\tinteger=%d\n", current_p_physics->total_mus[i], i); + // } + double mu_at_speed = Volumes[current_volume]->p_physics->my_a * (2200 / v_length); + for (int j = 0; j < current_p_physics->sampling_points; j++) { + current_p_physics->total_mus[j] += mu_at_speed * current_p_physics->dist; + } + double trans_prob; + for (int i = 0; i < current_p_physics->sampling_points; i++) { + trans_prob = exp (-current_p_physics->total_mus[i]); + if (i == 0) + current_p_physics->cumul_transmission_prob[i] = trans_prob; + else + current_p_physics->cumul_transmission_prob[i] = current_p_physics->cumul_transmission_prob[i - 1] * trans_prob; + } + + real_transmission_probability = current_p_physics->cumul_transmission_prob[current_p_physics->sampling_points - 1]; + } + // printf("Trans prop = %g\n", real_transmission_probability); if (Volumes[current_volume]->geometry.geometry_p_interact != 0) { mc_transmission_probability = (1.0 - Volumes[current_volume]->geometry.geometry_p_interact); if ((scattering_event = (rand01 () > mc_transmission_probability))) { @@ -1675,6 +1778,7 @@ TRACE } } else { // probability to scatter is the natural value + // printf("Real transmission prob %g\n", real_transmission_probability); scattering_event = rand01 () > real_transmission_probability; } } @@ -1689,6 +1793,31 @@ TRACE // Safety feature, alert in case of nonsense my results / negative absorption if (my_sum / my_sum_plus_abs > 1.0) printf ("WARNING: Absorption weight factor above 1! Should not happen! \n"); + // Select distance to scattering position + if (current_p_physics->sampling_points != 0) { + // Numerical integration happens, and therefore we must choose between the different samples + // We do this by drawing a random number between 0 and max cumul prob, + // and then seeing which cumul prob is the first to include it. + abs_weight_factor = 1; + double mu_at_speed = Volumes[current_volume]->p_physics->my_a * (2200 / v_length); + double pseudo_rand = rand01 () * (1 - current_p_physics->cumul_transmission_prob[current_p_physics->sampling_points - 1]); + for (int i = 0; i < current_p_physics->sampling_points; i++) { + // printf("\nCumul trans prob = %g\t pseudo rand = %g\n", current_p_physics->cumul_transmission_prob[i], pseudo_rand); + if (pseudo_rand >= 1 - current_p_physics->cumul_transmission_prob[i]) + continue; + selected_sampling = i; + break; + } + abs_weight_factor *= (current_p_physics->total_mus[selected_sampling] - mu_at_speed * current_p_physics->dist) / current_p_physics->total_mus[selected_sampling]; + + // printf("\nSelected_sampling = %d\n", selected_sampling); + + // printf("dist i = %g\tdist=%g\n", dist_i, dist); + double sampled_dist + = safety_distance + - log (1.0 - rand01 () * (1.0 - exp (-current_p_physics->total_mus[selected_sampling]))) / current_p_physics->total_mus[selected_sampling] * current_p_physics->dist; + length_to_scattering = current_p_physics->cumul_dists[selected_sampling] - current_p_physics->dist / 2 + sampled_dist; + } // Select process if (Volumes[current_volume]->p_physics->number_of_processes == 1) { // trivial case // Select the only available process, which will always have index 0 @@ -1724,42 +1853,56 @@ TRACE // Select a process based on their relative attenuations factors mc_prop = rand01 (); culmative_probability = 0; - for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { - culmative_probability += my_trace[iterator] / my_sum; - if (culmative_probability > mc_prop) { - selected_process = iterator; - break; + if (current_p_physics->sampling_points == 0) { + for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { + culmative_probability += my_trace[iterator] / my_sum; + if (culmative_probability > mc_prop) { + selected_process = iterator; + break; + } + } + } else { + for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { + culmative_probability += current_p_physics->mus[iterator][selected_sampling] / my_sum; + if (culmative_probability > mc_prop) { + selected_process = iterator; + break; + } } } } } - // Select distance to scattering position - if (Volumes[current_volume]->p_physics->p_scattering_array[selected_process].needs_cross_section_focus == 1) { - // Respect forced length to scattering chosen by process - length_to_scattering = forced_length_to_scattering; - // Drawing between 0 and L from constant s = 1/L and should have been q = A*exp(-kz). - // Normalizing A*exp(-kz) over 0 to L: A = k/(1-exp(-k*L)) - // Weight correction is ratio between s and q, L*A*exp(-kz) = L*k*exp(-kz)/(1-exp(-Lk)) - p *= length_to_boundary * my_sum_plus_abs * exp (-length_to_scattering * my_sum_plus_abs) / (1.0 - exp (-length_to_boundary * my_sum_plus_abs)); - #ifdef Union_trace_verbal_setting - printf ("Used forced length to scattering, %lf \n", length_to_scattering); - #endif + process = &Volumes[current_volume]->p_physics->p_scattering_array[selected_process]; + if (current_p_physics->sampling_points == 0) { + // No numerical integration is necessary. + if (process->needs_cross_section_focus == 1) { + // Respect forced length to scattering chosen by process + length_to_scattering = forced_length_to_scattering; + // Drawing between 0 and L from constant s = 1/L and should have been q = A*exp(-kz). + // Normalizing A*exp(-kz) over 0 to L: A = k/(1-exp(-k*L)) + // Weight correction is ratio between s and q, L*A*exp(-kz) = L*k*exp(-kz)/(1-exp(-Lk)) + p *= length_to_boundary * my_sum_plus_abs * exp (-length_to_scattering * my_sum_plus_abs) / (1.0 - exp (-length_to_boundary * my_sum_plus_abs)); + #ifdef Union_trace_verbal_setting + printf ("Used forced length to scattering, %lf \n", length_to_scattering); + #endif - } else { - // Decided the ray scatters, choose where on truncated exponential from safety_distance to length_to_boundary - safety_distance - length_to_scattering - = safety_distance - log (1.0 - rand0max ((1.0 - exp (-my_sum_plus_abs * (length_to_boundary - safety_distance2))))) / my_sum_plus_abs; - #ifdef Union_trace_verbal_setting - printf ("Sampled length to scattering, %lf \n", length_to_scattering); - #endif + } else { + // Decided the ray scatters, choose where on truncated exponential from safety_distance to length_to_boundary - safety_distance + length_to_scattering + = safety_distance - log (1.0 - rand0max (1.0 - exp (-my_sum_plus_abs * (length_to_boundary - safety_distance2)))) / my_sum_plus_abs; + #ifdef Union_trace_verbal_setting + printf ("Sampled length to scattering, %lf \n", length_to_scattering); + #endif + } } + } // Done handling sampling of position and process - } // Done handling scattering + } // Done handling scattering - } // Done handling situation where there are scattering processes in the material + } // Done handling situation where there are scattering processes in the material - } // Done checking for scttering event and in case of scattering selecting a process + // Done checking for scttering event and in case of scattering selecting a process // Record initial weight, absorption weight factor and initial position @@ -1769,6 +1912,9 @@ TRACE r_old[2] = r[2]; time_old = t; // Apply absorption + // if (abs_weight_factor != 1){ + // printf("Abs weight factor = %g\n", abs_weight_factor); + // } p *= abs_weight_factor; // Create event for absorption loggers @@ -1937,8 +2083,8 @@ TRACE printf("ERROR: Union_master: %s.Absorbed ray because scattering function returned 0 (error/absorb)\n",NAME_CURRENT_COMP); component_error_msg++; if (component_error_msg > 100) { - printf("To many errors encountered, exiting. \n"); - exit(1); + printf("To many errors encountered, exiting. \n"); + exit(1); } */ ABSORB; @@ -3001,7 +3147,7 @@ TRACE } else { ABSORB; // Absorb rays that didn't exit correctly for whatever reason - // Could error log here + // Could error log here } // Stores nubmer of scattering events in global master list so that another master with inherit_number_of_scattering_events can continue