} */
+ 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 0000000000..c2cbe1a308
--- /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 eac2f889af..755352aefd 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 1f7d2db17e..a1bdb00b11 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 0000000000..3e581ff684
--- /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_abs_logger_1D_space_event.comp b/mcstas-comps/union/Union_abs_logger_1D_space_event.comp
index 26769957c7..a099ad614c 100644
--- a/mcstas-comps/union/Union_abs_logger_1D_space_event.comp
+++ b/mcstas-comps/union/Union_abs_logger_1D_space_event.comp
@@ -497,7 +497,7 @@ INITIALIZE
strcat (this_abs_storage.Vars.option, option_string_part); // hardcoded min=0 for now
printf ("Monitor_nD option string: %s\n", this_abs_storage.Vars.option);
- Monitor_nD_Init (&this_abs_storage.DEFS, &this_abs_storage.Vars, 0.001, yheight, 0, 0, 0, 0, 0, 0, 0, 0); /* dims for mcdisplay */
+ Monitor_nD_Init (&this_abs_storage.DEFS, &this_abs_storage.Vars, 0.001, yheight, 0, 0, 0, 0, 0, 0, 0, 0, 0); /* dims for mcdisplay */
this_abs_storage.Vars.compcurpos = POS_A_CURRENT_COMP;
if (filename && strlen (filename) && strcmp (filename, "NULL") && strcmp (filename, "0"))
diff --git a/mcstas-comps/union/Union_abs_logger_event.comp b/mcstas-comps/union/Union_abs_logger_event.comp
index 3387a27bc6..e75ac3e63b 100644
--- a/mcstas-comps/union/Union_abs_logger_event.comp
+++ b/mcstas-comps/union/Union_abs_logger_event.comp
@@ -583,7 +583,7 @@ INITIALIZE
strcat (this_abs_storage.Vars.option, ", x y z vx vy vz t");
// strcat(this_abs_storage.Vars.option,", y t");
- Monitor_nD_Init (&this_abs_storage.DEFS, &this_abs_storage.Vars, 0.1, 0.1, 0, 0, 0, 0, 0, 0, 0, 0); /* dims for mcdisplay */
+ Monitor_nD_Init (&this_abs_storage.DEFS, &this_abs_storage.Vars, 0.1, 0.1, 0, 0, 0, 0, 0, 0, 0, 0, 0); /* dims for mcdisplay */
this_abs_storage.Vars.compcurpos = POS_A_CURRENT_COMP;
if (filename && strlen (filename) && strcmp (filename, "NULL") && strcmp (filename, "0"))
diff --git a/mcstas-comps/union/Union_abs_logger_nD.comp b/mcstas-comps/union/Union_abs_logger_nD.comp
index 9675d4bc3b..acdec0e39f 100644
--- a/mcstas-comps/union/Union_abs_logger_nD.comp
+++ b/mcstas-comps/union/Union_abs_logger_nD.comp
@@ -89,6 +89,7 @@
* restore_neutron: [0|1] Not functional for Union version
* geometry: [str] Name of an OFF file to specify a complex geometry detector
* nowritefile: [1] Not functional for Union version
+* nexus_bins: [1] NeXus mode only: store component BIN information
(-1 disable, 0 enable for list mode monitor, 1 enable for any montor)
*
* OUTPUT PARAMETERS:
*
@@ -108,7 +109,7 @@ SETTING PARAMETERS(string target_geometry="NULL",
xwidth=0, yheight=0, zdepth=0,
xmin=0, xmax=0, ymin=0, ymax=0, zmin=0, zmax=0,
int bins=0, min=-1e40, max=1e40, int restore_neutron=0, radius=0,
- string options="NULL", string filename="NULL",string geometry="NULL", int nowritefile=0,
+ string options="NULL", string filename="NULL",string geometry="NULL", int nowritefile=0, int nexus_bins=0,
string username1="NULL", string username2="NULL", string username3="NULL")
OUTPUT PARAMETERS ()
diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp
index 98dd5d6a5b..ae0f313459 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
diff --git a/mcstas-environment.yml b/mcstas-environment.yml
index 7989d20e43..191218cac3 100644
--- a/mcstas-environment.yml
+++ b/mcstas-environment.yml
@@ -24,9 +24,9 @@ dependencies:
- scipy
- pillow
- pyqtgraph
- - pyqt6
+ - pyqt
- qtpy
- #- qscintilla2
+ - qscintilla2
- nexusformat
- nexpy
- hdf5
diff --git a/mcxtrace-environment.yml b/mcxtrace-environment.yml
index e143a3fe59..e3cf61f939 100644
--- a/mcxtrace-environment.yml
+++ b/mcxtrace-environment.yml
@@ -24,9 +24,9 @@ dependencies:
- scipy
- pillow
- pyqtgraph
- - pyqt6
+ - pyqt
- qtpy
- #- qscintilla2
+ - qscintilla2
- nexusformat
- nexpy
- hdf5
diff --git a/meta-pkgs/windows/environment.yml b/meta-pkgs/windows/environment.yml
index 923a5b4953..651cbd1870 100644
--- a/meta-pkgs/windows/environment.yml
+++ b/meta-pkgs/windows/environment.yml
@@ -22,9 +22,9 @@ dependencies:
- scipy
- pillow
- pyqtgraph
- - pyqt6
+ - pyqt
- qtpy
- #- qscintilla2
+ - qscintilla2
- nexus
- nexusformat
- nexpy
diff --git a/meta-pkgs/windows/mcstas-environment.yml b/meta-pkgs/windows/mcstas-environment.yml
index 6f1c22e0eb..b478950345 100644
--- a/meta-pkgs/windows/mcstas-environment.yml
+++ b/meta-pkgs/windows/mcstas-environment.yml
@@ -21,9 +21,9 @@ dependencies:
- scipy
- pillow
- pyqtgraph
- - pyqt6
+ - pyqt
- qtpy
- #- qscintilla2
+ - qscintilla2
- nexus
- nexusformat
- nexpy
diff --git a/meta-pkgs/windows/mcxtrace-environment.yml b/meta-pkgs/windows/mcxtrace-environment.yml
index f6f234d2fb..94d7260699 100644
--- a/meta-pkgs/windows/mcxtrace-environment.yml
+++ b/meta-pkgs/windows/mcxtrace-environment.yml
@@ -21,9 +21,9 @@ dependencies:
- scipy
- pillow
- pyqtgraph
- - pyqt6
+ - pyqt
- qtpy
- #- qscintilla2
+ - qscintilla2
- nexus
- nexusformat
- nexpy
diff --git a/tools/Python/mccodelib/instrgeom.py b/tools/Python/mccodelib/instrgeom.py
index d6eb9eb1b7..910af79a25 100644
--- a/tools/Python/mccodelib/instrgeom.py
+++ b/tools/Python/mccodelib/instrgeom.py
@@ -578,7 +578,8 @@ def _get_points(self):
elif self.plane == 'yz':
return map(lambda p: cen.add(Vector3d(0, p.x, p.y)), square)
else:
- raise Exception('DrawCircle: invalid plane argument')
+ # Fallback to xy
+ return map(lambda p: cen.add(p), square)
def get_points_on_circle(self, steps=60):
''' returns points on the circle, transformed into the proper plane '''
diff --git a/tools/Python/mccodelib/mcplotloader.py b/tools/Python/mccodelib/mcplotloader.py
index 070ee63892..fe90a305fe 100644
--- a/tools/Python/mccodelib/mcplotloader.py
+++ b/tools/Python/mccodelib/mcplotloader.py
@@ -29,7 +29,12 @@ def __str__(self, *args, **kwargs):
class Data0D(DataMcCode):
# Bring back empty 'component' field, courtesy
# of the scan plotter (otherwise we fail abrubtly)
- DataMcCode.component=''
+ def __init__(self):
+ super(Data0D, self).__init__()
+ self.component=''
+ self.filename = ''
+ self.xlabel = ''
+ self.ylabel = ''
pass
@@ -329,13 +334,19 @@ def load(monfile):
if typ == 'array_0d':
print("load_monitor: Not loading 0d dataset %s" % monitorfile)
data = Data0D()
+ data.title='zero-dim monitor'
elif typ == 'array_1d':
data = _parse_1D_monitor(text)
elif typ == 'array_2d':
data = _parse_2D_monitor(text)
+ elif typ == 'list':
+ print('load_monitor: %s is an event list. loading suppressed' % f)
+ data = Data0D()
+ data.title='event list'
else:
print('load_monitor: unknown data format %s' % typ)
- data = None
+ data = Data0D()
+ data.title='Unknown format'
data.filepath = f
return data
else:
diff --git a/tools/Python/mccodelib/utils.py b/tools/Python/mccodelib/utils.py
index c2537c116c..0839a66d6b 100644
--- a/tools/Python/mccodelib/utils.py
+++ b/tools/Python/mccodelib/utils.py
@@ -918,6 +918,7 @@ def call_if_not_none(fct, *args):
''' shorthand utility for calling a function if it is defined, and otherwise ignoring it '''
if fct:
fct(*args)
+
if not cwd:
cwd = os.getcwd()
diff --git a/tools/Python/mcdisplay/webgl/package.json b/tools/Python/mcdisplay/webgl/package.json
index 558f43ce81..5a99abfdf5 100644
--- a/tools/Python/mcdisplay/webgl/package.json
+++ b/tools/Python/mcdisplay/webgl/package.json
@@ -25,7 +25,7 @@
"three": "^0.165.0",
"three-lut": "^79.0.2",
"three-orbitcontrols": "^2.110.3",
- "vite": "^5.2.13"
+ "vite": "^6.4.2"
},
"devDependencies": {
"typescript": "^5.4.5",
diff --git a/tools/Python/mcgui/mcgui.py b/tools/Python/mcgui/mcgui.py
index 0bd2816188..2591983f10 100755
--- a/tools/Python/mcgui/mcgui.py
+++ b/tools/Python/mcgui/mcgui.py
@@ -459,6 +459,23 @@ def run(self, fixed_params, params, inspect=None):
self.__emitter.status('Started simulation/trace in background shell...')
subprocess.Popen(runstr, shell=True)
+ def vis_run(self, tool=''):
+
+ runstr = tool + ' ' + os.path.basename(self.__instrFile) + ' -n100 -y'
+
+ # Add & for backgrounding on Unix systems
+ if not os.name == 'nt':
+ runstr = runstr + ' &'
+ else:
+ runstr = 'start ' + runstr
+
+ # Ensure assembled runstr is a string, not a QString
+ runstr = str(runstr)
+
+ self.__emitter.status('Starting visualisation in background shell...')
+ self.__emitter.message(runstr, gui=False)
+ subprocess.Popen(runstr, shell=True)
+
def __runFinished(self, process_returncode):
self.__fireSimStateUpdate()
if process_returncode == 0:
@@ -741,35 +758,22 @@ def handleMcDisplayWeb(self):
DISPLAY="mcdisplay"
else:
DISPLAY="mxdisplay"
- self.emitter.status('Running ' + DISPLAY + '-webgl...')
- try:
- cmd = DISPLAY+'-webgl --default -n100 ' + os.path.basename(self.state.getInstrumentFile()) + '&'
- self.emitter.message(cmd, gui=True)
- self.emitter.message('', gui=True)
-
- def messg(s): self.emitter.message(s)
- def messg_err(s): self.emitter.message(s, err_msg=True)
- utils.run_subtool_to_completion(cmd, stdout_cb=messg, stderr_cb=messg_err)
- finally:
- self.emitter.status('')
-
+ tool=DISPLAY + "-webgl-classic"
+
+ self.emitter.status('Running ' + tool)
+ self.state.vis_run(tool=tool)
+
def handleMcDisplay2D(self):
if mccode_config.configuration["MCCODE"]=="mcstas":
DISPLAY="mcdisplay"
else:
DISPLAY="mxdisplay"
- self.emitter.status('Running ' + DISPLAY + '-pyqtgraph...')
- try:
- cmd = DISPLAY+'-pyqtgraph --default -n100 ' + os.path.basename(self.state.getInstrumentFile()) + '&'
- self.emitter.message(cmd, gui=True)
- self.emitter.message('', gui=True)
-
- def messg(s): self.emitter.message(s)
- def messg_err(s): self.emitter.message(s, err_msg=True)
- utils.run_subtool_to_completion(cmd, stdout_cb=messg, stderr_cb=messg_err)
- finally:
- self.emitter.status('')
-
+ tool=DISPLAY + "-pyqtgraph"
+
+ self.emitter.status('Running ' + tool)
+ self.state.vis_run(tool=tool)
+
+
def handleHelpWeb(self):
# open the mcstas homepage
mcurl = 'http://www.'+mccode_config.configuration["MCCODE"]+'.org'
diff --git a/tools/Python/mcgui/mw.ui b/tools/Python/mcgui/mw.ui
index 8318961c45..aeaee83037 100644
--- a/tools/Python/mcgui/mw.ui
+++ b/tools/Python/mcgui/mw.ui
@@ -349,7 +349,7 @@
- Display-3D
+ Display-3D (classic)
diff --git a/tools/Python/mcgui/viewclasses.py b/tools/Python/mcgui/viewclasses.py
index 845d1178b4..52d97ee647 100644
--- a/tools/Python/mcgui/viewclasses.py
+++ b/tools/Python/mcgui/viewclasses.py
@@ -1043,19 +1043,22 @@ def _wEventFilter(subject, object, event):
# handle focus off
elif event.type() == QtCore.QEvent.Type.FocusOut:
- if edt.text() == '':
+ if edt.text() == '': # Empty
font = QtGui.QFont()
font.setItalic(True)
edt.setFont(font)
edt.setStyleSheet("color: grey;")
edt.setText(edt.defval)
- elif edt.text() == edt.defval:
+ elif edt.text() == edt.defval: # Default value
edt.setText(edt.defval)
font = QtGui.QFont()
font.setItalic(True)
edt.setFont(font)
edt.setStyleSheet("color: grey;")
-
+ else: # Something else
+ font = QtGui.QFont()
+ edt.setFont(font)
+ edt.setStyleSheet("color: black;")
return False
# init
diff --git a/tools/Python/mcgui/widgets.py b/tools/Python/mcgui/widgets.py
index f5dd00179e..925397ac46 100644
--- a/tools/Python/mcgui/widgets.py
+++ b/tools/Python/mcgui/widgets.py
@@ -263,7 +263,7 @@ def retranslateUi(self, MainWindow):
self.actionConfiguration.setToolTip(_translate("MainWindow", "mccode configuration"))
self.actionCompile_Instrument_MPI.setText(_translate("MainWindow", "Compile Instrument (MPI)"))
self.actionMcdoc.setText(_translate("MainWindow", "mcdoc Component Reference"))
- self.actionDisplay.setText(_translate("MainWindow", "Display-3D"))
+ self.actionDisplay.setText(_translate("MainWindow", "Display-3D (classic)"))
self.actionDisplay_2d.setText(_translate("MainWindow", "Display-2D"))
self.actionPlotOther.setText(_translate("MainWindow", "Plot Other Results"))
self.actionPlotOther.setShortcut(_translate("MainWindow", "Ctrl+Shift+P"))
diff --git a/tools/Python/mcplot/pyqtgraph/plotfuncs.py b/tools/Python/mcplot/pyqtgraph/plotfuncs.py
index 10afc0b6c6..13ce048b32 100644
--- a/tools/Python/mcplot/pyqtgraph/plotfuncs.py
+++ b/tools/Python/mcplot/pyqtgraph/plotfuncs.py
@@ -74,7 +74,7 @@ def paint(self, p, *args):
def plot_Data0D(data, plt, log=False, legend=True, icolormap=0, verbose=True, fontsize=10):
''' creates an empty plot '''
- plt.setTitle("zero-dim monitor")
+ plt.setTitle(data.title)
return plt.getViewBox()