Increase maximum number of allowed scatterings to 10000.

Particles that scatter more often are considered as "captured" in the statistics.
temken · Mar 31, 2023 · 07c9df0 · 07c9df0
1 parent e226aa4
commit 07c9df0
Show file tree

Hide file tree

Showing 3 changed files with 9 additions and 4 deletions.
diff --git a/include/Data_Generation.hpp b/include/Data_Generation.hpp
@@ -21,8 +21,8 @@ class Simulation_Data
 	unsigned int isoreflection_rings;
 	double initial_and_final_radius			   = 1.1 * libphysica::natural_units::rSun;
 	unsigned int minimum_number_of_scatterings = 1;
-	unsigned int maximum_number_of_scatterings = 1000;
-	unsigned long int maximum_free_time_steps  = 1e7;
+	unsigned int maximum_number_of_scatterings = 10000;
+	unsigned long int maximum_free_time_steps  = 1e8;
 
 	// Results
 	unsigned long int number_of_trajectories, number_of_free_particles, number_of_reflected_particles, number_of_captured_particles;
@@ -41,7 +41,7 @@ class Simulation_Data
 
 	Simulation_Data(unsigned int sample_size, double u_min = 0.0, unsigned int iso_rings = 1);
 
-	void Configure(double initial_radius, unsigned int min_scattering, unsigned int max_scattering, unsigned long int max_free_steps = 1e8);
+	void Configure(double initial_radius, unsigned int min_scattering, unsigned int max_scattering = 10000, unsigned long int max_free_steps = 1e8);
 
 	void Generate_Data(obscura::DM_Particle& DM, Solar_Model& solar_model, obscura::DM_Distribution& halo_model, unsigned int fixed_seed = 0);
 

diff --git a/src/Simulation_Trajectory.cpp b/src/Simulation_Trajectory.cpp
@@ -104,6 +104,12 @@ bool Trajectory_Simulator::Propagate_Freely(Event& current_event, obscura::DM_Pa
 					  << "\tAbort simulation." << std::endl;
 			return false;
 		}
+		else if(time_steps == maximum_time_steps)
+		{
+			std::cerr << "\nWarning in Propagate_Freely(): Number of time steps exceeds the maximum = " << maximum_time_steps << std::endl
+					  << "\tAbort simulation." << std::endl;
+			return false;
+		}
 
 		if(save_trajectories && time_steps % 20 == 0)
 		{

diff --git a/src/main.cpp b/src/main.cpp
@@ -52,7 +52,6 @@ int main(int argc, char* argv[])
 		double u_min = 0.0;
 		// double u_min = cfg.DM_detector->Minimum_DM_Speed(*cfg.DM);
 		Simulation_Data data_set(cfg.sample_size, u_min, cfg.isoreflection_rings);
-		data_set.Configure(1.1 * rSun, 1, 1000);
 		if(mpi_rank == 0)
 			std::cout << "\nDM parameters:" << std::endl
 					  << "\tm_DM [MeV]:\t" << libphysica::Round(In_Units(cfg.DM->mass, MeV)) << std::endl