Union refactoring of scattering while loop

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;
 }
 
@@ -8257,6 +8256,283 @@ 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.
+    if (!Volume->p_physics->number_of_processes && !Volume->p_physics->is_vacuum)
+      return 1;
+    return 0;
+  }
+
+  int
+  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;
+    }
+    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, 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* 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, 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* 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* physics, struct Volume_struct* Volume, double* real_transmission_probability,
+                                                double v_length) {
+    // Calculate the probabilities and then add them cumulatively
+    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 < physics->sampling_points; j++) {
+          physics->total_mus[j] += physics->mus[i][0];
+        }
+      } 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++) {
+      physics->total_mus[j] += mu_abs_at_speed;
+    }
+    double trans_prob;
+    for (int i = 0; i < physics->sampling_points; i++) {
+      trans_prob = exp (-physics->total_mus[i] * physics->dist);
+      if (i == 0)
+        physics->cumul_transmission_prob[i] = trans_prob;
+      else
+        physics->cumul_transmission_prob[i] = physics->cumul_transmission_prob[i - 1] * trans_prob;
+    }
+
+    *real_transmission_probability = physics->cumul_transmission_prob[physics->sampling_points - 1];
+  }
+
+  double
+  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
+    // 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_abs_at_speed = Volume->p_physics->my_a * (2200 / v_length);
+    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
+        *= (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 (-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* 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 += 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, _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) {
+      // 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, _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 ();
+    *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;
+      }
+    }
+  }
 struct union_state_struct *union_acquire_state(int master_present) {
   (void)master_present;
   union_state.refcount++;