Skip to content

Commit bd1c0c9

Browse files
JulienDoernerJulienDoerner
authored
Merge branch 'master' into ReflectiveShell
2 parents a9ab582 + 795ff09 commit bd1c0c9

File tree

13 files changed

+852
-96
lines changed

13 files changed

+852
-96
lines changed

‎.github/workflows/testing_OSX.yml

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -19,8 +19,8 @@ jobs:
1919
uses: actions/checkout@v4
2020
- name: Preinstall
2121
run: |
22-
brew install hdf5 fftw cfitsio muparser libomp swig
23-
pip install numpy==1.26
22+
brew install hdf5 fftw cfitsio muparser libomp swig python3
23+
pip install numpy==1.26.4 cmake
2424
- name: Set up the build
2525
env:
2626
CXX: ${{ matrix.config.cxx }}

‎CHANGELOG.md

Lines changed: 8 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -1,15 +1,20 @@
11
## CRPropa vNext
22

33
### Bug fixes:
4-
* Fixed sign for exponential decay of magn. field strength with Galactic height in LogarithmicSpiralField
5-
* Fixed r term in source distribution for SNR and Pulsar
6-
* Fixed wrong mass inheritance for secondaries other than nuclei or electron/positron
4+
* Fixed sign for exponential decay of magn. field strength with Galactic height in LogarithmicSpiralField
5+
* Fixed r term in source distribution for SNR and Pulsar
6+
* Fixed wrong mass inheritance for secondaries other than nuclei or electron/positron
7+
* Fixed wrong generation of interval ranges in ObserverTimeEvolution
78

89
### New features:
10+
911
* Added new backwards-compatible function particleMass that returns particle mass also for non-nuclei
1012
* Added the new Galactic magnetic field models from Unger&Farrar arXiv:2311.12120
1113
* Added EBL model from Finke et al. 2022
1214
* Added ReflectiveShell boundary condition
15+
* Added overridable getTime function to ObserverTimeEvolution which is called instead of detList.
16+
Old functionalities are preserved by adding old functions and setter/getter functions.
17+
1318

1419
### Interface changes:
1520

‎CMakeLists.txt

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -402,6 +402,7 @@ add_library(crpropa SHARED
402402
src/magneticField/turbulentField/SimpleGridTurbulence.cpp
403403
src/magneticField/TF17Field.cpp
404404
src/magneticField/UF23Field.cpp
405+
src/magneticField/KST24Field.cpp
405406
src/magneticField/CMZField.cpp
406407
src/advectionField/AdvectionField.cpp
407408
src/massDistribution/ConstantDensity.cpp

‎doc/pages/example_notebooks/Diffusion/MomentumDiffusion.ipynb

Lines changed: 29 additions & 42 deletions
Large diffs are not rendered by default.

‎include/CRPropa.h

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -65,6 +65,7 @@
6565
#include "crpropa/magneticField/QuimbyMagneticField.h"
6666
#include "crpropa/magneticField/TF17Field.h"
6767
#include "crpropa/magneticField/UF23Field.h"
68+
#include "crpropa/magneticField/KST24Field.h"
6869
#include "crpropa/magneticField/CMZField.h"
6970
#include "crpropa/magneticField/turbulentField/GridTurbulence.h"
7071
#include "crpropa/magneticField/turbulentField/HelicalGridTurbulence.h"

‎include/crpropa/ModuleList.h

Lines changed: 8 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,13 +1,19 @@
11
#ifndef CRPROPA_MODULE_LIST_H
22
#define CRPROPA_MODULE_LIST_H
33

4+
#include <algorithm>
5+
#include <csignal>
6+
#include <iostream>
7+
#include <vector>
8+
#include <exception>
9+
#include <sstream>
10+
#include <list>
11+
412
#include "crpropa/Candidate.h"
513
#include "crpropa/Module.h"
614
#include "crpropa/Source.h"
715
#include "crpropa/module/Output.h"
816

9-
#include <list>
10-
#include <sstream>
1117

1218
namespace crpropa {
1319

‎include/crpropa/magneticField/KST24Field.h

Lines changed: 154 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,154 @@
1+
#ifndef _KST24_GMF_H_
2+
#define _KST24_GMF_H_
3+
4+
#include <vector>
5+
#include "crpropa/magneticField/MagneticField.h"
6+
7+
namespace crpropa {
8+
9+
/*
10+
The C++ implementation of the GMF model KST24 (A.Korochkin, D.Semikoz, P.Tinyakov 2024)
11+
The model was presented in arXiv:2407.02148 and published in A&A
12+
If you use the model, please cite A&A, 693, A284 (2025)
13+
14+
In KST24 GMF model the position of the Solar System is at {-8.2 kpc, 0, 0}
15+
The Galactic Center is at {0, 0, 0}
16+
z-axis points to the North pole
17+
*/
18+
19+
20+
class KST24Field : public MagneticField
21+
{
22+
private:
23+
double north_tor_B_gauss;
24+
double north_tor_zmin_kpc;
25+
double north_tor_zmax_kpc;
26+
double north_tor_rmin_kpc;
27+
double north_tor_rmax_kpc;
28+
29+
double south_tor_B_gauss;
30+
double south_tor_zmin_kpc;
31+
double south_tor_zmax_kpc;
32+
double south_tor_rmin_kpc;
33+
double south_tor_rmax_kpc;
34+
35+
double Xfield_B_gauss;
36+
double Xfield_rmin_kpc;
37+
double Xfield_rmax_kpc;
38+
double Xfield_theta_deg;
39+
double Xfield_theta_rad;
40+
41+
double LB_B_gauss;
42+
double LB_lB_deg;
43+
double LB_bB_deg;
44+
double LB_rmin_kpc;
45+
double LB_dr_kpc;
46+
double LB_x0_kpc;
47+
double LB_y0_kpc;
48+
double LB_z0_kpc;
49+
double LB_Bdir_x;
50+
double LB_Bdir_y;
51+
double LB_Bdir_z;
52+
53+
double ScutumArm_B_gauss;
54+
double ScutumArm_pitch_deg;
55+
double ScutumArm_phi0_deg;
56+
double ScutumArm_x_shift_kpc;
57+
double ScutumArm_y_shift_kpc;
58+
double ScutumArm_arc_radius1_kpc;
59+
double ScutumArm_arc_radius2_kpc;
60+
double ScutumArm_arc_eps;
61+
double ScutumArm_arc_div_deg;
62+
double ScutumArm_rmin_kpc;
63+
double ScutumArm_rmax_kpc;
64+
double ScutumArm_zmin_kpc;
65+
double ScutumArm_zmax_kpc;
66+
67+
double CarinaSagittariusArm_B_gauss;
68+
double CarinaSagittariusArm_pitch_deg;
69+
double CarinaSagittariusArm_phi0_deg;
70+
double CarinaSagittariusArm_x_shift_kpc;
71+
double CarinaSagittariusArm_y_shift_kpc;
72+
double CarinaSagittariusArm_arc_radius1_kpc;
73+
double CarinaSagittariusArm_arc_radius2_kpc;
74+
double CarinaSagittariusArm_arc_eps;
75+
double CarinaSagittariusArm_arc_div_deg;
76+
double CarinaSagittariusArm_rmin_kpc;
77+
double CarinaSagittariusArm_rmax_kpc;
78+
double CarinaSagittariusArm_zmin_kpc;
79+
double CarinaSagittariusArm_zmax_kpc;
80+
81+
double LocalArm_B_gauss;
82+
double LocalArm_pitch_deg;
83+
double LocalArm_phi0_deg;
84+
double LocalArm_x_shift_kpc;
85+
double LocalArm_y_shift_kpc;
86+
double LocalArm_arc_radius1_kpc;
87+
double LocalArm_arc_radius2_kpc;
88+
double LocalArm_arc_eps;
89+
double LocalArm_arc_div_deg;
90+
double LocalArm_rmin_kpc;
91+
double LocalArm_rmax_kpc;
92+
double LocalArm_zmin_kpc;
93+
double LocalArm_zmax_kpc;
94+
95+
double PerseusArm1_B_gauss;
96+
double PerseusArm1_pitch_deg;
97+
double PerseusArm1_phi0_deg;
98+
double PerseusArm1_x_shift_kpc;
99+
double PerseusArm1_y_shift_kpc;
100+
double PerseusArm1_arc_radius1_kpc;
101+
double PerseusArm1_arc_radius2_kpc;
102+
double PerseusArm1_arc_eps;
103+
double PerseusArm1_arc_div_deg;
104+
double PerseusArm1_rmin_kpc;
105+
double PerseusArm1_rmax_kpc;
106+
double PerseusArm1_zmin_kpc;
107+
double PerseusArm1_zmax_kpc;
108+
109+
double PerseusArm2_B_gauss;
110+
double PerseusArm2_pitch_deg;
111+
double PerseusArm2_phi0_deg;
112+
double PerseusArm2_x_shift_kpc;
113+
double PerseusArm2_y_shift_kpc;
114+
double PerseusArm2_arc_radius1_kpc;
115+
double PerseusArm2_arc_radius2_kpc;
116+
double PerseusArm2_arc_eps;
117+
double PerseusArm2_arc_div_deg;
118+
double PerseusArm2_rmin_kpc;
119+
double PerseusArm2_rmax_kpc;
120+
double PerseusArm2_zmin_kpc;
121+
double PerseusArm2_zmax_kpc;
122+
123+
public:
124+
Vector3d getField(const Vector3d& pos) const;
125+
KST24Field();
126+
127+
Vector3d get_toroidal(const Vector3d pos_kpc, const double tor_B_gauss,
128+
const double tor_zmin_kpc, const double tor_zmax_kpc,
129+
const double tor_rmin_kpc, const double tor_rmax_kpc) const;
130+
131+
Vector3d get_Xfield(const Vector3d pos_kpc, const double Xfield_B_gauss,
132+
const double Xfield_rmin_kpc, const double Xfield_rmax_kpc,
133+
const double Xfield_theta_rad) const;
134+
135+
bool is_LB(const Vector3d pos_kpc, const double LB_rmin_kpc, const double LB_dr_kpc,
136+
const double LB_x0_kpc, const double LB_y0_kpc, const double LB_z0_kpc) const;
137+
138+
Vector3d get_LB(const Vector3d pos_kpc, const double LB_B_gauss,
139+
const double LB_lB_deg, const double LB_bB_deg,
140+
const double LB_rmin_kpc, const double LB_dr_kpc,
141+
const double LB_x0_kpc, const double LB_y0_kpc, const double LB_z0_kpc) const;
142+
143+
Vector3d get_logspiral(const Vector3d pos_kpc, const double B_gauss,
144+
const double pitch_deg, const double phi0_deg,
145+
const double x_shift_kpc, const double y_shift_kpc,
146+
const double arc_radius1_kpc, const double arc_radius2_kpc,
147+
const double arc_eps , const double arc_div_deg,
148+
const double rmin_kpc, const double rmax_kpc,
149+
const double zmin_kpc, const double zmax_kpc) const;
150+
};
151+
152+
}// namespace crpropa
153+
154+
#endif /* _KST24_GMF_H_ */

‎include/crpropa/module/Observer.h

Lines changed: 105 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -240,8 +240,24 @@ class ObserverParticleIdVeto: public ObserverFeature {
240240
This observer is very useful if the time evolution of the particle density is needed. It detects all candidates in lin-spaced, log-spaced, or user-defined time intervals and limits the nextStep of candidates to prevent overshooting of detection intervals.
241241
*/
242242
class ObserverTimeEvolution: public ObserverFeature {
243-
private:
243+
protected:
244+
int nIntervals; // number of time invervals
245+
bool isLogarithmicScaling = false; // enables or disables logarithmic scaling for the intervals
246+
bool doDetListConstruction = true; // enables the construction of detList in the relevant functions (addTime, addTimeRange)
247+
double minimum; // the minimum time
248+
double maximum; // the maximum time
249+
/** Vector containing all used times.
250+
It is only constructed by the user manually.
251+
If it is not empty, the vector will be used instead of the getTime function.
252+
(leave empty if you want to rather use functions)
253+
*/
244254
std::vector<double> detList;
255+
/**
256+
A temporary storage for detList, this enables the return of a List in getTimes
257+
without risking to modify detList
258+
*/
259+
mutable std::vector<double> tempDetList;
260+
245261
public:
246262
/** Default constructor
247263
*/
@@ -250,22 +266,105 @@ class ObserverTimeEvolution: public ObserverFeature {
250266
@param min minimum time
251267
@param dist time interval for detection
252268
@param numb number of time intervals
269+
270+
This constructor calculates the maximum from max = min + (numb - 1) * dist
253271
*/
254272
ObserverTimeEvolution(double min, double dist, double numb);
255273
/** Constructor
256274
@param min minimum time
257275
@param max maximum time
258276
@param numb number of time intervals
259277
@param log log (input: true) or lin (input: false) scaling between min and max with numb steps
278+
279+
This constructor sets the maximum directly and gets numb automatically.
280+
You need to set the log parameter, since an overload for the first three doubles exist.
260281
*/
261282
ObserverTimeEvolution(double min, double max, double numb, bool log);
262-
// Add a new time step to the detection time list of the observer
263-
void addTime(const double &position);
264-
// Using log or lin spacing of times in the range between min and
265-
// max for observing particles
283+
/** Constructor
284+
@param detList user defined vector<double> with times to check
285+
286+
This constructor uses a predefined vector containing the times that should be observed.
287+
The so created detList can then be modified via addTime, addTimeRange and setTimes.
288+
*/
289+
ObserverTimeEvolution(const std::vector<double> &detList);
290+
/** Destructor
291+
*/
292+
~ObserverTimeEvolution(){}
293+
294+
/** Function
295+
Generates the detList if it is empty when for example the
296+
ObserverTimeEvolution(const std::vector<double> &detList) constructor
297+
was not used.
298+
Use this function to create a detList with can then be modified.
299+
When detList is not empty its entries are used instead of a runtime calculation of the times.
300+
*/
301+
void constructDetListIfEmpty();
302+
/** Function
303+
@param candidate Candidate usally given by a module list
304+
305+
Checks whether to make a detection at the current step of candidate or not.
306+
This function is called in Observer.process with the simulated Candidate.
307+
*/
308+
DetectionState checkDetection(Candidate *candidate) const;
309+
/** Function
310+
@param enableConstruction if true, constructs detList from range of min, max, numb
311+
when calling addTime
312+
Clears the content of detList
313+
*/
314+
void clear();
315+
/** Function
316+
Checks if detList is empty
317+
*/
318+
bool empty(){return detList.empty();}
319+
320+
// setter functions:
321+
/** Function
322+
@param time Time that should be appended to detList
323+
324+
Makes a push_back on detList with the given time.
325+
*/
326+
void addTime(const double &time);
327+
/** Function
328+
@param min minimum time
329+
@param max maximum time
330+
@param numb number of time intervals
331+
@param log log (input: true) or lin (input: false) scaling between min and max with numb steps
332+
333+
Appends a linear or logarithmic time range to detList via repeatedly calling addTime
334+
*/
266335
void addTimeRange(double min, double max, double numb, bool log = false);
336+
/** Function
337+
@param detList vector<double> containing times when to observe
338+
339+
Sets this->detList to detList, with this it is possible to fully modify detList
340+
*/
341+
void setTimes(const std::vector<double> &detList);
342+
void setMinimum(double min);
343+
void setMaximum(double max){this->maximum = max;}
344+
void setNIntervals(int numb){this->nIntervals = numb;}
345+
void setIsLogarithmicScaling(bool log){this->isLogarithmicScaling = log;}
346+
347+
// getter functions:
348+
/** Function
349+
@param index index of the required time
350+
351+
Replaces the previous preconstructed detList and return the time for a specific index
352+
*/
353+
virtual double getTime(std::size_t index) const;
354+
double getMinimum() const {return minimum;}
355+
double getMaximum() const {return maximum;}
356+
int getNIntervals() const {return nIntervals;}
357+
bool getIsLogarithmicScaling() const {return isLogarithmicScaling;}
358+
/** Function
359+
Returns a vector<double> containing all times between min and max generated by getTime.
360+
This function does not return detList directly, it rather appends a new vector with
361+
getTime numb times.
362+
*/
267363
const std::vector<double>& getTimes() const;
268-
DetectionState checkDetection(Candidate *candidate) const;
364+
/** Function
365+
Returns a string containing a representation of all times.
366+
This function does not create a detList.
367+
*/
269368
std::string getDescription() const;
270369
};
271370

‎python/2_headers.i

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -346,6 +346,7 @@
346346
%include "crpropa/magneticField/PT11Field.h"
347347
%include "crpropa/magneticField/TF17Field.h"
348348
%include "crpropa/magneticField/UF23Field.h"
349+
%include "crpropa/magneticField/KST24Field.h"
349350
%include "crpropa/magneticField/ArchimedeanSpiralField.h"
350351
%include "crpropa/magneticField/CMZField.h"
351352
%include "crpropa/magneticField/turbulentField/TurbulentField.h"

‎src/ModuleList.cpp.in

Lines changed: 4 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,17 +1,13 @@
11
#include "crpropa/ModuleList.h"
22
#include "crpropa/ProgressBar.h"
33

4-
#if _OPENMP
5-
#include <omp.h>
6-
#define OMP_SCHEDULE @OMP_SCHEDULE@
7-
#endif
8-
9-
#include <algorithm>
10-
#include <csignal>
11-
#include <bits/stdc++.h>
124
#ifndef sighandler_t
135
typedef void (*sighandler_t)(int);
146
#endif
7+
#ifdef _OPENMP
8+
#include <omp.h>
9+
#define OMP_SCHEDULE @OMP_SCHEDULE@
10+
#endif
1511

1612
using namespace std;
1713

0 commit comments

Comments
 (0)