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 4 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
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];
Copy link
Member

Choose a reason for hiding this comment

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

alignment

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