Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Uf23 #490

Merged
merged 6 commits into from
Jun 3, 2024
Merged

Uf23 #490

Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
27 changes: 14 additions & 13 deletions CHANGELOG.md
Original file line number Diff line number Diff line change
@@ -1,19 +1,20 @@
## CRPropa vNext

### Bug fixes:
* Fixed sign for exponential decay of magn. field strength with Galactic height in LogarithmicSpiralField
* Fixed r term in source distribution for SNR and Pulsar
* Fixed sign for exponential decay of magn. field strength with Galactic height in LogarithmicSpiralField
* Fixed r term in source distribution for SNR and Pulsar
* Fixed wrong mass inheritance for secondaries other than nuclei or electron/positron

### New features:
* Added new backwards-compatible function particleMass that returns particle mass also for non-nuclei
* EBL model from Finke et al. 2022
* Added the new Galactic magnetic field models from Unger&Farrar arXiv:2311.12120
* Added EBL model from Finke et al. 2022

### Interface changes:

### Features that are deprecated and will be removed after this release

### Removed features
### Removed features
* AMRMagneticField - underlying library (saga) is no longer supported.
* ObserverPoint: Use Observer1D instead.

Expand Down Expand Up @@ -42,8 +43,8 @@
* ObserverPoint will be renamed into Observer1D.
* AMRMagneticField - underlying library (saga) is no longer supported.

### Removed features
* External extensions DINT and Eleca, which can be replaced with the
### Removed features
* External extensions DINT and Eleca, which can be replaced with the
EM*-modules combined with the thinning option for reasonable computation
times.

Expand All @@ -52,20 +53,20 @@
* grplinst
* monopole
* ROOTOutputPlugin


## CRPropa 3.2

### Bug fixes:
* Fix of reflective boundary condition for scalar- and vectorgrids
that showed asymmetry and discontinuities (See issue [#361]).
* Fix in EMTripletPairProduction
* Fix of the data files of the Hackstein EGMF models as well as the
* Fix of the data files of the Hackstein EGMF models as well as the
corresponding example notebook.
* Fix of axis normalization of getRotated in Vector3.h.
* Fix of secondary spectra in electromagnetic interactions
(EM*-modules), issue [#334] and pull request [#15] in crpropa-data.
* Fix weight inheritance for secondary particles; they are created with
* Fix weight inheritance for secondary particles; they are created with
their parents weights as intial weights now.

### New features:
Expand All @@ -75,15 +76,15 @@
and vectorgrids.
* Add the new PolarizedSingleModeMagneticField class for polarized/
helical single mode magnetic field models.
* Add a source feature for targeted emission, following the
* Add a source feature for targeted emission, following the
von-Mises-Fisher distribution
* Updates in SNR and pulsar source distributions
* Updates in SNR and pulsar source distributions

### Interface changes:
* Plane wave and grid turbulence models use same parameter convention now
* Plane wave and grid turbulence models use same parameter convention now

### Features that are deprecated and will be removed after this release
* External extensions DINT and Eleca, which can be replaced with the
* External extensions DINT and Eleca, which can be replaced with the
EM*-modules combined with the thinning option for reasonable computation
times.

Expand Down
5 changes: 3 additions & 2 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -313,7 +313,7 @@ if(DOWNLOAD_DATA)
message("-- Downloading data files from sciebo ~ 73 MB")
file(DOWNLOAD
https://ruhr-uni-bochum.sciebo.de/public.php/webdav/data-${CRPROPA_DATAFILE_VER}.tar.gz-CHECKSUM
${CMAKE_BINARY_DIR}/data-${CRPROPA_DATAFILE_VER}.tar.gz-CHECKSUM
${CMAKE_BINARY_DIR}/data-${CRPROPA_DATAFILE_VER}.tar.gz-CHECKSUM
USERPWD "3juW9sntQX2IWBS")
file(STRINGS ${CMAKE_BINARY_DIR}/data-${CRPROPA_DATAFILE_VER}.tar.gz-CHECKSUM DATA_CHECKSUM LIMIT_COUNT 1 LENGTH_MINIMUM 32 LENGTH_MAXIMUM 32)
file(DOWNLOAD
Expand Down Expand Up @@ -401,6 +401,7 @@ add_library(crpropa SHARED
src/magneticField/turbulentField/PlaneWaveTurbulence.cpp
src/magneticField/turbulentField/SimpleGridTurbulence.cpp
src/magneticField/TF17Field.cpp
src/magneticField/UF23Field.cpp
src/magneticField/CMZField.cpp
src/advectionField/AdvectionField.cpp
src/massDistribution/ConstantDensity.cpp
Expand Down Expand Up @@ -465,7 +466,7 @@ if(ENABLE_PYTHON AND Python_FOUND)
endif(Python_Development_FOUND)


# use Python_INSTALL_PACKAGE_DIR if provided; otherwise, install in Python_SITELIB
# use Python_INSTALL_PACKAGE_DIR if provided; otherwise, install in Python_SITELIB
set(Python_INSTALL_PACKAGE_DIR "${Python_SITELIB}" CACHE PATH "folder in which the python package is installed")
message(STATUS " package install directory: ${Python_INSTALL_PACKAGE_DIR}")

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -38,6 +38,7 @@
"\n",
"# magnetic field setup\n",
"B = JF12Field()\n",
"#B = UF23Field(UF23Field.base) # options: base,neCL,expX,spur,cre10,synCG,twistX,nebCor \n",
"#seed = 691342\n",
"#B.randomStriated(seed)\n",
"#B.randomTurbulent(seed)\n",
Expand Down
1 change: 1 addition & 0 deletions include/CRPropa.h
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@
#include "crpropa/magneticField/PT11Field.h"
#include "crpropa/magneticField/QuimbyMagneticField.h"
#include "crpropa/magneticField/TF17Field.h"
#include "crpropa/magneticField/UF23Field.h"
#include "crpropa/magneticField/CMZField.h"
#include "crpropa/magneticField/turbulentField/GridTurbulence.h"
#include "crpropa/magneticField/turbulentField/HelicalGridTurbulence.h"
Expand Down
153 changes: 153 additions & 0 deletions include/crpropa/magneticField/UF23Field.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,153 @@
#ifndef _UF23Field_h_
#define _UF23Field_h_

#include <vector>
#include "crpropa/magneticField/MagneticField.h"

namespace crpropa {
/**
* \addtogroup MagneticFields
* @{
*/


/**
@class UF23Field
@brief UF23Field Galactic magnetic field model

Implements the eight coherent magnetic field models of UF23
See: M. Unger and G.R. Farrar, arXiv:2311.12120

Assumes a galactocentric coordinate system with the Galactic center
at the origin, the x-axis pointing in the opposite direction of the
Sun, and the z-axis pointing towards Galactic North.

*/

class UF23Field : public MagneticField {
public:
/// model variations (see Tab.2 of UF23 paper)
enum ModelType {
base,
neCL,
expX,
spur,
cre10,
synCG,
twistX,
nebCor
};


public:
/**
@brief constructor
@param mt model type (see Tab.2 of UF23 paper)
@param maxRadiusInKpc maximum radius of field in kpc
*/
UF23Field(const ModelType mt);
/// no default constructor
UF23Field() = delete;
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You probably chose to not provide a default model because none of the models yielded a significantly better fit than the others, right?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think that not providing a default model makes the user aware that this is a collection of 8 models to choose from.


Vector3d getField(const Vector3d& pos) const;

private:

/**
@brief calculate coherent magnetic field at a given position
@param posInKpc position with components given in kpc
@return coherent field in microgauss
*/
Vector3d operator()(const Vector3d& posInKpc) const;

/// model parameters, see Table 3 of UF23 paper
enum EPar {
eDiskB1 = 0,
eDiskB2,
eDiskB3,
eDiskH,
eDiskPhase1,
eDiskPhase2,
eDiskPhase3,
eDiskPitch,
eDiskW,
ePoloidalA,
ePoloidalB,
ePoloidalP,
ePoloidalR,
ePoloidalW,
ePoloidalZ,
ePoloidalXi,
eSpurCenter,
eSpurLength,
eSpurWidth,
eStriation,
eToroidalBN,
eToroidalBS,
eToroidalR,
eToroidalW,
eToroidalZ,
eTwistingTime,
eNpar
};

/// model type given in constructor
const ModelType fModelType;
/// maximum galacto-centric radius beyond which B=0
const double fMaxRadiusSquared;

// parameters are stored in array
double fParameters[eNpar] = { 0 };
// references to parameters for convience
double& fDiskB1 = fParameters[eDiskB1];
double& fDiskB2 = fParameters[eDiskB2];
double& fDiskB3 = fParameters[eDiskB3];
double& fDiskH = fParameters[eDiskH];
double& fDiskPhase1 = fParameters[eDiskPhase1];
double& fDiskPhase2 = fParameters[eDiskPhase2];
double& fDiskPhase3 = fParameters[eDiskPhase3];
double& fDiskPitch = fParameters[eDiskPitch];
double& fDiskW = fParameters[eDiskW];
double& fPoloidalA = fParameters[ePoloidalA];
double& fPoloidalB = fParameters[ePoloidalB];
double& fPoloidalP = fParameters[ePoloidalP];
double& fPoloidalR = fParameters[ePoloidalR];
double& fPoloidalW = fParameters[ePoloidalW];
double& fPoloidalZ = fParameters[ePoloidalZ];
double& fPoloidalXi = fParameters[ePoloidalXi];
double& fSpurCenter = fParameters[eSpurCenter];
double& fSpurLength = fParameters[eSpurLength];
double& fSpurWidth = fParameters[eSpurWidth];
double& fStriation = fParameters[eStriation];
double& fToroidalBN = fParameters[eToroidalBN];
double& fToroidalBS = fParameters[eToroidalBS];
double& fToroidalR = fParameters[eToroidalR];
double& fToroidalW = fParameters[eToroidalW];
double& fToroidalZ = fParameters[eToroidalZ];
double& fTwistingTime = fParameters[eTwistingTime];

// some pre-calculated derived parameter values
double fSinPitch = 0;
double fCosPitch = 0;
double fTanPitch = 0;

/// major field components
Vector3d getDiskField(const Vector3d& pos) const;
Vector3d getHaloField(const Vector3d& pos) const;

/// sub-components depending on model type
/// -- Sec. 5.2.2
Vector3d getSpiralField(const double x, const double y, const double z) const;
/// -- Sec. 5.2.3
Vector3d getSpurField(const double x, const double y, const double z) const;
/// -- Sec. 5.3.1
Vector3d getToroidalHaloField(const double x, const double y, const double z) const;
/// -- Sec. 5.3.2
Vector3d getPoloidalHaloField(const double x, const double y, const double z) const;
/// -- Sec. 5.3.3
Vector3d getTwistedHaloField(const double x, const double y, const double z) const;

};
/** @} */
} // namespace crpropa
#endif
12 changes: 5 additions & 7 deletions python/2_headers.i
Original file line number Diff line number Diff line change
Expand Up @@ -340,6 +340,7 @@
%include "crpropa/magneticField/PolarizedSingleModeMagneticField.h"
%include "crpropa/magneticField/PT11Field.h"
%include "crpropa/magneticField/TF17Field.h"
%include "crpropa/magneticField/UF23Field.h"
%include "crpropa/magneticField/ArchimedeanSpiralField.h"
%include "crpropa/magneticField/CMZField.h"
%include "crpropa/magneticField/turbulentField/TurbulentField.h"
Expand Down Expand Up @@ -433,7 +434,7 @@
crpropa::ModuleList::iterator _end) :
cur(_cur), end(_end) {
}
ModuleListIterator* __iter__() {
ModuleListIterator* __iter__() {
return this;
}
crpropa::ModuleList::iterator cur;
Expand All @@ -454,7 +455,7 @@
ModuleListIterator __iter__() {
return ModuleListIterator($self->begin(), $self->end());
}

crpropa::ref_ptr<crpropa::Module> __getitem__(size_t i) {
if (i >= $self->size()) {
throw RangeError();
Expand All @@ -480,8 +481,8 @@
crpropa::ParticleCollector::iterator _end) :
cur(_cur), end(_end) {
}
ParticleCollectorIterator* __iter__() {
return this;
ParticleCollectorIterator* __iter__() {
return this;
}
crpropa::ParticleCollector::iterator cur;
crpropa::ParticleCollector::iterator end;
Expand Down Expand Up @@ -548,6 +549,3 @@
%template(StepLengthModifierRefPtr) crpropa::ref_ptr<crpropa::StepLengthModifier>;
%feature("director") crpropa::StepLengthModifier;
%include "crpropa/module/Acceleration.h"



Loading