From 1958da6514666b6cb2189cba4bd57eed13a342f7 Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 14:00:45 +0200 Subject: [PATCH 01/36] Simplify Volume is only an absorber but not vacuum path in union master --- mcstas-comps/union/Union_master.comp | 55 +++++++++++++++++++--------- 1 file changed, 38 insertions(+), 17 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index bb1b007b1a..b6f1abcf18 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -83,6 +83,40 @@ SHARE #ifndef MASTER_DETECTOR #define MASTER_DETECTOR dummy #endif + + int volume_is_only_absorber(struct Volume_struct* Volume){ + // This function returns true if a volume does not have any physical processes + // and if the volume is not a vacuum. + if (!Volume->p_physics->number_of_processes && !Volume->p_physics->is_vacuum) return 1; + return 0; + } + void adjust_abs_weight_factor(struct Volume_struct* Volume, + double* my_sum_plus_abs, + double* length_to_boundary, + double* v_length, + double* time_to_boundary, + double* abs_weight_factor, + int* abs_weight_factor_set){ + *my_sum_plus_abs = Volume->p_physics->my_a * (2200 / *v_length); + *length_to_boundary = *time_to_boundary * *v_length; + + *abs_weight_factor = exp (-Volume->p_physics->my_a * 2200 * *time_to_boundary); + *abs_weight_factor_set = 1; + + #ifdef Union_trace_verbal_setting + printf ("name of material: %s \n", Volumes->name); + printf ("length to boundery = %f\n", length_to_boundary); + printf ("absorption cross section = %f\n", Volumes->p_physics->my_a); + printf ("chance to get through this length of absorber: %f %%\n", + 100 * exp (-Volumes->p_physics->my_a * length_to_boundary)); + #endif + + } + + + + + %} DECLARE @@ -1544,23 +1578,10 @@ TRACE // Check if a scattering event should occur if (current_volume != 0) { // Volume 0 is always vacuum, and if this is the current volume, an event will not occur - if (Volumes[current_volume]->p_physics->number_of_processes == 0) { // If there are no processes, the volume could be vacuum or an absorber - if (Volumes[current_volume]->p_physics->is_vacuum == 0) { - // This volume does not have physical processes but does have an absorption cross section, so the ray weight is reduced accordingly - - my_sum_plus_abs = Volumes[current_volume]->p_physics->my_a * (2200 / v_length); - length_to_boundary = time_to_boundery * v_length; - - abs_weight_factor = exp (-Volumes[current_volume]->p_physics->my_a * 2200 * time_to_boundery); - abs_weight_factor_set = 1; - - #ifdef Union_trace_verbal_setting - printf ("name of material: %s \n", Volumes[current_volume]->name); - printf ("length to boundery = %f\n", length_to_boundary); - printf ("absorption cross section = %f\n", Volumes[current_volume]->p_physics->my_a); - printf ("chance to get through this length of absorber: %f %%\n", 100 * exp (-Volumes[current_volume]->p_physics->my_a * length_to_boundary)); - #endif - } + if (volume_is_only_absorber(Volumes[current_volume])) { // If there are no processes, the volume could be vacuum or an absorber + adjust_abs_weight_factor(Volumes[current_volume], &my_sum_plus_abs, + &length_to_boundary, &v_length, &time_to_boundery, + &abs_weight_factor, &abs_weight_factor_set); } 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; From 59fac42b75c3004668cd8c9fc37f3207a6bdc2be Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 14:33:06 +0200 Subject: [PATCH 02/36] Gather focus preprocessing into its own function --- mcstas-comps/union/Union_master.comp | 98 ++++++++++++++++++++-------- 1 file changed, 71 insertions(+), 27 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index b6f1abcf18..460e7e88c7 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -113,6 +113,43 @@ SHARE } + void choose_scattering_point_and_ray_direction(double* forced_length_to_scattering, + double* safety_distance, double* safety_distance2, + double* length_to_boundary, _class_particle *_particle, + Coords *ray_velocity, Coords* ray_position_geometry, + Coords *ray_position, + struct Volume_struct* Volume, struct focus_data_struct* this_focus_data){ + // Sample length_to_scattering in linear manner + *forced_length_to_scattering = *safety_distance + rand01 () * (*length_to_boundary - *safety_distance2); + + *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); // Test for root cause + // 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 scattering_displacement = coords_scalar_mult (direction, *forced_length_to_scattering); + Coords forced_ray_scattering_point = coords_add (*ray_position, scattering_displacement); + *ray_position_geometry + = coords_sub (forced_ray_scattering_point, Volume->geometry.center); // ray_position relative to geometry center + + // Calculate the aim for non isotropic processes + this_focus_data = &Volume->geometry.focus_data_array.elements[0]; + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, *ray_position_geometry); // Aim vector for this ray + + #ifdef Union_trace_verbal_setting + printf ("Prepared for focus in cross section calculation in volume: %s \n", Volume->name); + printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); + print_position (*ray_position, "ray_position"); + print_position (direction, "direction"); + print_position (scattering_displacement, "scattering_displacement"); + print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); + print_position (*ray_position_geometry, "ray_position_geometry"); + printf ("for isotropic processes this RayAim is used \n"); + print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + #endif + + + } + + @@ -1586,6 +1623,8 @@ TRACE // 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; + double forced_length_to_scattering; + Coords ray_position_geometry; my_sum = 0; k[0] = V2K * vx; k[1] = V2K * vy; @@ -1594,37 +1633,42 @@ TRACE wavevector = coords_set (k[0], k[1], k[2]); length_to_boundary = time_to_boundery * v_length; - double forced_length_to_scattering; - Coords ray_position_geometry; + // If any process in this material needs focusing, sample scattering position and update focus_data accordingly if (current_p_physics->any_process_needs_cross_section_focus == 1) { + choose_scattering_point_and_ray_direction(&forced_length_to_scattering, + &safety_distance, &safety_distance2, + &length_to_boundary, _particle, + &ray_velocity, &ray_position_geometry, + &ray_position, + Volumes[current_volume], this_focus_data); // Sample length_to_scattering in linear manner - forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); - - ray_velocity = coords_set (vx, vy, vz); // Test for root cause - // 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 scattering_displacement = coords_scalar_mult (direction, forced_length_to_scattering); - Coords forced_ray_scattering_point = coords_add (ray_position, scattering_displacement); - ray_position_geometry - = coords_sub (forced_ray_scattering_point, Volumes[current_volume]->geometry.center); // ray_position relative to geometry center - - // Calculate the aim for non isotropic processes - this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry); // Aim vector for this ray - - #ifdef Union_trace_verbal_setting - printf ("Prepared for focus in cross section calculation in volume: %s \n", Volumes[current_volume]->name); - printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); - print_position (ray_position, "ray_position"); - print_position (direction, "direction"); - print_position (scattering_displacement, "scattering_displacement"); - print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); - print_position (ray_position_geometry, "ray_position_geometry"); - printf ("for isotropic processes this RayAim is used \n"); - print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - #endif + // forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); + // + // ray_velocity = coords_set (vx, vy, vz); // Test for root cause + // // 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 scattering_displacement = coords_scalar_mult (direction, forced_length_to_scattering); + // Coords forced_ray_scattering_point = coords_add (ray_position, scattering_displacement); + // ray_position_geometry + // = coords_sub (forced_ray_scattering_point, Volumes[current_volume]->geometry.center); // ray_position relative to geometry center + // + // // Calculate the aim for non isotropic processes + // this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; + // this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry); // Aim vector for this ray + // + // #ifdef Union_trace_verbal_setting + // printf ("Prepared for focus in cross section calculation in volume: %s \n", Volumes[current_volume]->name); + // printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); + // print_position (ray_position, "ray_position"); + // print_position (direction, "direction"); + // print_position (scattering_displacement, "scattering_displacement"); + // print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); + // print_position (ray_position_geometry, "ray_position_geometry"); + // printf ("for isotropic processes this RayAim is used \n"); + // print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + // #endif /* // update focus data for this ray (could limit this to only update the necessary focus_data element, but there are typically very few) From be52fe54aef8de8513b93b61d7cd5c79599b5a87 Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 16:00:36 +0200 Subject: [PATCH 03/36] transfer non isotropic wavevector rotation into its own function --- mcstas-comps/union/Union_master.comp | 104 +++++++++++---------------- 1 file changed, 40 insertions(+), 64 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 460e7e88c7..3d20dfc5f5 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -148,7 +148,41 @@ SHARE } - + + + void transform_wavevector_into_local_coord_system(struct Volume_struct* Volume, + Coords * wavevector_rotated, double (*k_rotated)[3], + int *p_index, Coords* wavevector, Coords* ray_position_geometry, + struct focus_data_struct* this_focus_data) { + + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; + *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); + coords_get (*wavevector_rotated, k_rotated[0], k_rotated[1], k_rotated[2]); + + if (Volume->p_physics->p_scattering_array[*p_index].needs_cross_section_focus == 1) { + // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time + // focusing + Coords ray_position_geometry_rotated + = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *ray_position_geometry); + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray + #ifdef Union_trace_verbal_setting + printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", *p_index); + print_position (*ray_position_geometry, "ray_position_geometry"); + print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); + print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + #endif + } + } + + + int process_needs_multiple_sampling(struct physics_struct* current_p_physics, + struct scattering_process_struct* process){ + if (current_p_physics->sampling_points !=0 && process->needs_cross_section_focus){ + if (process->sampling_points != -1) + return 1; + } + return 0; + } @@ -1643,50 +1677,6 @@ TRACE &ray_velocity, &ray_position_geometry, &ray_position, Volumes[current_volume], this_focus_data); - // Sample length_to_scattering in linear manner - // forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); - // - // ray_velocity = coords_set (vx, vy, vz); // Test for root cause - // // 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 scattering_displacement = coords_scalar_mult (direction, forced_length_to_scattering); - // Coords forced_ray_scattering_point = coords_add (ray_position, scattering_displacement); - // ray_position_geometry - // = coords_sub (forced_ray_scattering_point, Volumes[current_volume]->geometry.center); // ray_position relative to geometry center - // - // // Calculate the aim for non isotropic processes - // this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; - // this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry); // Aim vector for this ray - // - // #ifdef Union_trace_verbal_setting - // printf ("Prepared for focus in cross section calculation in volume: %s \n", Volumes[current_volume]->name); - // printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); - // print_position (ray_position, "ray_position"); - // print_position (direction, "direction"); - // print_position (scattering_displacement, "scattering_displacement"); - // print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); - // print_position (ray_position_geometry, "ray_position_geometry"); - // printf ("for isotropic processes this RayAim is used \n"); - // print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - // #endif - - /* - // 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 - } - - printf("calculated forced_length_to_scattering = %lf, new RayAim \n", forced_length_to_scattering); - print_position(direction, "direction"); - print_position(scattering_displacement, "scattering_displacement"); - print_position(forced_ray_scattering_point, "forced_ray_scattering_point"); - print_position(ray_position_geometry, "ray_position_geometry"); - print_position(this_focus_data->RayAim, "this_focus_data->RayAim"); - */ - } else { forced_length_to_scattering = -1.0; // Signals that no forcing needed, could also if on the selected process struct } @@ -1700,24 +1690,10 @@ TRACE if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { // If the process is not isotropic, the wavevector is transformed into the local coordinate system of the process - int non_isotropic_rot_index = Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; - wavevector_rotated = rot_apply (Volumes[current_volume]->geometry.process_rot_matrix_array[non_isotropic_rot_index], wavevector); - coords_get (wavevector_rotated, &k_rotated[0], &k_rotated[1], &k_rotated[2]); - - if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].needs_cross_section_focus == 1) { - // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time - // focusing - Coords ray_position_geometry_rotated - = rot_apply (Volumes[current_volume]->geometry.process_rot_matrix_array[non_isotropic_rot_index], ray_position_geometry); - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray - #ifdef Union_trace_verbal_setting - printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", p_index); - print_position (ray_position_geometry, "ray_position_geometry"); - print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); - print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - #endif - } - + transform_wavevector_into_local_coord_system(Volumes[current_volume], + &wavevector_rotated, &k_rotated, + &p_index, &wavevector, &ray_position_geometry, + this_focus_data); } else { k_rotated[0] = k[0]; k_rotated[1] = k[1]; @@ -1732,7 +1708,7 @@ TRACE 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)) { + if (process_needs_multiple_sampling(current_p_physics, process)) { // Populate length and probability arrays: Coords original_position = coords_set (x, y, z); for (int i = 0; i < current_p_physics->sampling_points; i++) { From 8d0c5574682724661ba2c61367685d38ef9207d8 Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 16:53:11 +0200 Subject: [PATCH 04/36] Make inhomogenous sampling loop more transparent by moving most logic into a functions --- mcstas-comps/union/Union_master.comp | 81 +++++++++++++++++----------- 1 file changed, 51 insertions(+), 30 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 3d20dfc5f5..1bee050a35 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -152,8 +152,7 @@ SHARE void transform_wavevector_into_local_coord_system(struct Volume_struct* Volume, Coords * wavevector_rotated, double (*k_rotated)[3], - int *p_index, Coords* wavevector, Coords* ray_position_geometry, - struct focus_data_struct* this_focus_data) { + int *p_index, Coords* wavevector, Coords* ray_position_geometry) { int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); @@ -164,7 +163,11 @@ SHARE // focusing Coords ray_position_geometry_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *ray_position_geometry); + + int focus_data_index = Volume->geometry.focus_array_indices.elements[*p_index]; + struct focus_data_struct* this_focus_data = &Volume->geometry.focus_data_array.elements[focus_data_index]; this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray + #ifdef Union_trace_verbal_setting printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", *p_index); print_position (*ray_position_geometry, "ray_position_geometry"); @@ -175,7 +178,7 @@ SHARE } - int process_needs_multiple_sampling(struct physics_struct* current_p_physics, + int process_needs_inhomogenous_sampling(struct physics_struct* current_p_physics, struct scattering_process_struct* process){ if (current_p_physics->sampling_points !=0 && process->needs_cross_section_focus){ if (process->sampling_points != -1) @@ -185,6 +188,38 @@ SHARE } + void set_cumul_dist_array(struct physics_struct* current_p_physics, int 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; + } + + void move_and_aim_neutron(struct physics_struct* current_p_physics, int i, + struct scattering_process_struct* process, + _class_particle* _particle, + struct Volume_struct* Volume, int* p_index, + Coords* ray_velocity, Coords* ray_position, + struct focus_data_struct* this_focus_data + ){ + // Transport neutron to place inside geometry + *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->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, Volume->geometry.center); + // Also focus the ray at this point, if the component needs focusing + if (process->needs_cross_section_focus) { + if (Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index != -1) { + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; + sampling_point_geometry + = rot_apply (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, &_particle->x, &_particle->y, &_particle->z); + } + + @@ -1692,8 +1727,7 @@ TRACE // If the process is not isotropic, the wavevector is transformed into the local coordinate system of the process transform_wavevector_into_local_coord_system(Volumes[current_volume], &wavevector_rotated, &k_rotated, - &p_index, &wavevector, &ray_position_geometry, - this_focus_data); + &p_index, &wavevector, &ray_position_geometry); } else { k_rotated[0] = k[0]; k_rotated[1] = k[1]; @@ -1708,38 +1742,25 @@ TRACE double mu; current_p_physics->dist = length_to_boundary / current_p_physics->sampling_points - safety_distance; - if (process_needs_multiple_sampling(current_p_physics, process)) { - // Populate length and probability arrays: + if (process_needs_inhomogenous_sampling(current_p_physics, process)) { + // If mulitple sampling is required for this process, populate the neutron distance and probability arrays: Coords original_position = coords_set (x, y, z); + double mu; + *p_my_trace = 0; + this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; 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); - } + set_cumul_dist_array(current_p_physics, i); + move_and_aim_neutron(current_p_physics, i, process, _particle, + Volumes[current_volume], &p_index, + &ray_velocity, &ray_position, + this_focus_data); // 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; + *p_my_trace += mu / current_p_physics->sampling_points; } - 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) { From e45e91bec31451e4a5c799fb6636b63c7c922aac Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 17:36:05 +0200 Subject: [PATCH 05/36] Move inhomogenous distance sampling to a function --- mcstas-comps/union/Union_master.comp | 96 ++++++++++++++++------------ 1 file changed, 56 insertions(+), 40 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 1bee050a35..a17c0481e6 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -220,7 +220,55 @@ SHARE } + int mu_and_intersect_dist_safeguard(double mu_sum, double length_to_boundary, double safety_distance2, + int *scattering_event){ + if (mu_sum < 1E-18){ + scattering_event = 0; + return 0; + } + if (length_to_boundary < safety_distance2){ + scattering_event = 0; + return 0; + } + return 1; + } + void sample_inhomogenous_transmission_probability(struct physics_struct* current_p_physics, + struct Volume_struct* Volume, double* real_transmission_probability, + double v_length + ){ + // 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 < Volume->p_physics->number_of_processes; i++) { + struct scattering_process_struct* process_i = &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 = 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]; + } %} @@ -1701,9 +1749,6 @@ TRACE p_my_trace = my_trace; wavevector = coords_set (k[0], k[1], k[2]); length_to_boundary = time_to_boundery * v_length; - - - // If any process in this material needs focusing, sample scattering position and update focus_data accordingly if (current_p_physics->any_process_needs_cross_section_focus == 1) { choose_scattering_point_and_ray_direction(&forced_length_to_scattering, @@ -1761,7 +1806,8 @@ TRACE } coords_get (original_position, &x, &y, &z); - } else { + } + if (!process_needs_inhomogenous_sampling(current_p_physics, process)){ 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; @@ -1788,45 +1834,15 @@ TRACE // New flow:length_to_boundary // Calculate if scattering happens based on my_sub_plus_abs - if (my_sum < 1E-18) { - scattering_event = 0; - } else if (length_to_boundary < safety_distance2) { - scattering_event = 0; - } else { + if (mu_and_intersect_dist_safeguard(my_sum, length_to_boundary, safety_distance2, + &scattering_event)){ 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]; + sample_inhomogenous_transmission_probability(current_p_physics, + Volumes[current_volume], + &real_transmission_probability, + v_length); } // printf("Trans prop = %g\n", real_transmission_probability); if (Volumes[current_volume]->geometry.geometry_p_interact != 0) { From e322f23d4c97ca826b104dd8ff06508fefdaa1f4 Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 18:17:02 +0200 Subject: [PATCH 06/36] Move p_interact behaviour and inhomogenous scattering point sampling to their own functions --- mcstas-comps/union/Union_master.comp | 107 ++++++++++++++++++--------- 1 file changed, 72 insertions(+), 35 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index a17c0481e6..ae89f5c012 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -270,6 +270,61 @@ SHARE *real_transmission_probability = current_p_physics->cumul_transmission_prob[current_p_physics->sampling_points - 1]; } + double sample_scattering_point_inhomogenous(struct Volume_struct* Volume, + struct physics_struct* current_p_physics, + double* abs_weight_factor, + double* v_length, + double* safety_distance, + int* selected_sampling){ + + // 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 = 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; + return current_p_physics->cumul_dists[*selected_sampling] - current_p_physics->dist / 2 + sampled_dist; + } + + + int p_interact_is_set(struct Volume_struct* Volume){ + if (Volume->geometry.geometry_p_interact != 0){ + return 1; + } + return 0; + } + + void p_interact_check_scattering_event(struct Volume_struct* Volume, + int* scattering_event, + double* weight, + double real_transmission_prob){ + double mc_transmission_prob = 1 - Volume->geometry.geometry_p_interact; + *scattering_event = rand01() > mc_transmission_prob; + if (*scattering_event){ + // Scattering event happens, this is the correction for the weight + *weight *= (1.0 - real_transmission_prob)/(1.0 - mc_transmission_prob); + } + else { + // Scattering event does not happen, this is the appropriate correction + *weight *= real_transmission_prob / mc_transmission_prob; + } + } + %} @@ -1836,27 +1891,25 @@ TRACE // Calculate if scattering happens based on my_sub_plus_abs if (mu_and_intersect_dist_safeguard(my_sum, length_to_boundary, safety_distance2, &scattering_event)){ + // First calculate the transmission probability 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) { + } + + if (current_p_physics->sampling_points != 0) { sample_inhomogenous_transmission_probability(current_p_physics, Volumes[current_volume], &real_transmission_probability, v_length); } - // 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))) { - // Scattering event happens, this is the correction for the weight - p *= (1.0 - real_transmission_probability) / (1.0 - mc_transmission_probability); - } else { - // Scattering event does not happen, this is the appropriate correction - p *= real_transmission_probability / mc_transmission_probability; - } + + // Then check if we scatter. + if (p_interact_is_set(Volumes[current_volume])) { + p_interact_check_scattering_event(Volumes[current_volume], + &scattering_event, + &p, real_transmission_probability); } else { // probability to scatter is the natural value - // printf("Real transmission prob %g\n", real_transmission_probability); scattering_event = rand01 () > real_transmission_probability; } } @@ -1873,29 +1926,13 @@ TRACE 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; - } + sample_scattering_point_inhomogenous(Volumes[current_volume], + current_p_physics, + &abs_weight_factor, + &v_length, + &safety_distance, + &selected_sampling); + } // 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 From 3128fbb0300bf0dbf72cb0c81ca201830087e919 Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 19:03:25 +0200 Subject: [PATCH 07/36] Move process choice after scattering into functions for p_interact and inhomogenous --- mcstas-comps/union/Union_master.comp | 108 ++++++++++++++++++--------- 1 file changed, 73 insertions(+), 35 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index ae89f5c012..d15ef32306 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -326,6 +326,59 @@ SHARE } + void p_interact_select_process(struct Volume_struct* Volume, + double* my_trace_fraction_control, + double* my_trace, + double* total_process_interact, + double* culmative_probability, + double* mc_prop, + double* weight, + double* my_sum, + int* selected_process + ){ + // Interact_fraction is used to influence the choice of process in this material + *mc_prop = rand01 (); + *culmative_probability = 0; + *total_process_interact = 1.0; + + // If any of the processes have probability 0, they are excluded from the selection + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + if (my_trace[i] < 1E-18) { + // When this happens, the total force probability is corrected and the probability for this particular instance is set to 0 + *total_process_interact -= Volume->p_physics->p_scattering_array[i].process_p_interact; + my_trace_fraction_control[i] = 0; + // In cases where my_trace is not zero, the forced fraction is still used. + } else + my_trace_fraction_control[i] = Volume->p_physics->p_scattering_array[i].process_p_interact; + } + // Randomly select a process using the weights stored in my_trace_fraction_control divided by total_process_interact + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + *culmative_probability += my_trace_fraction_control[i] / *total_process_interact; + if (*culmative_probability > *mc_prop) { + *selected_process = i; + *weight *= (my_trace[i] / *my_sum) * (*total_process_interact / my_trace_fraction_control[i]); + break; + } + } + } + + void inhomogenous_choose_process(struct physics_struct* current_p_physics, + struct Volume_struct* Volume, + double* culmative_probability, + double* mc_prop, + double* my_sum, + int* selected_sampling, + int* selected_process){ + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + *culmative_probability += current_p_physics->mus[i][*selected_sampling] / *my_sum; + if (*culmative_probability > *mc_prop) { + *selected_process = i; + break; + } + } + } + + %} DECLARE @@ -1932,39 +1985,13 @@ TRACE &v_length, &safety_distance, &selected_sampling); - } + } // 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 selected_process = 0; } else { - if (Volumes[current_volume]->p_physics->interact_control == 1) { - // Interact_fraction is used to influence the choice of process in this material - mc_prop = rand01 (); - culmative_probability = 0; - total_process_interact = 1.0; - - // If any of the processes have probability 0, they are excluded from the selection - for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { - if (my_trace[iterator] < 1E-18) { - // When this happens, the total force probability is corrected and the probability for this particular instance is set to 0 - total_process_interact -= Volumes[current_volume]->p_physics->p_scattering_array[iterator].process_p_interact; - my_trace_fraction_control[iterator] = 0; - // In cases where my_trace is not zero, the forced fraction is still used. - } else - my_trace_fraction_control[iterator] = Volumes[current_volume]->p_physics->p_scattering_array[iterator].process_p_interact; - } - // Randomly select a process using the weights stored in my_trace_fraction_control divided by total_process_interact - for (iterator = 0; iterator < Volumes[current_volume]->p_physics->number_of_processes; iterator++) { - culmative_probability += my_trace_fraction_control[iterator] / total_process_interact; - if (culmative_probability > mc_prop) { - selected_process = iterator; - p *= (my_trace[iterator] / my_sum) * (total_process_interact / my_trace_fraction_control[iterator]); - break; - } - } - - } else { + if (Volumes[current_volume]->p_physics->interact_control != 1) { // Select a process based on their relative attenuations factors mc_prop = rand01 (); culmative_probability = 0; @@ -1977,14 +2004,25 @@ TRACE } } } 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; - } - } + inhomogenous_choose_process(current_p_physics, + Volumes[current_volume], + &culmative_probability, + &mc_prop, + &my_sum, + &selected_sampling, + &selected_process); + } + } else { + p_interact_select_process(Volumes[current_volume], + my_trace_fraction_control, + my_trace, + &total_process_interact, + &culmative_probability, + &mc_prop, + &p, + &my_sum, + &selected_process); } } From 78189d5b2698eac8741c198612fa9c22ed7543ad Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 19:12:24 +0200 Subject: [PATCH 08/36] Move using forced length to scatter from to a separate function --- mcstas-comps/union/Union_master.comp | 32 +++++++++++++++++++--------- 1 file changed, 22 insertions(+), 10 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index d15ef32306..d3b5a9dea1 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -378,6 +378,25 @@ SHARE } } + void dir_focus_set_scat_length(double* length_to_scattering, + double* forced_length_to_scattering, + double* weight, + double* length_to_boundary, + double* my_sum_plus_abs){ + // 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 + + } + %} @@ -2030,16 +2049,9 @@ TRACE 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 - + dir_focus_set_scat_length(&length_to_scattering, &forced_length_to_scattering, + &p, &length_to_boundary, + &my_sum_plus_abs); } else { // Decided the ray scatters, choose where on truncated exponential from safety_distance to length_to_boundary - safety_distance length_to_scattering From 6550419b85fb25933261bf2150ec3f3d949904c6 Mon Sep 17 00:00:00 2001 From: Diablo Date: Tue, 14 Jul 2026 19:14:50 +0200 Subject: [PATCH 09/36] fix p being used instead of weight in dir length function --- mcstas-comps/union/Union_master.comp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index d3b5a9dea1..020915ec08 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -388,7 +388,7 @@ SHARE // 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 + *weight *= *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 From 4265eef4e0f35ca01da8c4c49a6be85337d98d8e Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 15 Jul 2026 12:44:56 +0200 Subject: [PATCH 10/36] Fix pointer magic that resulted in only writing to the first value of k_rotated in transform_wavevector function --- mcstas-comps/union/Union_master.comp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 020915ec08..f4152f5103 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -156,7 +156,7 @@ SHARE int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); - coords_get (*wavevector_rotated, k_rotated[0], k_rotated[1], k_rotated[2]); + coords_get (*wavevector_rotated, &(*k_rotated)[0], &(*k_rotated)[1], &(*k_rotated)[2]); if (Volume->p_physics->p_scattering_array[*p_index].needs_cross_section_focus == 1) { // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time From 17a1be5b32bbed9f1615cbfd6ae162ed09512913 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 15 Jul 2026 13:50:32 +0200 Subject: [PATCH 11/36] Assign scattering length for inhomogenous process --- mcstas-comps/union/Union_master.comp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index f4152f5103..1b7b9599b8 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -1998,7 +1998,7 @@ TRACE printf ("WARNING: Absorption weight factor above 1! Should not happen! \n"); // Select distance to scattering position if (current_p_physics->sampling_points != 0) { - sample_scattering_point_inhomogenous(Volumes[current_volume], + length_to_scattering = sample_scattering_point_inhomogenous(Volumes[current_volume], current_p_physics, &abs_weight_factor, &v_length, From 4b687fcb00b6a9a6d16e86e1fc8c2c554f5879bd Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 15 Jul 2026 13:54:23 +0200 Subject: [PATCH 12/36] Apply mccode-clangformat to Union_master.comp --- mcstas-comps/union/Union_master.comp | 340 +++++++++++---------------- 1 file changed, 131 insertions(+), 209 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 1b7b9599b8..b79b7e80b1 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -84,41 +84,35 @@ SHARE #define MASTER_DETECTOR dummy #endif - int volume_is_only_absorber(struct Volume_struct* Volume){ + int + volume_is_only_absorber (struct Volume_struct* Volume) { // This function returns true if a volume does not have any physical processes // and if the volume is not a vacuum. - if (!Volume->p_physics->number_of_processes && !Volume->p_physics->is_vacuum) return 1; + if (!Volume->p_physics->number_of_processes && !Volume->p_physics->is_vacuum) + return 1; return 0; } - void adjust_abs_weight_factor(struct Volume_struct* Volume, - double* my_sum_plus_abs, - double* length_to_boundary, - double* v_length, - double* time_to_boundary, - double* abs_weight_factor, - int* abs_weight_factor_set){ - *my_sum_plus_abs = Volume->p_physics->my_a * (2200 / *v_length); - *length_to_boundary = *time_to_boundary * *v_length; - - *abs_weight_factor = exp (-Volume->p_physics->my_a * 2200 * *time_to_boundary); - *abs_weight_factor_set = 1; - - #ifdef Union_trace_verbal_setting - printf ("name of material: %s \n", Volumes->name); - printf ("length to boundery = %f\n", length_to_boundary); - printf ("absorption cross section = %f\n", Volumes->p_physics->my_a); - printf ("chance to get through this length of absorber: %f %%\n", - 100 * exp (-Volumes->p_physics->my_a * length_to_boundary)); - #endif + void + adjust_abs_weight_factor (struct Volume_struct* Volume, double* my_sum_plus_abs, double* length_to_boundary, double* v_length, double* time_to_boundary, + double* abs_weight_factor, int* abs_weight_factor_set) { + *my_sum_plus_abs = Volume->p_physics->my_a * (2200 / *v_length); + *length_to_boundary = *time_to_boundary * *v_length; + + *abs_weight_factor = exp (-Volume->p_physics->my_a * 2200 * *time_to_boundary); + *abs_weight_factor_set = 1; + #ifdef Union_trace_verbal_setting + printf ("name of material: %s \n", Volumes->name); + printf ("length to boundery = %f\n", length_to_boundary); + printf ("absorption cross section = %f\n", Volumes->p_physics->my_a); + printf ("chance to get through this length of absorber: %f %%\n", 100 * exp (-Volumes->p_physics->my_a * length_to_boundary)); + #endif } - void choose_scattering_point_and_ray_direction(double* forced_length_to_scattering, - double* safety_distance, double* safety_distance2, - double* length_to_boundary, _class_particle *_particle, - Coords *ray_velocity, Coords* ray_position_geometry, - Coords *ray_position, - struct Volume_struct* Volume, struct focus_data_struct* this_focus_data){ + void + choose_scattering_point_and_ray_direction (double* forced_length_to_scattering, double* safety_distance, double* safety_distance2, double* length_to_boundary, + _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, + struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { // Sample length_to_scattering in linear manner *forced_length_to_scattering = *safety_distance + rand01 () * (*length_to_boundary - *safety_distance2); @@ -127,8 +121,7 @@ SHARE Coords direction = coords_scalar_mult (*ray_velocity, 1.0 / length_of_position_vector (*ray_velocity)); Coords scattering_displacement = coords_scalar_mult (direction, *forced_length_to_scattering); Coords forced_ray_scattering_point = coords_add (*ray_position, scattering_displacement); - *ray_position_geometry - = coords_sub (forced_ray_scattering_point, Volume->geometry.center); // ray_position relative to geometry center + *ray_position_geometry = coords_sub (forced_ray_scattering_point, Volume->geometry.center); // ray_position relative to geometry center // Calculate the aim for non isotropic processes this_focus_data = &Volume->geometry.focus_data_array.elements[0]; @@ -145,98 +138,86 @@ SHARE printf ("for isotropic processes this RayAim is used \n"); print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); #endif + } + void + transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double (*k_rotated)[3], int* p_index, + Coords* wavevector, Coords* ray_position_geometry) { - } + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; + *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); + coords_get (*wavevector_rotated, &(*k_rotated)[0], &(*k_rotated)[1], &(*k_rotated)[2]); + if (Volume->p_physics->p_scattering_array[*p_index].needs_cross_section_focus == 1) { + // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time + // focusing + Coords ray_position_geometry_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *ray_position_geometry); - void transform_wavevector_into_local_coord_system(struct Volume_struct* Volume, - Coords * wavevector_rotated, double (*k_rotated)[3], - int *p_index, Coords* wavevector, Coords* ray_position_geometry) { - - int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; - *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); - coords_get (*wavevector_rotated, &(*k_rotated)[0], &(*k_rotated)[1], &(*k_rotated)[2]); - - if (Volume->p_physics->p_scattering_array[*p_index].needs_cross_section_focus == 1) { - // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time - // focusing - Coords ray_position_geometry_rotated - = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *ray_position_geometry); - - int focus_data_index = Volume->geometry.focus_array_indices.elements[*p_index]; - struct focus_data_struct* this_focus_data = &Volume->geometry.focus_data_array.elements[focus_data_index]; - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray - - #ifdef Union_trace_verbal_setting - printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", *p_index); - print_position (*ray_position_geometry, "ray_position_geometry"); - print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); - print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - #endif - } - } + int focus_data_index = Volume->geometry.focus_array_indices.elements[*p_index]; + struct focus_data_struct* this_focus_data = &Volume->geometry.focus_data_array.elements[focus_data_index]; + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray + #ifdef Union_trace_verbal_setting + printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", *p_index); + print_position (*ray_position_geometry, "ray_position_geometry"); + print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); + print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + #endif + } + } - int process_needs_inhomogenous_sampling(struct physics_struct* current_p_physics, - struct scattering_process_struct* process){ - if (current_p_physics->sampling_points !=0 && process->needs_cross_section_focus){ + int + process_needs_inhomogenous_sampling (struct physics_struct* current_p_physics, struct scattering_process_struct* process) { + if (current_p_physics->sampling_points != 0 && process->needs_cross_section_focus) { if (process->sampling_points != -1) return 1; } return 0; } - - void set_cumul_dist_array(struct physics_struct* current_p_physics, int 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; + void + set_cumul_dist_array (struct physics_struct* current_p_physics, int 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; } - void move_and_aim_neutron(struct physics_struct* current_p_physics, int i, - struct scattering_process_struct* process, - _class_particle* _particle, - struct Volume_struct* Volume, int* p_index, - Coords* ray_velocity, Coords* ray_position, - struct focus_data_struct* this_focus_data - ){ - // Transport neutron to place inside geometry - *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->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, Volume->geometry.center); - // Also focus the ray at this point, if the component needs focusing - if (process->needs_cross_section_focus) { - if (Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index != -1) { - int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; - sampling_point_geometry - = rot_apply (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); + void + move_and_aim_neutron (struct physics_struct* current_p_physics, int i, struct scattering_process_struct* process, _class_particle* _particle, + struct Volume_struct* Volume, int* p_index, Coords* ray_velocity, Coords* ray_position, struct focus_data_struct* this_focus_data) { + // Transport neutron to place inside geometry + *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->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, Volume->geometry.center); + // Also focus the ray at this point, if the component needs focusing + if (process->needs_cross_section_focus) { + if (Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index != -1) { + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; + sampling_point_geometry = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], sampling_point_geometry); } - // Calculate mu and probability - coords_get (sampling_point_geometry, &_particle->x, &_particle->y, &_particle->z); + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, sampling_point_geometry); + } + // Calculate mu and probability + coords_get (sampling_point_geometry, &_particle->x, &_particle->y, &_particle->z); } - - int mu_and_intersect_dist_safeguard(double mu_sum, double length_to_boundary, double safety_distance2, - int *scattering_event){ - if (mu_sum < 1E-18){ + int + mu_and_intersect_dist_safeguard (double mu_sum, double length_to_boundary, double safety_distance2, int* scattering_event) { + if (mu_sum < 1E-18) { scattering_event = 0; - return 0; + return 0; } - if (length_to_boundary < safety_distance2){ + if (length_to_boundary < safety_distance2) { scattering_event = 0; return 0; } return 1; } - void sample_inhomogenous_transmission_probability(struct physics_struct* current_p_physics, - struct Volume_struct* Volume, double* real_transmission_probability, - double v_length - ){ + void + sample_inhomogenous_transmission_probability (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* real_transmission_probability, + double v_length) { // Calculate the probabilities and then add them cumulatively memset (current_p_physics->total_mus, 0, sizeof (double) * current_p_physics->sampling_points); @@ -270,12 +251,9 @@ SHARE *real_transmission_probability = current_p_physics->cumul_transmission_prob[current_p_physics->sampling_points - 1]; } - double sample_scattering_point_inhomogenous(struct Volume_struct* Volume, - struct physics_struct* current_p_physics, - double* abs_weight_factor, - double* v_length, - double* safety_distance, - int* selected_sampling){ + double + sample_scattering_point_inhomogenous (struct Volume_struct* Volume, struct physics_struct* current_p_physics, double* abs_weight_factor, double* v_length, + double* safety_distance, int* selected_sampling) { // 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, @@ -290,52 +268,42 @@ SHARE *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]; + *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; + 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; return current_p_physics->cumul_dists[*selected_sampling] - current_p_physics->dist / 2 + sampled_dist; } - - int p_interact_is_set(struct Volume_struct* Volume){ - if (Volume->geometry.geometry_p_interact != 0){ + int + p_interact_is_set (struct Volume_struct* Volume) { + if (Volume->geometry.geometry_p_interact != 0) { return 1; } return 0; } - void p_interact_check_scattering_event(struct Volume_struct* Volume, - int* scattering_event, - double* weight, - double real_transmission_prob){ + void + p_interact_check_scattering_event (struct Volume_struct* Volume, int* scattering_event, double* weight, double real_transmission_prob) { double mc_transmission_prob = 1 - Volume->geometry.geometry_p_interact; - *scattering_event = rand01() > mc_transmission_prob; - if (*scattering_event){ + *scattering_event = rand01 () > mc_transmission_prob; + if (*scattering_event) { // Scattering event happens, this is the correction for the weight - *weight *= (1.0 - real_transmission_prob)/(1.0 - mc_transmission_prob); - } - else { + *weight *= (1.0 - real_transmission_prob) / (1.0 - mc_transmission_prob); + } else { // Scattering event does not happen, this is the appropriate correction *weight *= real_transmission_prob / mc_transmission_prob; } } - - void p_interact_select_process(struct Volume_struct* Volume, - double* my_trace_fraction_control, - double* my_trace, - double* total_process_interact, - double* culmative_probability, - double* mc_prop, - double* weight, - double* my_sum, - int* selected_process - ){ + void + p_interact_select_process (struct Volume_struct* Volume, double* my_trace_fraction_control, double* my_trace, double* total_process_interact, + double* culmative_probability, double* mc_prop, double* weight, double* my_sum, int* selected_process) { // Interact_fraction is used to influence the choice of process in this material *mc_prop = rand01 (); *culmative_probability = 0; @@ -362,13 +330,9 @@ SHARE } } - void inhomogenous_choose_process(struct physics_struct* current_p_physics, - struct Volume_struct* Volume, - double* culmative_probability, - double* mc_prop, - double* my_sum, - int* selected_sampling, - int* selected_process){ + void + inhomogenous_choose_process (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* culmative_probability, double* mc_prop, + double* my_sum, int* selected_sampling, int* selected_process) { for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { *culmative_probability += current_p_physics->mus[i][*selected_sampling] / *my_sum; if (*culmative_probability > *mc_prop) { @@ -378,26 +342,19 @@ SHARE } } - void dir_focus_set_scat_length(double* length_to_scattering, - double* forced_length_to_scattering, - double* weight, - double* length_to_boundary, - double* my_sum_plus_abs){ + void + dir_focus_set_scat_length (double* length_to_scattering, double* forced_length_to_scattering, double* weight, double* length_to_boundary, + double* my_sum_plus_abs) { // 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)) - *weight *= *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)); + *weight *= *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 - } - - %} DECLARE @@ -1858,11 +1815,10 @@ TRACE scattering_event = 0; // Assume a scattering event will not occur // Check if a scattering event should occur - if (current_volume != 0) { // Volume 0 is always vacuum, and if this is the current volume, an event will not occur - if (volume_is_only_absorber(Volumes[current_volume])) { // If there are no processes, the volume could be vacuum or an absorber - adjust_abs_weight_factor(Volumes[current_volume], &my_sum_plus_abs, - &length_to_boundary, &v_length, &time_to_boundery, - &abs_weight_factor, &abs_weight_factor_set); + if (current_volume != 0) { // Volume 0 is always vacuum, and if this is the current volume, an event will not occur + if (volume_is_only_absorber (Volumes[current_volume])) { // If there are no processes, the volume could be vacuum or an absorber + adjust_abs_weight_factor (Volumes[current_volume], &my_sum_plus_abs, &length_to_boundary, &v_length, &time_to_boundery, &abs_weight_factor, + &abs_weight_factor_set); } 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; @@ -1878,12 +1834,8 @@ TRACE length_to_boundary = time_to_boundery * v_length; // If any process in this material needs focusing, sample scattering position and update focus_data accordingly if (current_p_physics->any_process_needs_cross_section_focus == 1) { - choose_scattering_point_and_ray_direction(&forced_length_to_scattering, - &safety_distance, &safety_distance2, - &length_to_boundary, _particle, - &ray_velocity, &ray_position_geometry, - &ray_position, - Volumes[current_volume], this_focus_data); + choose_scattering_point_and_ray_direction (&forced_length_to_scattering, &safety_distance, &safety_distance2, &length_to_boundary, _particle, + &ray_velocity, &ray_position_geometry, &ray_position, Volumes[current_volume], this_focus_data); } else { forced_length_to_scattering = -1.0; // Signals that no forcing needed, could also if on the selected process struct } @@ -1897,9 +1849,8 @@ TRACE if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { // If the process is not isotropic, the wavevector is transformed into the local coordinate system of the process - transform_wavevector_into_local_coord_system(Volumes[current_volume], - &wavevector_rotated, &k_rotated, - &p_index, &wavevector, &ray_position_geometry); + transform_wavevector_into_local_coord_system (Volumes[current_volume], &wavevector_rotated, &k_rotated, &p_index, &wavevector, + &ray_position_geometry); } else { k_rotated[0] = k[0]; k_rotated[1] = k[1]; @@ -1914,27 +1865,23 @@ TRACE double mu; current_p_physics->dist = length_to_boundary / current_p_physics->sampling_points - safety_distance; - if (process_needs_inhomogenous_sampling(current_p_physics, process)) { + if (process_needs_inhomogenous_sampling (current_p_physics, process)) { // If mulitple sampling is required for this process, populate the neutron distance and probability arrays: Coords original_position = coords_set (x, y, z); double mu; *p_my_trace = 0; this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; for (int i = 0; i < current_p_physics->sampling_points; i++) { - set_cumul_dist_array(current_p_physics, i); - move_and_aim_neutron(current_p_physics, i, process, _particle, - Volumes[current_volume], &p_index, - &ray_velocity, &ray_position, - this_focus_data); + set_cumul_dist_array (current_p_physics, i); + move_and_aim_neutron (current_p_physics, i, process, _particle, Volumes[current_volume], &p_index, &ray_velocity, &ray_position, this_focus_data); // Calculate mu and probability physics_my (process->eProcess, &mu, k_rotated, process->data_transfer, this_focus_data, _particle); current_p_physics->mus[p_index][i] = mu; *p_my_trace += mu / current_p_physics->sampling_points; } coords_get (original_position, &x, &y, &z); - - } - if (!process_needs_inhomogenous_sampling(current_p_physics, process)){ + } + if (!process_needs_inhomogenous_sampling (current_p_physics, process)) { 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; @@ -1961,25 +1908,19 @@ TRACE // New flow:length_to_boundary // Calculate if scattering happens based on my_sub_plus_abs - if (mu_and_intersect_dist_safeguard(my_sum, length_to_boundary, safety_distance2, - &scattering_event)){ + if (mu_and_intersect_dist_safeguard (my_sum, length_to_boundary, safety_distance2, &scattering_event)) { // First calculate the transmission probability if (current_p_physics->sampling_points == 0) { real_transmission_probability = exp (-length_to_boundary * my_sum_plus_abs); - } + } if (current_p_physics->sampling_points != 0) { - sample_inhomogenous_transmission_probability(current_p_physics, - Volumes[current_volume], - &real_transmission_probability, - v_length); + sample_inhomogenous_transmission_probability (current_p_physics, Volumes[current_volume], &real_transmission_probability, v_length); } // Then check if we scatter. - if (p_interact_is_set(Volumes[current_volume])) { - p_interact_check_scattering_event(Volumes[current_volume], - &scattering_event, - &p, real_transmission_probability); + if (p_interact_is_set (Volumes[current_volume])) { + p_interact_check_scattering_event (Volumes[current_volume], &scattering_event, &p, real_transmission_probability); } else { // probability to scatter is the natural value scattering_event = rand01 () > real_transmission_probability; @@ -1998,12 +1939,8 @@ TRACE printf ("WARNING: Absorption weight factor above 1! Should not happen! \n"); // Select distance to scattering position if (current_p_physics->sampling_points != 0) { - length_to_scattering = sample_scattering_point_inhomogenous(Volumes[current_volume], - current_p_physics, - &abs_weight_factor, - &v_length, - &safety_distance, - &selected_sampling); + length_to_scattering = sample_scattering_point_inhomogenous (Volumes[current_volume], current_p_physics, &abs_weight_factor, &v_length, + &safety_distance, &selected_sampling); } // Select process if (Volumes[current_volume]->p_physics->number_of_processes == 1) { // trivial case @@ -2023,25 +1960,12 @@ TRACE } } } else { - inhomogenous_choose_process(current_p_physics, - Volumes[current_volume], - &culmative_probability, - &mc_prop, - &my_sum, - &selected_sampling, - &selected_process); - + inhomogenous_choose_process (current_p_physics, Volumes[current_volume], &culmative_probability, &mc_prop, &my_sum, &selected_sampling, + &selected_process); } } else { - p_interact_select_process(Volumes[current_volume], - my_trace_fraction_control, - my_trace, - &total_process_interact, - &culmative_probability, - &mc_prop, - &p, - &my_sum, - &selected_process); + p_interact_select_process (Volumes[current_volume], my_trace_fraction_control, my_trace, &total_process_interact, &culmative_probability, + &mc_prop, &p, &my_sum, &selected_process); } } @@ -2049,9 +1973,7 @@ TRACE if (current_p_physics->sampling_points == 0) { // No numerical integration is necessary. if (process->needs_cross_section_focus == 1) { - dir_focus_set_scat_length(&length_to_scattering, &forced_length_to_scattering, - &p, &length_to_boundary, - &my_sum_plus_abs); + dir_focus_set_scat_length (&length_to_scattering, &forced_length_to_scattering, &p, &length_to_boundary, &my_sum_plus_abs); } else { // Decided the ray scatters, choose where on truncated exponential from safety_distance to length_to_boundary - safety_distance length_to_scattering @@ -2114,7 +2036,7 @@ TRACE // Logging for detector components assosiated with this volume for (log_index = 0; log_index < Volumes[current_volume]->abs_loggers.num_elements; log_index++) { // This function calls a logger function which in turn stores some data among the passed, and possibly performs some basic data analysis - // Position and k_new given in master coordinates, the abs_logger must transform to its coordnate system if required + // Position and k_new given in master coordinates, the abs_logger must transform to its coordnate system if required Volumes[current_volume]->abs_loggers.p_abs_logger[log_index]->function_pointers.active_record_function ( &abs_position, k_new, initial_weight * (1.0 - abs_weight_factor), t + t_abs_propagation, scattered_flag[current_volume], number_of_scattering_events, Volumes[current_volume]->abs_loggers.p_abs_logger[log_index], &abs_loggers_with_data_array); From 12041c4123ddde0be29b84ad6abac997217c671f Mon Sep 17 00:00:00 2001 From: Diablo Date: Sun, 2 Aug 2026 11:28:48 +0200 Subject: [PATCH 13/36] Rename choose scattering point and direction, to reflect that it is only for the focusing processes --- mcstas-comps/union/Union_master.comp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index b79b7e80b1..1d73f0d485 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -110,7 +110,7 @@ SHARE } void - choose_scattering_point_and_ray_direction (double* forced_length_to_scattering, double* safety_distance, double* safety_distance2, double* length_to_boundary, + focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double* safety_distance, double* safety_distance2, double* length_to_boundary, _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { // Sample length_to_scattering in linear manner @@ -1834,7 +1834,7 @@ TRACE length_to_boundary = time_to_boundery * v_length; // If any process in this material needs focusing, sample scattering position and update focus_data accordingly if (current_p_physics->any_process_needs_cross_section_focus == 1) { - choose_scattering_point_and_ray_direction (&forced_length_to_scattering, &safety_distance, &safety_distance2, &length_to_boundary, _particle, + focus_in_cross_section_set_forced_point_and_dir (&forced_length_to_scattering, &safety_distance, &safety_distance2, &length_to_boundary, _particle, &ray_velocity, &ray_position_geometry, &ray_position, Volumes[current_volume], this_focus_data); } else { forced_length_to_scattering = -1.0; // Signals that no forcing needed, could also if on the selected process struct From abc75144c0317a35cf48d6b4c51a7b0639ff7e04 Mon Sep 17 00:00:00 2001 From: Diablo Date: Sun, 2 Aug 2026 11:43:26 +0200 Subject: [PATCH 14/36] fix erroneus if statements inside process_needs_inhomogenous sampling --- mcstas-comps/union/Union_master.comp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 1d73f0d485..3560dbaa13 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -168,8 +168,8 @@ SHARE int process_needs_inhomogenous_sampling (struct physics_struct* current_p_physics, struct scattering_process_struct* process) { - if (current_p_physics->sampling_points != 0 && process->needs_cross_section_focus) { - if (process->sampling_points != -1) + if (current_p_physics->sampling_points != 0) { + if (process->needs_cross_section_focus || process->sampling_points != -1) return 1; } return 0; From c25937ca97fe9740712b04d06a6a00141eba8a63 Mon Sep 17 00:00:00 2001 From: Diablo Date: Sun, 2 Aug 2026 11:44:39 +0200 Subject: [PATCH 15/36] Add pointer dereference to scattering_event inside safeguard function --- .../Tests_union/Test_inhomogenous_process/Union_master.comp | 1 + 1 file changed, 1 insertion(+) create mode 120000 mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Union_master.comp diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Union_master.comp b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Union_master.comp new file mode 120000 index 0000000000..2c17b361c4 --- /dev/null +++ b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Union_master.comp @@ -0,0 +1 @@ +../../../union/Union_master.comp \ No newline at end of file From 0e155de48b31d6991baa355f4d1edfdd314d56b9 Mon Sep 17 00:00:00 2001 From: Diablo Date: Sun, 2 Aug 2026 11:46:03 +0200 Subject: [PATCH 16/36] Remove erroneously placed symbolic link to Union master --- .../Tests_union/Test_inhomogenous_process/Union_master.comp | 1 - 1 file changed, 1 deletion(-) delete mode 120000 mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Union_master.comp diff --git a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Union_master.comp b/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Union_master.comp deleted file mode 120000 index 2c17b361c4..0000000000 --- a/mcstas-comps/examples/Tests_union/Test_inhomogenous_process/Union_master.comp +++ /dev/null @@ -1 +0,0 @@ -../../../union/Union_master.comp \ No newline at end of file From 3f8a4f08c972b097f68fa2efcc5d3b9dd10455ec Mon Sep 17 00:00:00 2001 From: Diablo Date: Sun, 2 Aug 2026 11:46:29 +0200 Subject: [PATCH 17/36] Add pointer dereference to scattering_event inside safeguard function. Now in correct file --- mcstas-comps/union/Union_master.comp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 3560dbaa13..27c051e56e 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -205,11 +205,11 @@ SHARE int mu_and_intersect_dist_safeguard (double mu_sum, double length_to_boundary, double safety_distance2, int* scattering_event) { if (mu_sum < 1E-18) { - scattering_event = 0; + *scattering_event = 0; return 0; } if (length_to_boundary < safety_distance2) { - scattering_event = 0; + *scattering_event = 0; return 0; } return 1; From e45afc6805929765e83fff17b8d5179a7c69893e Mon Sep 17 00:00:00 2001 From: Diablo Date: Sun, 2 Aug 2026 11:50:32 +0200 Subject: [PATCH 18/36] Change debug statement from using Volumes array to using correct Volume struct --- mcstas-comps/union/Union_master.comp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 27c051e56e..224f4b93e0 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -102,10 +102,10 @@ SHARE *abs_weight_factor_set = 1; #ifdef Union_trace_verbal_setting - printf ("name of material: %s \n", Volumes->name); + printf ("name of material: %s \n", Volume->name); printf ("length to boundery = %f\n", length_to_boundary); - printf ("absorption cross section = %f\n", Volumes->p_physics->my_a); - printf ("chance to get through this length of absorber: %f %%\n", 100 * exp (-Volumes->p_physics->my_a * length_to_boundary)); + printf ("absorption cross section = %f\n", Volume->p_physics->my_a); + printf ("chance to get through this length of absorber: %f %%\n", 100 * exp (-Volume->p_physics->my_a * length_to_boundary)); #endif } From c29d82f5d66eaaaa4477f2d3d6c491fe4f8d6136 Mon Sep 17 00:00:00 2001 From: Diablo Date: Sun, 2 Aug 2026 12:10:04 +0200 Subject: [PATCH 19/36] regularize naming conventions across the functions for the logic path. Also move the functions toward their relevant others, i.e inhomogenous with inhomogenous --- mcstas-comps/union/Union_master.comp | 155 ++++++++++++++------------- 1 file changed, 82 insertions(+), 73 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 224f4b93e0..fdf012a956 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -92,6 +92,16 @@ SHARE return 1; return 0; } + + int + process_needs_inhomogenous_sampling (struct physics_struct* current_p_physics, struct scattering_process_struct* process) { + if (current_p_physics->sampling_points != 0) { + if (process->needs_cross_section_focus || process->sampling_points != -1) + return 1; + } + return 0; + } + void adjust_abs_weight_factor (struct Volume_struct* Volume, double* my_sum_plus_abs, double* length_to_boundary, double* v_length, double* time_to_boundary, double* abs_weight_factor, int* abs_weight_factor_set) { @@ -109,36 +119,6 @@ SHARE #endif } - void - focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double* safety_distance, double* safety_distance2, double* length_to_boundary, - _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, - struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { - // Sample length_to_scattering in linear manner - *forced_length_to_scattering = *safety_distance + rand01 () * (*length_to_boundary - *safety_distance2); - - *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); // Test for root cause - // 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 scattering_displacement = coords_scalar_mult (direction, *forced_length_to_scattering); - Coords forced_ray_scattering_point = coords_add (*ray_position, scattering_displacement); - *ray_position_geometry = coords_sub (forced_ray_scattering_point, Volume->geometry.center); // ray_position relative to geometry center - - // Calculate the aim for non isotropic processes - this_focus_data = &Volume->geometry.focus_data_array.elements[0]; - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, *ray_position_geometry); // Aim vector for this ray - - #ifdef Union_trace_verbal_setting - printf ("Prepared for focus in cross section calculation in volume: %s \n", Volume->name); - printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); - print_position (*ray_position, "ray_position"); - print_position (direction, "direction"); - print_position (scattering_displacement, "scattering_displacement"); - print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); - print_position (*ray_position_geometry, "ray_position_geometry"); - printf ("for isotropic processes this RayAim is used \n"); - print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - #endif - } void transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double (*k_rotated)[3], int* p_index, @@ -166,19 +146,7 @@ SHARE } } - int - process_needs_inhomogenous_sampling (struct physics_struct* current_p_physics, struct scattering_process_struct* process) { - if (current_p_physics->sampling_points != 0) { - if (process->needs_cross_section_focus || process->sampling_points != -1) - return 1; - } - return 0; - } - void - set_cumul_dist_array (struct physics_struct* current_p_physics, int 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; - } void move_and_aim_neutron (struct physics_struct* current_p_physics, int i, struct scattering_process_struct* process, _class_particle* _particle, @@ -215,8 +183,60 @@ SHARE return 1; } + void - sample_inhomogenous_transmission_probability (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* real_transmission_probability, + focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double* safety_distance, double* safety_distance2, double* length_to_boundary, + _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, + struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { + // Sample length_to_scattering in linear manner + *forced_length_to_scattering = *safety_distance + rand01 () * (*length_to_boundary - *safety_distance2); + + *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); // Test for root cause + // 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 scattering_displacement = coords_scalar_mult (direction, *forced_length_to_scattering); + Coords forced_ray_scattering_point = coords_add (*ray_position, scattering_displacement); + *ray_position_geometry = coords_sub (forced_ray_scattering_point, Volume->geometry.center); // ray_position relative to geometry center + + // Calculate the aim for non isotropic processes + this_focus_data = &Volume->geometry.focus_data_array.elements[0]; + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, *ray_position_geometry); // Aim vector for this ray + + #ifdef Union_trace_verbal_setting + printf ("Prepared for focus in cross section calculation in volume: %s \n", Volume->name); + printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); + print_position (*ray_position, "ray_position"); + print_position (direction, "direction"); + print_position (scattering_displacement, "scattering_displacement"); + print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); + print_position (*ray_position_geometry, "ray_position_geometry"); + printf ("for isotropic processes this RayAim is used \n"); + print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + #endif + } + + void + focus_in_cross_section_set_scat_length (double* length_to_scattering, double* forced_length_to_scattering, double* weight, double* length_to_boundary, + double* my_sum_plus_abs) { + // 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)) + *weight *= *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 + } + + + void + inhomogenous_set_cumul_dist_array (struct physics_struct* current_p_physics, int 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; + } + + void + inhomogenous_sample_transmission_probability (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* real_transmission_probability, double v_length) { // Calculate the probabilities and then add them cumulatively memset (current_p_physics->total_mus, 0, sizeof (double) * current_p_physics->sampling_points); @@ -252,7 +272,7 @@ SHARE } double - sample_scattering_point_inhomogenous (struct Volume_struct* Volume, struct physics_struct* current_p_physics, double* abs_weight_factor, double* v_length, + inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* current_p_physics, double* abs_weight_factor, double* v_length, double* safety_distance, int* selected_sampling) { // Numerical integration happens, and therefore we must choose between the different samples @@ -280,6 +300,20 @@ SHARE return current_p_physics->cumul_dists[*selected_sampling] - current_p_physics->dist / 2 + sampled_dist; } + + void + inhomogenous_choose_process (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* culmative_probability, double* mc_prop, + double* my_sum, int* selected_sampling, int* selected_process) { + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + *culmative_probability += current_p_physics->mus[i][*selected_sampling] / *my_sum; + if (*culmative_probability > *mc_prop) { + *selected_process = i; + break; + } + } + } + + int p_interact_is_set (struct Volume_struct* Volume) { if (Volume->geometry.geometry_p_interact != 0) { @@ -330,31 +364,6 @@ SHARE } } - void - inhomogenous_choose_process (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* culmative_probability, double* mc_prop, - double* my_sum, int* selected_sampling, int* selected_process) { - for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { - *culmative_probability += current_p_physics->mus[i][*selected_sampling] / *my_sum; - if (*culmative_probability > *mc_prop) { - *selected_process = i; - break; - } - } - } - - void - dir_focus_set_scat_length (double* length_to_scattering, double* forced_length_to_scattering, double* weight, double* length_to_boundary, - double* my_sum_plus_abs) { - // 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)) - *weight *= *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 - } %} DECLARE @@ -1872,7 +1881,7 @@ TRACE *p_my_trace = 0; this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; for (int i = 0; i < current_p_physics->sampling_points; i++) { - set_cumul_dist_array (current_p_physics, i); + inhomogenous_set_cumul_dist_array (current_p_physics, i); move_and_aim_neutron (current_p_physics, i, process, _particle, Volumes[current_volume], &p_index, &ray_velocity, &ray_position, this_focus_data); // Calculate mu and probability physics_my (process->eProcess, &mu, k_rotated, process->data_transfer, this_focus_data, _particle); @@ -1915,7 +1924,7 @@ TRACE } if (current_p_physics->sampling_points != 0) { - sample_inhomogenous_transmission_probability (current_p_physics, Volumes[current_volume], &real_transmission_probability, v_length); + inhomogenous_sample_transmission_probability (current_p_physics, Volumes[current_volume], &real_transmission_probability, v_length); } // Then check if we scatter. @@ -1939,7 +1948,7 @@ TRACE printf ("WARNING: Absorption weight factor above 1! Should not happen! \n"); // Select distance to scattering position if (current_p_physics->sampling_points != 0) { - length_to_scattering = sample_scattering_point_inhomogenous (Volumes[current_volume], current_p_physics, &abs_weight_factor, &v_length, + length_to_scattering = inhomogenous_sample_scattering_point (Volumes[current_volume], current_p_physics, &abs_weight_factor, &v_length, &safety_distance, &selected_sampling); } // Select process @@ -1973,7 +1982,7 @@ TRACE if (current_p_physics->sampling_points == 0) { // No numerical integration is necessary. if (process->needs_cross_section_focus == 1) { - dir_focus_set_scat_length (&length_to_scattering, &forced_length_to_scattering, &p, &length_to_boundary, &my_sum_plus_abs); + focus_in_cross_section_set_scat_length (&length_to_scattering, &forced_length_to_scattering, &p, &length_to_boundary, &my_sum_plus_abs); } else { // Decided the ray scatters, choose where on truncated exponential from safety_distance to length_to_boundary - safety_distance length_to_scattering From 8f0a6063a0feb5d4938ddc1caf7936536224fceb Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 16 Sep 2026 13:38:07 +0200 Subject: [PATCH 20/36] Move the helper functions from Union master to union-lib.c --- mcstas-comps/share/union-lib.c | 283 +++++++++++++++++++++++++++ mcstas-comps/union/Union_master.comp | 279 -------------------------- 2 files changed, 283 insertions(+), 279 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 8e6e6176df..dce01ed98b 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -9221,3 +9221,286 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } } +//============================================================================== +//=========== Helper functions to Union_master ================================= +//============================================================================== + int + volume_is_only_absorber (struct Volume_struct* Volume) { + // This function returns true if a volume does not have any physical processes + // and if the volume is not a vacuum. + printf("TESTING TESTING \n"); + if (!Volume->p_physics->number_of_processes && !Volume->p_physics->is_vacuum) + return 1; + return 0; + } + + int + process_needs_inhomogenous_sampling (struct physics_struct* current_p_physics, struct scattering_process_struct* process) { + if (current_p_physics->sampling_points != 0) { + if (process->needs_cross_section_focus || process->sampling_points != -1) + return 1; + } + return 0; + } + + void + adjust_abs_weight_factor (struct Volume_struct* Volume, double* my_sum_plus_abs, double* length_to_boundary, double* v_length, double* time_to_boundary, + double* abs_weight_factor, int* abs_weight_factor_set) { + *my_sum_plus_abs = Volume->p_physics->my_a * (2200 / *v_length); + *length_to_boundary = *time_to_boundary * *v_length; + + *abs_weight_factor = exp (-Volume->p_physics->my_a * 2200 * *time_to_boundary); + *abs_weight_factor_set = 1; + + #ifdef Union_trace_verbal_setting + printf ("name of material: %s \n", Volume->name); + printf ("length to boundery = %f\n", length_to_boundary); + printf ("absorption cross section = %f\n", Volume->p_physics->my_a); + printf ("chance to get through this length of absorber: %f %%\n", 100 * exp (-Volume->p_physics->my_a * length_to_boundary)); + #endif + } + + + void + transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double (*k_rotated)[3], int* p_index, + Coords* wavevector, Coords* ray_position_geometry) { + + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; + *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); + coords_get (*wavevector_rotated, &(*k_rotated)[0], &(*k_rotated)[1], &(*k_rotated)[2]); + + if (Volume->p_physics->p_scattering_array[*p_index].needs_cross_section_focus == 1) { + // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time + // focusing + Coords ray_position_geometry_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *ray_position_geometry); + + int focus_data_index = Volume->geometry.focus_array_indices.elements[*p_index]; + struct focus_data_struct* this_focus_data = &Volume->geometry.focus_data_array.elements[focus_data_index]; + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray + + #ifdef Union_trace_verbal_setting + printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", *p_index); + print_position (*ray_position_geometry, "ray_position_geometry"); + print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); + print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + #endif + } + } + + + + void + move_and_aim_neutron (struct physics_struct* current_p_physics, int i, struct scattering_process_struct* process, _class_particle* _particle, + struct Volume_struct* Volume, int* p_index, Coords* ray_velocity, Coords* ray_position, struct focus_data_struct* this_focus_data) { + // Transport neutron to place inside geometry + *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->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, Volume->geometry.center); + // Also focus the ray at this point, if the component needs focusing + if (process->needs_cross_section_focus) { + if (Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index != -1) { + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; + sampling_point_geometry = rot_apply (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, &_particle->x, &_particle->y, &_particle->z); + } + + int + mu_and_intersect_dist_safeguard (double mu_sum, double length_to_boundary, double safety_distance2, int* scattering_event) { + if (mu_sum < 1E-18) { + *scattering_event = 0; + return 0; + } + if (length_to_boundary < safety_distance2) { + *scattering_event = 0; + return 0; + } + return 1; + } + + + void + focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double* safety_distance, double* safety_distance2, double* length_to_boundary, + _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, + struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { + // Sample length_to_scattering in linear manner + *forced_length_to_scattering = *safety_distance + rand01 () * (*length_to_boundary - *safety_distance2); + + *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); // Test for root cause + // 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 scattering_displacement = coords_scalar_mult (direction, *forced_length_to_scattering); + Coords forced_ray_scattering_point = coords_add (*ray_position, scattering_displacement); + *ray_position_geometry = coords_sub (forced_ray_scattering_point, Volume->geometry.center); // ray_position relative to geometry center + + // Calculate the aim for non isotropic processes + this_focus_data = &Volume->geometry.focus_data_array.elements[0]; + this_focus_data->RayAim = coords_sub (this_focus_data->Aim, *ray_position_geometry); // Aim vector for this ray + + #ifdef Union_trace_verbal_setting + printf ("Prepared for focus in cross section calculation in volume: %s \n", Volume->name); + printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); + print_position (*ray_position, "ray_position"); + print_position (direction, "direction"); + print_position (scattering_displacement, "scattering_displacement"); + print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); + print_position (*ray_position_geometry, "ray_position_geometry"); + printf ("for isotropic processes this RayAim is used \n"); + print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); + #endif + } + + void + focus_in_cross_section_set_scat_length (double* length_to_scattering, double* forced_length_to_scattering, double* weight, double* length_to_boundary, + double* my_sum_plus_abs) { + // 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)) + *weight *= *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 + } + + + void + inhomogenous_set_cumul_dist_array (struct physics_struct* current_p_physics, int 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; + } + + void + inhomogenous_sample_transmission_probability (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* real_transmission_probability, + double v_length) { + // 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 < Volume->p_physics->number_of_processes; i++) { + struct scattering_process_struct* process_i = &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 = 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]; + } + + double + inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* current_p_physics, double* abs_weight_factor, double* v_length, + double* safety_distance, int* selected_sampling) { + + // 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 = 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; + return current_p_physics->cumul_dists[*selected_sampling] - current_p_physics->dist / 2 + sampled_dist; + } + + + void + inhomogenous_choose_process (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* culmative_probability, double* mc_prop, + double* my_sum, int* selected_sampling, int* selected_process) { + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + *culmative_probability += current_p_physics->mus[i][*selected_sampling] / *my_sum; + if (*culmative_probability > *mc_prop) { + *selected_process = i; + break; + } + } + } + + + int + p_interact_is_set (struct Volume_struct* Volume) { + if (Volume->geometry.geometry_p_interact != 0) { + return 1; + } + return 0; + } + + void + p_interact_check_scattering_event (struct Volume_struct* Volume, int* scattering_event, double* weight, double real_transmission_prob) { + double mc_transmission_prob = 1 - Volume->geometry.geometry_p_interact; + *scattering_event = rand01 () > mc_transmission_prob; + if (*scattering_event) { + // Scattering event happens, this is the correction for the weight + *weight *= (1.0 - real_transmission_prob) / (1.0 - mc_transmission_prob); + } else { + // Scattering event does not happen, this is the appropriate correction + *weight *= real_transmission_prob / mc_transmission_prob; + } + } + + void + p_interact_select_process (struct Volume_struct* Volume, double* my_trace_fraction_control, double* my_trace, double* total_process_interact, + double* culmative_probability, double* mc_prop, double* weight, double* my_sum, int* selected_process) { + // Interact_fraction is used to influence the choice of process in this material + *mc_prop = rand01 (); + *culmative_probability = 0; + *total_process_interact = 1.0; + + // If any of the processes have probability 0, they are excluded from the selection + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + if (my_trace[i] < 1E-18) { + // When this happens, the total force probability is corrected and the probability for this particular instance is set to 0 + *total_process_interact -= Volume->p_physics->p_scattering_array[i].process_p_interact; + my_trace_fraction_control[i] = 0; + // In cases where my_trace is not zero, the forced fraction is still used. + } else + my_trace_fraction_control[i] = Volume->p_physics->p_scattering_array[i].process_p_interact; + } + // Randomly select a process using the weights stored in my_trace_fraction_control divided by total_process_interact + for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { + *culmative_probability += my_trace_fraction_control[i] / *total_process_interact; + if (*culmative_probability > *mc_prop) { + *selected_process = i; + *weight *= (my_trace[i] / *my_sum) * (*total_process_interact / my_trace_fraction_control[i]); + break; + } + } + } diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index fdf012a956..1f930b697d 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -84,285 +84,6 @@ SHARE #define MASTER_DETECTOR dummy #endif - int - volume_is_only_absorber (struct Volume_struct* Volume) { - // This function returns true if a volume does not have any physical processes - // and if the volume is not a vacuum. - if (!Volume->p_physics->number_of_processes && !Volume->p_physics->is_vacuum) - return 1; - return 0; - } - - int - process_needs_inhomogenous_sampling (struct physics_struct* current_p_physics, struct scattering_process_struct* process) { - if (current_p_physics->sampling_points != 0) { - if (process->needs_cross_section_focus || process->sampling_points != -1) - return 1; - } - return 0; - } - - void - adjust_abs_weight_factor (struct Volume_struct* Volume, double* my_sum_plus_abs, double* length_to_boundary, double* v_length, double* time_to_boundary, - double* abs_weight_factor, int* abs_weight_factor_set) { - *my_sum_plus_abs = Volume->p_physics->my_a * (2200 / *v_length); - *length_to_boundary = *time_to_boundary * *v_length; - - *abs_weight_factor = exp (-Volume->p_physics->my_a * 2200 * *time_to_boundary); - *abs_weight_factor_set = 1; - - #ifdef Union_trace_verbal_setting - printf ("name of material: %s \n", Volume->name); - printf ("length to boundery = %f\n", length_to_boundary); - printf ("absorption cross section = %f\n", Volume->p_physics->my_a); - printf ("chance to get through this length of absorber: %f %%\n", 100 * exp (-Volume->p_physics->my_a * length_to_boundary)); - #endif - } - - - void - transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double (*k_rotated)[3], int* p_index, - Coords* wavevector, Coords* ray_position_geometry) { - - int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; - *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); - coords_get (*wavevector_rotated, &(*k_rotated)[0], &(*k_rotated)[1], &(*k_rotated)[2]); - - if (Volume->p_physics->p_scattering_array[*p_index].needs_cross_section_focus == 1) { - // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time - // focusing - Coords ray_position_geometry_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *ray_position_geometry); - - int focus_data_index = Volume->geometry.focus_array_indices.elements[*p_index]; - struct focus_data_struct* this_focus_data = &Volume->geometry.focus_data_array.elements[focus_data_index]; - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray - - #ifdef Union_trace_verbal_setting - printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", *p_index); - print_position (*ray_position_geometry, "ray_position_geometry"); - print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); - print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - #endif - } - } - - - - void - move_and_aim_neutron (struct physics_struct* current_p_physics, int i, struct scattering_process_struct* process, _class_particle* _particle, - struct Volume_struct* Volume, int* p_index, Coords* ray_velocity, Coords* ray_position, struct focus_data_struct* this_focus_data) { - // Transport neutron to place inside geometry - *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->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, Volume->geometry.center); - // Also focus the ray at this point, if the component needs focusing - if (process->needs_cross_section_focus) { - if (Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index != -1) { - int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; - sampling_point_geometry = rot_apply (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, &_particle->x, &_particle->y, &_particle->z); - } - - int - mu_and_intersect_dist_safeguard (double mu_sum, double length_to_boundary, double safety_distance2, int* scattering_event) { - if (mu_sum < 1E-18) { - *scattering_event = 0; - return 0; - } - if (length_to_boundary < safety_distance2) { - *scattering_event = 0; - return 0; - } - return 1; - } - - - void - focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double* safety_distance, double* safety_distance2, double* length_to_boundary, - _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, - struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { - // Sample length_to_scattering in linear manner - *forced_length_to_scattering = *safety_distance + rand01 () * (*length_to_boundary - *safety_distance2); - - *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); // Test for root cause - // 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 scattering_displacement = coords_scalar_mult (direction, *forced_length_to_scattering); - Coords forced_ray_scattering_point = coords_add (*ray_position, scattering_displacement); - *ray_position_geometry = coords_sub (forced_ray_scattering_point, Volume->geometry.center); // ray_position relative to geometry center - - // Calculate the aim for non isotropic processes - this_focus_data = &Volume->geometry.focus_data_array.elements[0]; - this_focus_data->RayAim = coords_sub (this_focus_data->Aim, *ray_position_geometry); // Aim vector for this ray - - #ifdef Union_trace_verbal_setting - printf ("Prepared for focus in cross section calculation in volume: %s \n", Volume->name); - printf ("forced_length_to_scattering =%lf \n", forced_length_to_scattering); - print_position (*ray_position, "ray_position"); - print_position (direction, "direction"); - print_position (scattering_displacement, "scattering_displacement"); - print_position (forced_ray_scattering_point, "forced_ray_scattering_point"); - print_position (*ray_position_geometry, "ray_position_geometry"); - printf ("for isotropic processes this RayAim is used \n"); - print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); - #endif - } - - void - focus_in_cross_section_set_scat_length (double* length_to_scattering, double* forced_length_to_scattering, double* weight, double* length_to_boundary, - double* my_sum_plus_abs) { - // 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)) - *weight *= *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 - } - - - void - inhomogenous_set_cumul_dist_array (struct physics_struct* current_p_physics, int 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; - } - - void - inhomogenous_sample_transmission_probability (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* real_transmission_probability, - double v_length) { - // 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 < Volume->p_physics->number_of_processes; i++) { - struct scattering_process_struct* process_i = &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 = 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]; - } - - double - inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* current_p_physics, double* abs_weight_factor, double* v_length, - double* safety_distance, int* selected_sampling) { - - // 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 = 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; - return current_p_physics->cumul_dists[*selected_sampling] - current_p_physics->dist / 2 + sampled_dist; - } - - - void - inhomogenous_choose_process (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* culmative_probability, double* mc_prop, - double* my_sum, int* selected_sampling, int* selected_process) { - for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { - *culmative_probability += current_p_physics->mus[i][*selected_sampling] / *my_sum; - if (*culmative_probability > *mc_prop) { - *selected_process = i; - break; - } - } - } - - - int - p_interact_is_set (struct Volume_struct* Volume) { - if (Volume->geometry.geometry_p_interact != 0) { - return 1; - } - return 0; - } - - void - p_interact_check_scattering_event (struct Volume_struct* Volume, int* scattering_event, double* weight, double real_transmission_prob) { - double mc_transmission_prob = 1 - Volume->geometry.geometry_p_interact; - *scattering_event = rand01 () > mc_transmission_prob; - if (*scattering_event) { - // Scattering event happens, this is the correction for the weight - *weight *= (1.0 - real_transmission_prob) / (1.0 - mc_transmission_prob); - } else { - // Scattering event does not happen, this is the appropriate correction - *weight *= real_transmission_prob / mc_transmission_prob; - } - } - - void - p_interact_select_process (struct Volume_struct* Volume, double* my_trace_fraction_control, double* my_trace, double* total_process_interact, - double* culmative_probability, double* mc_prop, double* weight, double* my_sum, int* selected_process) { - // Interact_fraction is used to influence the choice of process in this material - *mc_prop = rand01 (); - *culmative_probability = 0; - *total_process_interact = 1.0; - - // If any of the processes have probability 0, they are excluded from the selection - for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { - if (my_trace[i] < 1E-18) { - // When this happens, the total force probability is corrected and the probability for this particular instance is set to 0 - *total_process_interact -= Volume->p_physics->p_scattering_array[i].process_p_interact; - my_trace_fraction_control[i] = 0; - // In cases where my_trace is not zero, the forced fraction is still used. - } else - my_trace_fraction_control[i] = Volume->p_physics->p_scattering_array[i].process_p_interact; - } - // Randomly select a process using the weights stored in my_trace_fraction_control divided by total_process_interact - for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { - *culmative_probability += my_trace_fraction_control[i] / *total_process_interact; - if (*culmative_probability > *mc_prop) { - *selected_process = i; - *weight *= (my_trace[i] / *my_sum) * (*total_process_interact / my_trace_fraction_control[i]); - break; - } - } - } %} From 25f865961072d941b2136727fa3893e44588a8d7 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 16 Sep 2026 13:58:39 +0200 Subject: [PATCH 21/36] Use copying of values instead of pointer passing for doubles and ints that are not changed in the functions --- mcstas-comps/share/union-lib.c | 58 ++++++++++++++-------------- mcstas-comps/union/Union_master.comp | 18 ++++----- 2 files changed, 38 insertions(+), 38 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index dce01ed98b..de62bb61eb 100755 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -9228,7 +9228,6 @@ void overwrite_if_empty(char *input_string, char *overwrite) { volume_is_only_absorber (struct Volume_struct* Volume) { // This function returns true if a volume does not have any physical processes // and if the volume is not a vacuum. - printf("TESTING TESTING \n"); if (!Volume->p_physics->number_of_processes && !Volume->p_physics->is_vacuum) return 1; return 0; @@ -9244,12 +9243,13 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } void - adjust_abs_weight_factor (struct Volume_struct* Volume, double* my_sum_plus_abs, double* length_to_boundary, double* v_length, double* time_to_boundary, + adjust_abs_weight_factor (struct Volume_struct* Volume, double* my_sum_plus_abs, + double* length_to_boundary, double v_length, double time_to_boundary, double* abs_weight_factor, int* abs_weight_factor_set) { - *my_sum_plus_abs = Volume->p_physics->my_a * (2200 / *v_length); - *length_to_boundary = *time_to_boundary * *v_length; + *my_sum_plus_abs = Volume->p_physics->my_a * (2200 / v_length); + *length_to_boundary = time_to_boundary * v_length; - *abs_weight_factor = exp (-Volume->p_physics->my_a * 2200 * *time_to_boundary); + *abs_weight_factor = exp (-Volume->p_physics->my_a * 2200 * time_to_boundary); *abs_weight_factor_set = 1; #ifdef Union_trace_verbal_setting @@ -9262,24 +9262,24 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void - transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double (*k_rotated)[3], int* p_index, + transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double (*k_rotated)[3], int p_index, Coords* wavevector, Coords* ray_position_geometry) { - int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); coords_get (*wavevector_rotated, &(*k_rotated)[0], &(*k_rotated)[1], &(*k_rotated)[2]); - if (Volume->p_physics->p_scattering_array[*p_index].needs_cross_section_focus == 1) { + if (Volume->p_physics->p_scattering_array[p_index].needs_cross_section_focus == 1) { // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time // focusing Coords ray_position_geometry_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *ray_position_geometry); - int focus_data_index = Volume->geometry.focus_array_indices.elements[*p_index]; + int focus_data_index = Volume->geometry.focus_array_indices.elements[p_index]; struct focus_data_struct* this_focus_data = &Volume->geometry.focus_data_array.elements[focus_data_index]; this_focus_data->RayAim = coords_sub (this_focus_data->Aim, ray_position_geometry_rotated); // Aim vector for this ray #ifdef Union_trace_verbal_setting - printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", *p_index); + printf ("Checking process number : %d, it was not isotropic, so RayAim updated \n", p_index); print_position (*ray_position_geometry, "ray_position_geometry"); print_position (ray_position_geometry_rotated, "ray_position_geometry_rotated"); print_position (this_focus_data->RayAim, "this_focus_data->RayAim"); @@ -9291,7 +9291,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void move_and_aim_neutron (struct physics_struct* current_p_physics, int i, struct scattering_process_struct* process, _class_particle* _particle, - struct Volume_struct* Volume, int* p_index, Coords* ray_velocity, Coords* ray_position, struct focus_data_struct* this_focus_data) { + struct Volume_struct* Volume, int p_index, Coords* ray_velocity, Coords* ray_position, struct focus_data_struct* this_focus_data) { // Transport neutron to place inside geometry *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); // Find location of scattering point in master coordinate system without changing main position / velocity variables @@ -9301,8 +9301,8 @@ void overwrite_if_empty(char *input_string, char *overwrite) { Coords sampling_point_geometry = coords_sub (sampling_point, Volume->geometry.center); // Also focus the ray at this point, if the component needs focusing if (process->needs_cross_section_focus) { - if (Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index != -1) { - int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[*p_index].non_isotropic_rot_index; + if (Volume->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { + int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; sampling_point_geometry = rot_apply (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); @@ -9326,11 +9326,11 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void - focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double* safety_distance, double* safety_distance2, double* length_to_boundary, + focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double safety_distance, double safety_distance2, double length_to_boundary, _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { // Sample length_to_scattering in linear manner - *forced_length_to_scattering = *safety_distance + rand01 () * (*length_to_boundary - *safety_distance2); + *forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->vz); // Test for root cause // Find location of scattering point in master coordinate system without changing main position / velocity variables @@ -9357,14 +9357,14 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } void - focus_in_cross_section_set_scat_length (double* length_to_scattering, double* forced_length_to_scattering, double* weight, double* length_to_boundary, - double* my_sum_plus_abs) { + focus_in_cross_section_set_scat_length (double* length_to_scattering, double forced_length_to_scattering, double* weight, double length_to_boundary, + double my_sum_plus_abs) { // Respect forced length to scattering chosen by process - *length_to_scattering = *forced_length_to_scattering; + *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)) - *weight *= *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)); + *weight *= 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 @@ -9413,14 +9413,14 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } double - inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* current_p_physics, double* abs_weight_factor, double* v_length, - double* safety_distance, int* selected_sampling) { + inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* current_p_physics, double* abs_weight_factor, double v_length, + double safety_distance, int* selected_sampling) { // 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 = Volume->p_physics->my_a * (2200 / *v_length); + double mu_at_speed = 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); @@ -9435,7 +9435,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { // printf("\nSelected_sampling = %d\n", selected_sampling); // printf("dist i = %g\tdist=%g\n", dist_i, dist); - double sampled_dist = *safety_distance + 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; return current_p_physics->cumul_dists[*selected_sampling] - current_p_physics->dist / 2 + sampled_dist; @@ -9443,11 +9443,11 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void - inhomogenous_choose_process (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* culmative_probability, double* mc_prop, - double* my_sum, int* selected_sampling, int* selected_process) { + inhomogenous_choose_process (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* culmative_probability, double mc_prop, + double my_sum, int selected_sampling, int* selected_process) { for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { - *culmative_probability += current_p_physics->mus[i][*selected_sampling] / *my_sum; - if (*culmative_probability > *mc_prop) { + *culmative_probability += current_p_physics->mus[i][selected_sampling] / my_sum; + if (*culmative_probability > mc_prop) { *selected_process = i; break; } @@ -9478,7 +9478,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void p_interact_select_process (struct Volume_struct* Volume, double* my_trace_fraction_control, double* my_trace, double* total_process_interact, - double* culmative_probability, double* mc_prop, double* weight, double* my_sum, int* selected_process) { + double* culmative_probability, double* mc_prop, double* weight, double my_sum, int* selected_process) { // Interact_fraction is used to influence the choice of process in this material *mc_prop = rand01 (); *culmative_probability = 0; @@ -9499,7 +9499,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { *culmative_probability += my_trace_fraction_control[i] / *total_process_interact; if (*culmative_probability > *mc_prop) { *selected_process = i; - *weight *= (my_trace[i] / *my_sum) * (*total_process_interact / my_trace_fraction_control[i]); + *weight *= (my_trace[i] / my_sum) * (*total_process_interact / my_trace_fraction_control[i]); break; } } diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 1f930b697d..b79e73a92b 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -1547,7 +1547,7 @@ TRACE // Check if a scattering event should occur if (current_volume != 0) { // Volume 0 is always vacuum, and if this is the current volume, an event will not occur if (volume_is_only_absorber (Volumes[current_volume])) { // If there are no processes, the volume could be vacuum or an absorber - adjust_abs_weight_factor (Volumes[current_volume], &my_sum_plus_abs, &length_to_boundary, &v_length, &time_to_boundery, &abs_weight_factor, + adjust_abs_weight_factor (Volumes[current_volume], &my_sum_plus_abs, &length_to_boundary, v_length, time_to_boundery, &abs_weight_factor, &abs_weight_factor_set); } else { // Since there is a non-zero number of processes in this material, all the scattering cross section for these are calculated @@ -1564,7 +1564,7 @@ TRACE length_to_boundary = time_to_boundery * v_length; // If any process in this material needs focusing, sample scattering position and update focus_data accordingly if (current_p_physics->any_process_needs_cross_section_focus == 1) { - focus_in_cross_section_set_forced_point_and_dir (&forced_length_to_scattering, &safety_distance, &safety_distance2, &length_to_boundary, _particle, + focus_in_cross_section_set_forced_point_and_dir (&forced_length_to_scattering, safety_distance, safety_distance2, length_to_boundary, _particle, &ray_velocity, &ray_position_geometry, &ray_position, Volumes[current_volume], this_focus_data); } else { forced_length_to_scattering = -1.0; // Signals that no forcing needed, could also if on the selected process struct @@ -1579,7 +1579,7 @@ TRACE if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { // If the process is not isotropic, the wavevector is transformed into the local coordinate system of the process - transform_wavevector_into_local_coord_system (Volumes[current_volume], &wavevector_rotated, &k_rotated, &p_index, &wavevector, + transform_wavevector_into_local_coord_system (Volumes[current_volume], &wavevector_rotated, &k_rotated, p_index, &wavevector, &ray_position_geometry); } else { k_rotated[0] = k[0]; @@ -1603,7 +1603,7 @@ TRACE this_focus_data = &Volumes[current_volume]->geometry.focus_data_array.elements[0]; for (int i = 0; i < current_p_physics->sampling_points; i++) { inhomogenous_set_cumul_dist_array (current_p_physics, i); - move_and_aim_neutron (current_p_physics, i, process, _particle, Volumes[current_volume], &p_index, &ray_velocity, &ray_position, this_focus_data); + move_and_aim_neutron (current_p_physics, i, process, _particle, Volumes[current_volume], p_index, &ray_velocity, &ray_position, this_focus_data); // Calculate mu and probability physics_my (process->eProcess, &mu, k_rotated, process->data_transfer, this_focus_data, _particle); current_p_physics->mus[p_index][i] = mu; @@ -1669,8 +1669,8 @@ TRACE printf ("WARNING: Absorption weight factor above 1! Should not happen! \n"); // Select distance to scattering position if (current_p_physics->sampling_points != 0) { - length_to_scattering = inhomogenous_sample_scattering_point (Volumes[current_volume], current_p_physics, &abs_weight_factor, &v_length, - &safety_distance, &selected_sampling); + length_to_scattering = inhomogenous_sample_scattering_point (Volumes[current_volume], current_p_physics, &abs_weight_factor, v_length, + safety_distance, &selected_sampling); } // Select process if (Volumes[current_volume]->p_physics->number_of_processes == 1) { // trivial case @@ -1690,12 +1690,12 @@ TRACE } } } else { - inhomogenous_choose_process (current_p_physics, Volumes[current_volume], &culmative_probability, &mc_prop, &my_sum, &selected_sampling, + inhomogenous_choose_process (current_p_physics, Volumes[current_volume], &culmative_probability, mc_prop, my_sum, selected_sampling, &selected_process); } } else { p_interact_select_process (Volumes[current_volume], my_trace_fraction_control, my_trace, &total_process_interact, &culmative_probability, - &mc_prop, &p, &my_sum, &selected_process); + &mc_prop, &p, my_sum, &selected_process); } } @@ -1703,7 +1703,7 @@ TRACE if (current_p_physics->sampling_points == 0) { // No numerical integration is necessary. if (process->needs_cross_section_focus == 1) { - focus_in_cross_section_set_scat_length (&length_to_scattering, &forced_length_to_scattering, &p, &length_to_boundary, &my_sum_plus_abs); + focus_in_cross_section_set_scat_length (&length_to_scattering, forced_length_to_scattering, &p, length_to_boundary, my_sum_plus_abs); } else { // Decided the ray scatters, choose where on truncated exponential from safety_distance to length_to_boundary - safety_distance length_to_scattering From 8b2bb9b3364735e5560178c56aa9712124f7d06d Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 7 Oct 2026 15:38:54 +0200 Subject: [PATCH 22/36] Set initialization of sampling points to -1 both in physics struct and processes. Remove the barely used need_numerical_integration as it was only a pseudonym for sampling_points != -1 --- mcstas-comps/share/union-lib.c | 5 ++--- mcstas-comps/union/Union_master.comp | 15 +++++++-------- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 1ef130c34d..31be13763d 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -25,7 +25,6 @@ void scattering_process_struct_init(struct scattering_process_struct *sps) 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; } @@ -8271,7 +8270,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { int process_needs_inhomogenous_sampling (struct physics_struct* current_p_physics, struct scattering_process_struct* process) { - if (current_p_physics->sampling_points != 0) { + if (current_p_physics->sampling_points != -1) { if (process->needs_cross_section_focus || process->sampling_points != -1) return 1; } @@ -8420,7 +8419,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { struct scattering_process_struct* process_i = &Volume->p_physics->p_scattering_array[i]; - if (process_i->needs_numerical_integration != 1) + if (process_i->sampling_points != -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; } diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 54bc8e93aa..c1feb1836a 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -853,10 +853,9 @@ INITIALIZE 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; + physics->sampling_points = -1; // Initialize n points to avoid undefined behaviour 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; } @@ -1611,7 +1610,7 @@ TRACE } if (!process_needs_inhomogenous_sampling (current_p_physics, process)) { UNION_PHYSICS_MY (physics_output, process, p_my_trace, k_rotated, this_focus_data, _particle); - if (current_p_physics->sampling_points != 0) { + if (current_p_physics->sampling_points != -1) { current_p_physics->mus[p_index][0] = *p_my_trace; } } @@ -1638,11 +1637,11 @@ TRACE // Calculate if scattering happens based on my_sub_plus_abs if (mu_and_intersect_dist_safeguard (my_sum, length_to_boundary, safety_distance2, &scattering_event)) { // First calculate the transmission probability - if (current_p_physics->sampling_points == 0) { + if (current_p_physics->sampling_points == -1) { real_transmission_probability = exp (-length_to_boundary * my_sum_plus_abs); } - if (current_p_physics->sampling_points != 0) { + if (current_p_physics->sampling_points != -1) { inhomogenous_sample_transmission_probability (current_p_physics, Volumes[current_volume], &real_transmission_probability, v_length); } @@ -1666,7 +1665,7 @@ TRACE 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) { + if (current_p_physics->sampling_points != -1) { length_to_scattering = inhomogenous_sample_scattering_point (Volumes[current_volume], current_p_physics, &abs_weight_factor, v_length, safety_distance, &selected_sampling); } @@ -1679,7 +1678,7 @@ TRACE // Select a process based on their relative attenuations factors mc_prop = rand01 (); culmative_probability = 0; - if (current_p_physics->sampling_points == 0) { + if (current_p_physics->sampling_points == -1) { 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) { @@ -1698,7 +1697,7 @@ TRACE } process = &Volumes[current_volume]->p_physics->p_scattering_array[selected_process]; - if (current_p_physics->sampling_points == 0) { + if (current_p_physics->sampling_points == -1) { // No numerical integration is necessary. if (process->needs_cross_section_focus == 1) { focus_in_cross_section_set_scat_length (&length_to_scattering, forced_length_to_scattering, &p, length_to_boundary, my_sum_plus_abs); From c59743791253cf2589296d43c5f7045a36645b54 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 7 Oct 2026 16:00:13 +0200 Subject: [PATCH 23/36] Simplify the k_rotated in accordance with review --- mcstas-comps/share/union-lib.c | 4 ++-- mcstas-comps/share/union-lib.h | 2 +- mcstas-comps/union/Union_master.comp | 2 +- 3 files changed, 4 insertions(+), 4 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 31be13763d..835bca3940 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8297,12 +8297,12 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void - transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double (*k_rotated)[3], int p_index, + transform_wavevector_into_local_coord_system (struct Volume_struct* Volume, Coords* wavevector_rotated, double* k_rotated, int p_index, Coords* wavevector, Coords* ray_position_geometry) { int non_isotropic_rot_index = Volume->p_physics->p_scattering_array[p_index].non_isotropic_rot_index; *wavevector_rotated = rot_apply (Volume->geometry.process_rot_matrix_array[non_isotropic_rot_index], *wavevector); - coords_get (*wavevector_rotated, &(*k_rotated)[0], &(*k_rotated)[1], &(*k_rotated)[2]); + coords_get (*wavevector_rotated, &k_rotated[0], &k_rotated[1], &k_rotated[2]); if (Volume->p_physics->p_scattering_array[p_index].needs_cross_section_focus == 1) { // Prepare focus data using ray_position_geometry of forced scattering point which will be prepared if any process needs cross_section time diff --git a/mcstas-comps/share/union-lib.h b/mcstas-comps/share/union-lib.h index 6b0cbfd832..8043a16c01 100755 --- a/mcstas-comps/share/union-lib.h +++ b/mcstas-comps/share/union-lib.h @@ -1478,7 +1478,7 @@ void adjust_abs_weight_factor (struct Volume_struct *Volume, double *my_sum_plus double *abs_weight_factor, int *abs_weight_factor_set); void transform_wavevector_into_local_coord_system (struct Volume_struct *Volume, Coords *wavevector_rotated, - double (*k_rotated)[3], int p_index, Coords *wavevector, + double *k_rotated, int p_index, Coords *wavevector, Coords *ray_position_geometry); void move_and_aim_neutron (struct physics_struct *current_p_physics, int i, diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index c1feb1836a..3f3f189e44 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -1576,7 +1576,7 @@ TRACE if (Volumes[current_volume]->p_physics->p_scattering_array[p_index].non_isotropic_rot_index != -1) { // If the process is not isotropic, the wavevector is transformed into the local coordinate system of the process - transform_wavevector_into_local_coord_system (Volumes[current_volume], &wavevector_rotated, &k_rotated, p_index, &wavevector, + transform_wavevector_into_local_coord_system (Volumes[current_volume], &wavevector_rotated, k_rotated, p_index, &wavevector, &ray_position_geometry); } else { k_rotated[0] = k[0]; From 7340f63ab52c2409e3fa1129a8d82fe59483a1f3 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 7 Oct 2026 16:04:39 +0200 Subject: [PATCH 24/36] Rename mu_at_speed to mu_abs_at_speed --- mcstas-comps/share/union-lib.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 835bca3940..042640be5f 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8363,7 +8363,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double safety_distance, double safety_distance2, double length_to_boundary, _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, - struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { + struct Volume_struct* Volume, struct focus_data_struct* this_focus_data, ) { // Sample length_to_scattering in linear manner *forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); @@ -8431,9 +8431,9 @@ void overwrite_if_empty(char *input_string, char *overwrite) { // for (int i =0;isampling_points;i++){ // printf("\nTotalmu=%g\tinteger=%d\n", current_p_physics->total_mus[i], i); // } - double mu_at_speed = Volume->p_physics->my_a * (2200 / v_length); + double mu_abs_at_speed = 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; + current_p_physics->total_mus[j] += mu_abs_at_speed * current_p_physics->dist; } double trans_prob; for (int i = 0; i < current_p_physics->sampling_points; i++) { @@ -8455,7 +8455,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { // 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 = Volume->p_physics->my_a * (2200 / v_length); + double mu_abs_at_speed = 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); @@ -8465,7 +8465,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { 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]; + *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed * current_p_physics->dist) / current_p_physics->total_mus[*selected_sampling]; // printf("\nSelected_sampling = %d\n", selected_sampling); From 4bb4c3e5fa012cd31d3e4b635d299899b44a3bf5 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 7 Oct 2026 16:08:30 +0200 Subject: [PATCH 25/36] Pass _particle to inhomogeneous_sample_scattering point --- mcstas-comps/share/union-lib.c | 4 ++-- mcstas-comps/share/union-lib.h | 1 + mcstas-comps/union/Union_master.comp | 2 +- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 042640be5f..25f491cda1 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8363,7 +8363,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void focus_in_cross_section_set_forced_point_and_dir (double* forced_length_to_scattering, double safety_distance, double safety_distance2, double length_to_boundary, _class_particle* _particle, Coords* ray_velocity, Coords* ray_position_geometry, Coords* ray_position, - struct Volume_struct* Volume, struct focus_data_struct* this_focus_data, ) { + struct Volume_struct* Volume, struct focus_data_struct* this_focus_data) { // Sample length_to_scattering in linear manner *forced_length_to_scattering = safety_distance + rand01 () * (length_to_boundary - safety_distance2); @@ -8448,7 +8448,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } double - inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* current_p_physics, double* abs_weight_factor, double v_length, + inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* current_p_physics, _class_particle* _particle, double* abs_weight_factor, double v_length, double safety_distance, int* selected_sampling) { // Numerical integration happens, and therefore we must choose between the different samples diff --git a/mcstas-comps/share/union-lib.h b/mcstas-comps/share/union-lib.h index 8043a16c01..7e6053d535 100755 --- a/mcstas-comps/share/union-lib.h +++ b/mcstas-comps/share/union-lib.h @@ -1507,6 +1507,7 @@ void inhomogenous_sample_transmission_probability (struct physics_struct *curren double inhomogenous_sample_scattering_point (struct Volume_struct *Volume, struct physics_struct *current_p_physics, + _class_particle* _particle, double *abs_weight_factor, double v_length, double safety_distance, int *selected_sampling); diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 3f3f189e44..4c5e678ee8 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -1666,7 +1666,7 @@ TRACE printf ("WARNING: Absorption weight factor above 1! Should not happen! \n"); // Select distance to scattering position if (current_p_physics->sampling_points != -1) { - length_to_scattering = inhomogenous_sample_scattering_point (Volumes[current_volume], current_p_physics, &abs_weight_factor, v_length, + length_to_scattering = inhomogenous_sample_scattering_point (Volumes[current_volume], current_p_physics, _particle, &abs_weight_factor, v_length, safety_distance, &selected_sampling); } // Select process From 373b056a5238d17862c119af763c91a17b5cb2ab Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 7 Oct 2026 16:11:53 +0200 Subject: [PATCH 26/36] invert logic in inhomogeneous position sampling for loop --- mcstas-comps/share/union-lib.c | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 25f491cda1..6da8eed704 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8459,10 +8459,9 @@ void overwrite_if_empty(char *input_string, char *overwrite) { 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; + if (pseudo_rand < 1 - current_p_physics->cumul_transmission_prob[i]) + *selected_sampling = i; + break; } *abs_weight_factor *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed * current_p_physics->dist) / current_p_physics->total_mus[*selected_sampling]; From 261df9c7467743600875ff536eba929465dee253 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 7 Oct 2026 16:14:48 +0200 Subject: [PATCH 27/36] Remove erroneous multiplication by dist in weight correction for sampling scattering in inhomogeneous logic path --- mcstas-comps/share/union-lib.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 6da8eed704..b5f6a21ed9 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8464,7 +8464,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { break; } *abs_weight_factor - *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed * current_p_physics->dist) / current_p_physics->total_mus[*selected_sampling]; + *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed / current_p_physics->total_mus[*selected_sampling]; // printf("\nSelected_sampling = %d\n", selected_sampling); From a89b3678d9d60d6c95dc08cbb7fce9cad6512940 Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 7 Oct 2026 16:19:04 +0200 Subject: [PATCH 28/36] Pass the particle to the p_interact logic paths, and fix a missing parenthesis in weight correction for inhomogenous processes --- mcstas-comps/share/union-lib.c | 6 +++--- mcstas-comps/share/union-lib.h | 4 ++-- mcstas-comps/union/Union_master.comp | 4 ++-- 3 files changed, 7 insertions(+), 7 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index b5f6a21ed9..f671a1ff2e 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8464,7 +8464,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { break; } *abs_weight_factor - *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed / current_p_physics->total_mus[*selected_sampling]; + *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed ) / current_p_physics->total_mus[*selected_sampling]; // printf("\nSelected_sampling = %d\n", selected_sampling); @@ -8498,7 +8498,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } void - p_interact_check_scattering_event (struct Volume_struct* Volume, int* scattering_event, double* weight, double real_transmission_prob) { + p_interact_check_scattering_event (struct Volume_struct* Volume, _class_particle* _particle, int* scattering_event, double* weight, double real_transmission_prob) { double mc_transmission_prob = 1 - Volume->geometry.geometry_p_interact; *scattering_event = rand01 () > mc_transmission_prob; if (*scattering_event) { @@ -8511,7 +8511,7 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } void - p_interact_select_process (struct Volume_struct* Volume, double* my_trace_fraction_control, double* my_trace, double* total_process_interact, + p_interact_select_process (struct Volume_struct* Volume, _class_particle* _particle, double* my_trace_fraction_control, double* my_trace, double* total_process_interact, double* culmative_probability, double* mc_prop, double* weight, double my_sum, int* selected_process) { // Interact_fraction is used to influence the choice of process in this material *mc_prop = rand01 (); diff --git a/mcstas-comps/share/union-lib.h b/mcstas-comps/share/union-lib.h index 7e6053d535..420d3dd4ea 100755 --- a/mcstas-comps/share/union-lib.h +++ b/mcstas-comps/share/union-lib.h @@ -1517,10 +1517,10 @@ void inhomogenous_choose_process (struct physics_struct *current_p_physics, stru int p_interact_is_set (struct Volume_struct *Volume); -void p_interact_check_scattering_event (struct Volume_struct *Volume, int *scattering_event, +void p_interact_check_scattering_event (struct Volume_struct *Volume, _class_particle* _particle, int *scattering_event, double *weight, double real_transmission_prob); -void p_interact_select_process (struct Volume_struct *Volume, double *my_trace_fraction_control, +void p_interact_select_process (struct Volume_struct *Volume, _class_particle* _particle, double *my_trace_fraction_control, double *my_trace, double *total_process_interact, double *culmative_probability, double *mc_prop, double *weight, double my_sum, int *selected_process); diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index 4c5e678ee8..f6db6a9b3d 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -1647,7 +1647,7 @@ TRACE // Then check if we scatter. if (p_interact_is_set (Volumes[current_volume])) { - p_interact_check_scattering_event (Volumes[current_volume], &scattering_event, &p, real_transmission_probability); + p_interact_check_scattering_event (Volumes[current_volume], _particle, &scattering_event, &p, real_transmission_probability); } else { // probability to scatter is the natural value scattering_event = rand01 () > real_transmission_probability; @@ -1691,7 +1691,7 @@ TRACE &selected_process); } } else { - p_interact_select_process (Volumes[current_volume], my_trace_fraction_control, my_trace, &total_process_interact, &culmative_probability, + p_interact_select_process (Volumes[current_volume], _particle, my_trace_fraction_control, my_trace, &total_process_interact, &culmative_probability, &mc_prop, &p, my_sum, &selected_process); } } From f2c073c7d5d0b9c2014210e0256415f13fe99aae Mon Sep 17 00:00:00 2001 From: Diablo Date: Wed, 7 Oct 2026 16:22:59 +0200 Subject: [PATCH 29/36] Remove unnecessary if statement, add comment instead --- mcstas-comps/union/Union_master.comp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index f6db6a9b3d..df8c07f57b 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -1543,7 +1543,7 @@ TRACE // Check if a scattering event should occur if (current_volume != 0) { // Volume 0 is always vacuum, and if this is the current volume, an event will not occur - if (volume_is_only_absorber (Volumes[current_volume])) { // If there are no processes, the volume could be vacuum or an absorber + if (volume_is_only_absorber (Volumes[current_volume])) { // If there are no processes, check if the volume is a pure absorber and not a vacuum adjust_abs_weight_factor (Volumes[current_volume], &my_sum_plus_abs, &length_to_boundary, v_length, time_to_boundery, &abs_weight_factor, &abs_weight_factor_set); } else { @@ -1608,7 +1608,7 @@ TRACE } coords_get (original_position, &x, &y, &z); } - if (!process_needs_inhomogenous_sampling (current_p_physics, process)) { + else {// process does not need inhomogenous sampling UNION_PHYSICS_MY (physics_output, process, p_my_trace, k_rotated, this_focus_data, _particle); if (current_p_physics->sampling_points != -1) { current_p_physics->mus[p_index][0] = *p_my_trace; From f11f560f44060ee1dc5c941af7a75ec88755e67f Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 8 Oct 2026 14:42:28 +0200 Subject: [PATCH 30/36] Remove needs_numerical integration from union-lib.h since the parameter is no longer in use --- mcstas-comps/share/union-lib.h | 1 - 1 file changed, 1 deletion(-) diff --git a/mcstas-comps/share/union-lib.h b/mcstas-comps/share/union-lib.h index 420d3dd4ea..cd56efa1b2 100755 --- a/mcstas-comps/share/union-lib.h +++ b/mcstas-comps/share/union-lib.h @@ -637,7 +637,6 @@ struct scattering_process_struct 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. From abc9457e2ad74d4bdd9c0ccfc03b022e46e4ea48 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 8 Oct 2026 14:42:52 +0200 Subject: [PATCH 31/36] Add missing brackets to if statement with break included --- mcstas-comps/share/union-lib.c | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index f671a1ff2e..325de9ea76 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8458,17 +8458,14 @@ void overwrite_if_empty(char *input_string, char *overwrite) { double mu_abs_at_speed = 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]) + if (pseudo_rand < 1 - current_p_physics->cumul_transmission_prob[i]){ *selected_sampling = i; break; + } } *abs_weight_factor *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed ) / 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; From 1f4629bbadab5b03cc34952af63f88c4c169458a Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 8 Oct 2026 14:43:22 +0200 Subject: [PATCH 32/36] Remove needs_numerical_integration as it is no longer in use --- mcstas-comps/union/Inhomogenous_incoherent_process.comp | 1 - 1 file changed, 1 deletion(-) diff --git a/mcstas-comps/union/Inhomogenous_incoherent_process.comp b/mcstas-comps/union/Inhomogenous_incoherent_process.comp index a796d94f3c..bdab502635 100755 --- a/mcstas-comps/union/Inhomogenous_incoherent_process.comp +++ b/mcstas-comps/union/Inhomogenous_incoherent_process.comp @@ -369,7 +369,6 @@ INITIALIZE 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. From 5afae3e9f388135a3b75e6bde13828040ce2db32 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 8 Oct 2026 15:18:31 +0200 Subject: [PATCH 33/36] total_mus was implicitly using the cross section times the distance, thereby being misnamed. Fixed now --- mcstas-comps/share/union-lib.c | 15 ++++++--------- 1 file changed, 6 insertions(+), 9 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 325de9ea76..e6d05d6601 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8421,23 +8421,20 @@ void overwrite_if_empty(char *input_string, char *overwrite) { struct scattering_process_struct* process_i = &Volume->p_physics->p_scattering_array[i]; if (process_i->sampling_points != -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; + current_p_physics->total_mus[j] += current_p_physics->mus[i][0]; } 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; + current_p_physics->total_mus[j] += current_p_physics->mus[i][j]; } } - // for (int i =0;isampling_points;i++){ - // printf("\nTotalmu=%g\tinteger=%d\n", current_p_physics->total_mus[i], i); - // } double mu_abs_at_speed = 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_abs_at_speed * current_p_physics->dist; + current_p_physics->total_mus[j] += mu_abs_at_speed; } double trans_prob; for (int i = 0; i < current_p_physics->sampling_points; i++) { - trans_prob = exp (-current_p_physics->total_mus[i]); + trans_prob = exp (-current_p_physics->total_mus[i] * current_p_physics->dist); if (i == 0) current_p_physics->cumul_transmission_prob[i] = trans_prob; else @@ -8467,8 +8464,8 @@ void overwrite_if_empty(char *input_string, char *overwrite) { *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed ) / current_p_physics->total_mus[*selected_sampling]; 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; + - log (1.0 - rand01 () * (1.0 - exp (-current_p_physics->total_mus[*selected_sampling] * current_p_physics->dist))) + / current_p_physics->total_mus[*selected_sampling]; return current_p_physics->cumul_dists[*selected_sampling] - current_p_physics->dist / 2 + sampled_dist; } From 33149a55e6ff49ec05f8f1e52c3fd9c33eae79dc Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 8 Oct 2026 16:09:24 +0200 Subject: [PATCH 34/36] Remove unnecessary if statement --- mcstas-comps/union/Union_master.comp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/mcstas-comps/union/Union_master.comp b/mcstas-comps/union/Union_master.comp index df8c07f57b..82192f5e90 100755 --- a/mcstas-comps/union/Union_master.comp +++ b/mcstas-comps/union/Union_master.comp @@ -1639,9 +1639,7 @@ TRACE // First calculate the transmission probability if (current_p_physics->sampling_points == -1) { real_transmission_probability = exp (-length_to_boundary * my_sum_plus_abs); - } - - if (current_p_physics->sampling_points != -1) { + } else { inhomogenous_sample_transmission_probability (current_p_physics, Volumes[current_volume], &real_transmission_probability, v_length); } From 5bc7408d3574a9ceb0008cf3a423ccb11ce74f20 Mon Sep 17 00:00:00 2001 From: Diablo Date: Thu, 8 Oct 2026 16:11:26 +0200 Subject: [PATCH 35/36] Remove unnecessary current_p in front of physics as it is inside a function where the context makes it obvious --- mcstas-comps/share/union-lib.c | 58 +++++++++++++++++----------------- 1 file changed, 29 insertions(+), 29 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index e6d05d6601..9267a6a856 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8269,8 +8269,8 @@ void overwrite_if_empty(char *input_string, char *overwrite) { } int - process_needs_inhomogenous_sampling (struct physics_struct* current_p_physics, struct scattering_process_struct* process) { - if (current_p_physics->sampling_points != -1) { + process_needs_inhomogenous_sampling (struct physics_struct* physics, struct scattering_process_struct* process) { + if (physics->sampling_points != -1) { if (process->needs_cross_section_focus || process->sampling_points != -1) return 1; } @@ -8325,13 +8325,13 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void - move_and_aim_neutron (struct physics_struct* current_p_physics, int i, struct scattering_process_struct* process, _class_particle* _particle, + move_and_aim_neutron (struct physics_struct* physics, int i, struct scattering_process_struct* process, _class_particle* _particle, struct Volume_struct* Volume, int p_index, Coords* ray_velocity, Coords* ray_position, struct focus_data_struct* this_focus_data) { // Transport neutron to place inside geometry *ray_velocity = coords_set (_particle->vx, _particle->vy, _particle->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_displacement = coords_scalar_mult (direction, physics->cumul_dists[i]); Coords sampling_point = coords_add (*ray_position, sampling_displacement); Coords sampling_point_geometry = coords_sub (sampling_point, Volume->geometry.center); // Also focus the ray at this point, if the component needs focusing @@ -8407,45 +8407,45 @@ void overwrite_if_empty(char *input_string, char *overwrite) { void - inhomogenous_set_cumul_dist_array (struct physics_struct* current_p_physics, int 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; + inhomogenous_set_cumul_dist_array (struct physics_struct* physics, int i) { + physics->cumul_dists[i] = (i > 0) ? physics->cumul_dists[i - 1] + physics->dist : physics->dist / 2; } void - inhomogenous_sample_transmission_probability (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* real_transmission_probability, + inhomogenous_sample_transmission_probability (struct physics_struct* physics, struct Volume_struct* Volume, double* real_transmission_probability, double v_length) { // Calculate the probabilities and then add them cumulatively - memset (current_p_physics->total_mus, 0, sizeof (double) * current_p_physics->sampling_points); + memset (physics->total_mus, 0, sizeof (double) * physics->sampling_points); for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { struct scattering_process_struct* process_i = &Volume->p_physics->p_scattering_array[i]; if (process_i->sampling_points != -1) - for (int j = 0; j < current_p_physics->sampling_points; j++) { - current_p_physics->total_mus[j] += current_p_physics->mus[i][0]; + for (int j = 0; j < physics->sampling_points; j++) { + physics->total_mus[j] += physics->mus[i][0]; } else - for (int j = 0; j < current_p_physics->sampling_points; j++) { - current_p_physics->total_mus[j] += current_p_physics->mus[i][j]; + for (int j = 0; j < physics->sampling_points; j++) { + physics->total_mus[j] += physics->mus[i][j]; } } double mu_abs_at_speed = 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_abs_at_speed; + for (int j = 0; j < physics->sampling_points; j++) { + physics->total_mus[j] += mu_abs_at_speed; } double trans_prob; - for (int i = 0; i < current_p_physics->sampling_points; i++) { - trans_prob = exp (-current_p_physics->total_mus[i] * current_p_physics->dist); + for (int i = 0; i < physics->sampling_points; i++) { + trans_prob = exp (-physics->total_mus[i] * physics->dist); if (i == 0) - current_p_physics->cumul_transmission_prob[i] = trans_prob; + 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; + physics->cumul_transmission_prob[i] = physics->cumul_transmission_prob[i - 1] * trans_prob; } - *real_transmission_probability = current_p_physics->cumul_transmission_prob[current_p_physics->sampling_points - 1]; + *real_transmission_probability = physics->cumul_transmission_prob[physics->sampling_points - 1]; } double - inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* current_p_physics, _class_particle* _particle, double* abs_weight_factor, double v_length, + inhomogenous_sample_scattering_point (struct Volume_struct* Volume, struct physics_struct* physics, _class_particle* _particle, double* abs_weight_factor, double v_length, double safety_distance, int* selected_sampling) { // Numerical integration happens, and therefore we must choose between the different samples @@ -8453,28 +8453,28 @@ void overwrite_if_empty(char *input_string, char *overwrite) { // and then seeing which cumul prob is the first to include it. *abs_weight_factor = 1; double mu_abs_at_speed = 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++) { - if (pseudo_rand < 1 - current_p_physics->cumul_transmission_prob[i]){ + double pseudo_rand = rand01 () * (1 - physics->cumul_transmission_prob[physics->sampling_points - 1]); + for (int i = 0; i < physics->sampling_points; i++) { + if (pseudo_rand < 1 - physics->cumul_transmission_prob[i]){ *selected_sampling = i; break; } } *abs_weight_factor - *= (current_p_physics->total_mus[*selected_sampling] - mu_abs_at_speed ) / current_p_physics->total_mus[*selected_sampling]; + *= (physics->total_mus[*selected_sampling] - mu_abs_at_speed ) / physics->total_mus[*selected_sampling]; double sampled_dist = safety_distance - - log (1.0 - rand01 () * (1.0 - exp (-current_p_physics->total_mus[*selected_sampling] * current_p_physics->dist))) - / current_p_physics->total_mus[*selected_sampling]; - return current_p_physics->cumul_dists[*selected_sampling] - current_p_physics->dist / 2 + sampled_dist; + - log (1.0 - rand01 () * (1.0 - exp (-physics->total_mus[*selected_sampling] * physics->dist))) + / physics->total_mus[*selected_sampling]; + return physics->cumul_dists[*selected_sampling] - physics->dist / 2 + sampled_dist; } void - inhomogenous_choose_process (struct physics_struct* current_p_physics, struct Volume_struct* Volume, double* culmative_probability, double mc_prop, + inhomogenous_choose_process (struct physics_struct* physics, struct Volume_struct* Volume, double* culmative_probability, double mc_prop, double my_sum, int selected_sampling, int* selected_process) { for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { - *culmative_probability += current_p_physics->mus[i][selected_sampling] / my_sum; + *culmative_probability += physics->mus[i][selected_sampling] / my_sum; if (*culmative_probability > mc_prop) { *selected_process = i; break; From 6d4606a64cea4ea1994410770a12b68822ec76e4 Mon Sep 17 00:00:00 2001 From: Diablo Date: Fri, 9 Oct 2026 07:57:04 +0200 Subject: [PATCH 36/36] Flip process inhomogeneous condition in probability function --- mcstas-comps/share/union-lib.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/mcstas-comps/share/union-lib.c b/mcstas-comps/share/union-lib.c index 9267a6a856..2ff1696fc3 100644 --- a/mcstas-comps/share/union-lib.c +++ b/mcstas-comps/share/union-lib.c @@ -8419,14 +8419,15 @@ void overwrite_if_empty(char *input_string, char *overwrite) { for (int i = 0; i < Volume->p_physics->number_of_processes; i++) { struct scattering_process_struct* process_i = &Volume->p_physics->p_scattering_array[i]; - if (process_i->sampling_points != -1) + if (process_i->sampling_points == -1){ for (int j = 0; j < physics->sampling_points; j++) { physics->total_mus[j] += physics->mus[i][0]; } - else + } else { for (int j = 0; j < physics->sampling_points; j++) { physics->total_mus[j] += physics->mus[i][j]; } + } } double mu_abs_at_speed = Volume->p_physics->my_a * (2200 / v_length); for (int j = 0; j < physics->sampling_points; j++) {