diff --git a/CHANGELOG.md b/CHANGELOG.md index ed7d14f33..129e75637 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -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. @@ -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. @@ -52,7 +53,7 @@ * grplinst * monopole * ROOTOutputPlugin - + ## CRPropa 3.2 @@ -60,12 +61,12 @@ * 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: @@ -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. diff --git a/CMakeLists.txt b/CMakeLists.txt index 5431a940a..43483efaf 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -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 @@ -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 @@ -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}") diff --git a/doc/pages/example_notebooks/galactic_backtracking/galactic_backtracking.ipynb b/doc/pages/example_notebooks/galactic_backtracking/galactic_backtracking.ipynb index 2058d32c1..cef88981a 100644 --- a/doc/pages/example_notebooks/galactic_backtracking/galactic_backtracking.ipynb +++ b/doc/pages/example_notebooks/galactic_backtracking/galactic_backtracking.ipynb @@ -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", diff --git a/include/CRPropa.h b/include/CRPropa.h index abd06532a..bac540bb9 100644 --- a/include/CRPropa.h +++ b/include/CRPropa.h @@ -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" diff --git a/include/crpropa/magneticField/UF23Field.h b/include/crpropa/magneticField/UF23Field.h new file mode 100644 index 000000000..462ed3142 --- /dev/null +++ b/include/crpropa/magneticField/UF23Field.h @@ -0,0 +1,153 @@ +#ifndef _UF23Field_h_ +#define _UF23Field_h_ + +#include +#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; + + 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 diff --git a/python/2_headers.i b/python/2_headers.i index 6aaad3ac7..d34665838 100644 --- a/python/2_headers.i +++ b/python/2_headers.i @@ -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" @@ -433,7 +434,7 @@ crpropa::ModuleList::iterator _end) : cur(_cur), end(_end) { } - ModuleListIterator* __iter__() { + ModuleListIterator* __iter__() { return this; } crpropa::ModuleList::iterator cur; @@ -454,7 +455,7 @@ ModuleListIterator __iter__() { return ModuleListIterator($self->begin(), $self->end()); } - + crpropa::ref_ptr __getitem__(size_t i) { if (i >= $self->size()) { throw RangeError(); @@ -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; @@ -548,6 +549,3 @@ %template(StepLengthModifierRefPtr) crpropa::ref_ptr; %feature("director") crpropa::StepLengthModifier; %include "crpropa/module/Acceleration.h" - - - diff --git a/src/magneticField/UF23Field.cpp b/src/magneticField/UF23Field.cpp new file mode 100644 index 000000000..6a65b1b80 --- /dev/null +++ b/src/magneticField/UF23Field.cpp @@ -0,0 +1,558 @@ +#include "crpropa/magneticField/UF23Field.h" + +#include +#include +#include +#include + +// local helper functions and constants +namespace uf23 { + + template + crpropa::Vector3d CylToCart(const T v, const double cosPhi, const double sinPhi) + { + return crpropa::Vector3d(v[0] * cosPhi - v[1] * sinPhi, + v[0] * sinPhi + v[1] * cosPhi, + v[2]); + } + + template + crpropa::Vector3d CartToCyl(const T v, const double cosPhi, const double sinPhi) + { + return crpropa::Vector3d(v[0] * cosPhi + v[1] * sinPhi, + -v[0] * sinPhi + v[1] * cosPhi, + v[2]); + } + + // logistic sigmoid function + inline + double + Sigmoid(const double x, const double x0, const double w) + { + return 1 / (1 + exp(-(x-x0)/w)); + } + + // angle between v0 = (cos(phi0), sin(phi0)) and v1 = (cos(phi1), sin(phi1)) + inline + double + DeltaPhi(const double phi0, const double phi1) + { + return acos(cos(phi1)*cos(phi0) + sin(phi1)*sin(phi0)); + } + + // Interal units used in this code. + // Convert to crpropa with e.g. uf23::kpc / crpropa::kpc. + // The conversion is, however, only needed in the single non-private + // getField() method, all other functions use uf23 units. + const double kPi = 3.1415926535897932384626; + const double kTwoPi = 2*kPi; + const double degree = kPi/180.; + const double kpc = 1; + const double microgauss = 1; + const double megayear = 1; + const double Gpc = 1e6*kpc; + const double pc = 1e-3*kpc; + const double second = megayear / (1e6*60*60*24*365.25); + const double kilometer = kpc / 3.0856775807e+16; +} + +namespace crpropa { +UF23Field::UF23Field(const ModelType mt) : + fModelType(mt), + fMaxRadiusSquared(pow(30*uf23::kpc, 2)) +{ + + // all but expX model have a-->\infty, Eq.(38) + fPoloidalA = 1 * Gpc; + + switch (fModelType) { + // --------------------------------------------- + case base: { + fDiskB1 = 1.0878565e+00 * uf23::microgauss; + fDiskB2 = 2.6605034e+00 * uf23::microgauss; + fDiskB3 = 3.1166311e+00 * uf23::microgauss; + fDiskH = 7.9408965e-01 * uf23::kpc; + fDiskPhase1 = 2.6316589e+02 * uf23::degree; + fDiskPhase2 = 9.7782269e+01 * uf23::degree; + fDiskPhase3 = 3.5112281e+01 * uf23::degree; + fDiskPitch = 1.0106900e+01 * uf23::degree; + fDiskW = 1.0720909e-01 * uf23::kpc; + fPoloidalB = 9.7775487e-01 * uf23::microgauss; + fPoloidalP = 1.4266186e+00 * uf23::kpc; + fPoloidalR = 7.2925417e+00 * uf23::kpc; + fPoloidalW = 1.1188158e-01 * uf23::kpc; + fPoloidalZ = 4.4597373e+00 * uf23::kpc; + fStriation = 3.4557571e-01; + fToroidalBN = 3.2556760e+00 * uf23::microgauss; + fToroidalBS = -3.0914569e+00 * uf23::microgauss; + fToroidalR = 1.0193815e+01 * uf23::kpc; + fToroidalW = 1.6936993e+00 * uf23::kpc; + fToroidalZ = 4.0242749e+00 * uf23::kpc; + break; + } + case cre10: { + // --------------------------------------------- + fDiskB1 = 1.2035697e+00 * uf23::microgauss; + fDiskB2 = 2.7478490e+00 * uf23::microgauss; + fDiskB3 = 3.2104342e+00 * uf23::microgauss; + fDiskH = 8.0844932e-01 * uf23::kpc; + fDiskPhase1 = 2.6515882e+02 * uf23::degree; + fDiskPhase2 = 9.8211313e+01 * uf23::degree; + fDiskPhase3 = 3.5944588e+01 * uf23::degree; + fDiskPitch = 1.0162759e+01 * uf23::degree; + fDiskW = 1.0824003e-01 * uf23::kpc; + fPoloidalB = 9.6938453e-01 * uf23::microgauss; + fPoloidalP = 1.4150957e+00 * uf23::kpc; + fPoloidalR = 7.2987296e+00 * uf23::kpc; + fPoloidalW = 1.0923051e-01 * uf23::kpc; + fPoloidalZ = 4.5748332e+00 * uf23::kpc; + fStriation = 2.4950386e-01; + fToroidalBN = 3.7308133e+00 * uf23::microgauss; + fToroidalBS = -3.5039958e+00 * uf23::microgauss; + fToroidalR = 1.0407507e+01 * uf23::kpc; + fToroidalW = 1.7398375e+00 * uf23::kpc; + fToroidalZ = 2.9272800e+00 * uf23::kpc; + break; + } + case nebCor: { + // --------------------------------------------- + fDiskB1 = 1.4081935e+00 * uf23::microgauss; + fDiskB2 = 3.5292400e+00 * uf23::microgauss; + fDiskB3 = 4.1290147e+00 * uf23::microgauss; + fDiskH = 8.1151971e-01 * uf23::kpc; + fDiskPhase1 = 2.6447529e+02 * uf23::degree; + fDiskPhase2 = 9.7572660e+01 * uf23::degree; + fDiskPhase3 = 3.6403798e+01 * uf23::degree; + fDiskPitch = 1.0151183e+01 * uf23::degree; + fDiskW = 1.1863734e-01 * uf23::kpc; + fPoloidalB = 1.3485916e+00 * uf23::microgauss; + fPoloidalP = 1.3414395e+00 * uf23::kpc; + fPoloidalR = 7.2473841e+00 * uf23::kpc; + fPoloidalW = 1.4318227e-01 * uf23::kpc; + fPoloidalZ = 4.8242603e+00 * uf23::kpc; + fStriation = 3.8610837e-10; + fToroidalBN = 4.6491142e+00 * uf23::microgauss; + fToroidalBS = -4.5006610e+00 * uf23::microgauss; + fToroidalR = 1.0205288e+01 * uf23::kpc; + fToroidalW = 1.7004868e+00 * uf23::kpc; + fToroidalZ = 3.5557767e+00 * uf23::kpc; + break; + } + case neCL: { + // --------------------------------------------- + fDiskB1 = 1.4259645e+00 * uf23::microgauss; + fDiskB2 = 1.3543223e+00 * uf23::microgauss; + fDiskB3 = 3.4390669e+00 * uf23::microgauss; + fDiskH = 6.7405199e-01 * uf23::kpc; + fDiskPhase1 = 1.9961898e+02 * uf23::degree; + fDiskPhase2 = 1.3541461e+02 * uf23::degree; + fDiskPhase3 = 6.4909767e+01 * uf23::degree; + fDiskPitch = 1.1867859e+01 * uf23::degree; + fDiskW = 6.1162799e-02 * uf23::kpc; + fPoloidalB = 9.8387831e-01 * uf23::microgauss; + fPoloidalP = 1.6773615e+00 * uf23::kpc; + fPoloidalR = 7.4084361e+00 * uf23::kpc; + fPoloidalW = 1.4168192e-01 * uf23::kpc; + fPoloidalZ = 3.6521188e+00 * uf23::kpc; + fStriation = 3.3600213e-01; + fToroidalBN = 2.6256593e+00 * uf23::microgauss; + fToroidalBS = -2.5699466e+00 * uf23::microgauss; + fToroidalR = 1.0134257e+01 * uf23::kpc; + fToroidalW = 1.1547728e+00 * uf23::kpc; + fToroidalZ = 4.5585463e+00 * uf23::kpc; + break; + } + case spur: { + // --------------------------------------------- + fDiskB1 = -4.2993328e+00 * uf23::microgauss; + fDiskH = 7.5019749e-01 * uf23::kpc; + fDiskPhase1 = 1.5589875e+02 * uf23::degree; + fDiskPitch = 1.2074432e+01 * uf23::degree; + fDiskW = 1.2263120e-01 * uf23::kpc; + fPoloidalB = 9.9302987e-01 * uf23::microgauss; + fPoloidalP = 1.3982374e+00 * uf23::kpc; + fPoloidalR = 7.1973387e+00 * uf23::kpc; + fPoloidalW = 1.2262244e-01 * uf23::kpc; + fPoloidalZ = 4.4853270e+00 * uf23::kpc; + fSpurCenter = 1.5718686e+02 * uf23::degree; + fSpurLength = 3.1839577e+01 * uf23::degree; + fSpurWidth = 1.0318114e+01 * uf23::degree; + fStriation = 3.3022369e-01; + fToroidalBN = 2.9286724e+00 * uf23::microgauss; + fToroidalBS = -2.5979895e+00 * uf23::microgauss; + fToroidalR = 9.7536425e+00 * uf23::kpc; + fToroidalW = 1.4210055e+00 * uf23::kpc; + fToroidalZ = 6.0941229e+00 * uf23::kpc; + break; + } + case synCG: { + // --------------------------------------------- + fDiskB1 = 8.1386878e-01 * uf23::microgauss; + fDiskB2 = 2.0586930e+00 * uf23::microgauss; + fDiskB3 = 2.9437335e+00 * uf23::microgauss; + fDiskH = 6.2172353e-01 * uf23::kpc; + fDiskPhase1 = 2.2988551e+02 * uf23::degree; + fDiskPhase2 = 9.7388282e+01 * uf23::degree; + fDiskPhase3 = 3.2927367e+01 * uf23::degree; + fDiskPitch = 9.9034844e+00 * uf23::degree; + fDiskW = 6.6517521e-02 * uf23::kpc; + fPoloidalB = 8.0883734e-01 * uf23::microgauss; + fPoloidalP = 1.5820957e+00 * uf23::kpc; + fPoloidalR = 7.4625235e+00 * uf23::kpc; + fPoloidalW = 1.5003765e-01 * uf23::kpc; + fPoloidalZ = 3.5338550e+00 * uf23::kpc; + fStriation = 6.3434763e-01; + fToroidalBN = 2.3991193e+00 * uf23::microgauss; + fToroidalBS = -2.0919944e+00 * uf23::microgauss; + fToroidalR = 9.4227834e+00 * uf23::kpc; + fToroidalW = 9.1608418e-01 * uf23::kpc; + fToroidalZ = 5.5844594e+00 * uf23::kpc; + break; + } + case twistX: { + // --------------------------------------------- + fDiskB1 = 1.3741995e+00 * uf23::microgauss; + fDiskB2 = 2.0089881e+00 * uf23::microgauss; + fDiskB3 = 1.5212463e+00 * uf23::microgauss; + fDiskH = 9.3806180e-01 * uf23::kpc; + fDiskPhase1 = 2.3560316e+02 * uf23::degree; + fDiskPhase2 = 1.0189856e+02 * uf23::degree; + fDiskPhase3 = 5.6187572e+01 * uf23::degree; + fDiskPitch = 1.2100979e+01 * uf23::degree; + fDiskW = 1.4933338e-01 * uf23::kpc; + fPoloidalB = 6.2793114e-01 * uf23::microgauss; + fPoloidalP = 2.3292519e+00 * uf23::kpc; + fPoloidalR = 7.9212358e+00 * uf23::kpc; + fPoloidalW = 2.9056201e-01 * uf23::kpc; + fPoloidalZ = 2.6274437e+00 * uf23::kpc; + fStriation = 7.7616317e-01; + fTwistingTime = 5.4733549e+01 * uf23::megayear; + break; + } + case expX: { + // --------------------------------------------- + fDiskB1 = 9.9258148e-01 * uf23::microgauss; + fDiskB2 = 2.1821124e+00 * uf23::microgauss; + fDiskB3 = 3.1197345e+00 * uf23::microgauss; + fDiskH = 7.1508681e-01 * uf23::kpc; + fDiskPhase1 = 2.4745741e+02 * uf23::degree; + fDiskPhase2 = 9.8578879e+01 * uf23::degree; + fDiskPhase3 = 3.4884485e+01 * uf23::degree; + fDiskPitch = 1.0027070e+01 * uf23::degree; + fDiskW = 9.8524736e-02 * uf23::kpc; + fPoloidalA = 6.1938701e+00 * uf23::kpc; + fPoloidalB = 5.8357990e+00 * uf23::microgauss; + fPoloidalP = 1.9510779e+00 * uf23::kpc; + fPoloidalR = 2.4994376e+00 * uf23::kpc; + // internally, xi is fitted and z = tan(xi)*a + fPoloidalXi = 2.0926122e+01 * uf23::degree; + fPoloidalZ = fPoloidalA*tan(fPoloidalXi); + fStriation = 5.1440500e-01; + fToroidalBN = 2.7077434e+00 * uf23::microgauss; + fToroidalBS = -2.5677104e+00 * uf23::microgauss; + fToroidalR = 1.0134022e+01 * uf23::kpc; + fToroidalW = 2.0956159e+00 * uf23::kpc; + fToroidalZ = 5.4564991e+00 * uf23::kpc; + break; + } + default: { + throw std::runtime_error("unknown field model"); + break; + } + } + + fSinPitch = sin(fDiskPitch); + fCosPitch = cos(fDiskPitch); + fTanPitch = tan(fDiskPitch); + +} + +Vector3d +UF23Field::getField(const Vector3d &position) + const +{ + Vector3d posInKpc = position / kpc; + return (this->operator()(posInKpc)) / uf23::microgauss * microgauss; +} + + +Vector3d +UF23Field::operator()(const Vector3d& posInKpc) + const +{ + const auto pos = posInKpc * uf23::kpc; + if (pos.getR2() > fMaxRadiusSquared) + return Vector3d(0, 0, 0); + else { + const auto diskField = getDiskField(pos); + const auto haloField = getHaloField(pos); + return (diskField + haloField) / uf23::microgauss; + } +} + +Vector3d +UF23Field::getDiskField(const Vector3d& pos) + const +{ + if (fModelType == spur) + return getSpurField(pos.x, pos.y, pos.z); + else + return getSpiralField(pos.x, pos.y, pos.z); +} + + +Vector3d +UF23Field::getHaloField(const Vector3d& pos) + const +{ + if (fModelType == twistX) + return getTwistedHaloField(pos.x, pos.y, pos.z); + else + return + getToroidalHaloField(pos.x, pos.y, pos.z) + + getPoloidalHaloField(pos.x, pos.y, pos.z); +} + + +Vector3d +UF23Field::getTwistedHaloField(const double x, const double y, const double z) + const +{ + const double r = sqrt(x*x + y*y); + const double cosPhi = r > std::numeric_limits::min() ? x / r : 1; + const double sinPhi = r > std::numeric_limits::min() ? y / r : 0; + + const Vector3d bXCart = getPoloidalHaloField(x, y, z); + const double bXCartTmp[3] = {bXCart.x, bXCart.y, bXCart.z}; + const Vector3d bXCyl = uf23::CartToCyl(bXCartTmp, cosPhi, sinPhi); + + const double bZ = bXCyl.z; + const double bR = bXCyl.x; + + double bPhi = 0; + + if (fTwistingTime != 0 && r != 0) { + // radial rotation curve parameters (fit to Reid et al 2014) + const double v0 = -240 * uf23::kilometer/uf23::second; + const double r0 = 1.6 * uf23::kpc; + // vertical gradient (Levine+08) + const double z0 = 10 * uf23::kpc; + + // Eq.(43) + const double fr = 1 - exp(-r/r0); + // Eq.(44) + const double t0 = exp(2*std::abs(z)/z0); + const double gz = 2 / (1 + t0); + + // Eq. (46) + const double signZ = z < 0 ? -1 : 1; + const double deltaZ = -signZ * v0 * fr / z0 * t0 * pow(gz, 2); + // Eq. (47) + const double deltaR = v0 * ((1-fr)/r0 - fr/r) * gz; + + // Eq.(45) + bPhi = (bZ * deltaZ + bR * deltaR) * fTwistingTime; + + } + const double bCylX[3] = {bR, bPhi , bZ}; + return uf23::CylToCart(bCylX, cosPhi, sinPhi); +} + +Vector3d +UF23Field::getToroidalHaloField(const double x, const double y, const double z) + const +{ + const double r2 = x*x + y*y; + const double r = sqrt(r2); + const double absZ = std::abs(z); + + const double b0 = z >= 0 ? fToroidalBN : fToroidalBS; + const double rh = fToroidalR; + const double z0 = fToroidalZ; + const double fwh = fToroidalW; + const double sigmoidR = uf23::Sigmoid(r, rh, fwh); + const double sigmoidZ = uf23::Sigmoid(absZ, fDiskH, fDiskW); + + // Eq. (21) + const double bPhi = b0 * (1. - sigmoidR) * sigmoidZ * exp(-absZ/z0); + + const double bCyl[3] = {0, bPhi, 0}; + const double cosPhi = r > std::numeric_limits::min() ? x / r : 1; + const double sinPhi = r > std::numeric_limits::min() ? y / r : 0; + return uf23::CylToCart(bCyl, cosPhi, sinPhi); +} + +Vector3d +UF23Field::getPoloidalHaloField(const double x, const double y, const double z) + const +{ + const double r2 = x*x + y*y; + const double r = sqrt(r2); + + const double c = pow(fPoloidalA/fPoloidalZ, fPoloidalP); + const double a0p = pow(fPoloidalA, fPoloidalP); + const double rp = pow(r, fPoloidalP); + const double abszp = pow(std::abs(z), fPoloidalP); + const double cabszp = c*abszp; + + /* + since $\sqrt{a^2 + b} - a$ is numerical unstable for $b\ll a$, + we use $(\sqrt{a^2 + b} - a) \frac{\sqrt{a^2 + b} + a}{\sqrt{a^2 + + b} + a} = \frac{b}{\sqrt{a^2 + b} + a}$} + */ + + const double t0 = a0p + cabszp - rp; + const double t1 = sqrt(pow(t0, 2) + 4*a0p*rp); + const double ap = 2*a0p*rp / (t1 + t0); + + double a = 0; + if (ap < 0) { + if (r > std::numeric_limits::min()) { + // this should never happen + throw std::runtime_error("ap = " + std::to_string(ap)); + } + else + a = 0; + } + else + a = pow(ap, 1/fPoloidalP); + + // Eq.(29) and Eq.(32) + const double radialDependence = + fModelType == expX ? + exp(-a/fPoloidalR) : + 1 - uf23::Sigmoid(a, fPoloidalR, fPoloidalW); + + // Eq.(28) + const double Bzz = fPoloidalB * radialDependence; + + // (r/a) + const double rOverA = 1 / pow(2*a0p / (t1 + t0), 1/fPoloidalP); + + // Eq.(35) for p=n + const double signZ = z < 0 ? -1 : 1; + const double Br = + Bzz * c * a / rOverA * signZ * pow(std::abs(z), fPoloidalP - 1) / t1; + + // Eq.(36) for p=n + const double Bz = Bzz * pow(rOverA, fPoloidalP-2) * (ap + a0p) / t1; + + if (r < std::numeric_limits::min()) + return Vector3d(0, 0, Bz); + else { + const double bCylX[3] = {Br, 0 , Bz}; + const double cosPhi = x / r; + const double sinPhi = y / r; + return uf23::CylToCart(bCylX, cosPhi, sinPhi); + } +} + +Vector3d +UF23Field::getSpurField(const double x, const double y, const double z) + const +{ + // reference approximately at solar radius + const double rRef = 8.2*uf23::kpc; + + // cylindrical coordinates + const double r2 = x*x + y*y; + const double r = sqrt(r2); + if (r < std::numeric_limits::min()) + return Vector3d(0, 0, 0); + + double phi = atan2(y, x); + if (phi < 0) + phi += uf23::kTwoPi; + + const double phiRef = fDiskPhase1; + int iBest = -2; + double bestDist = -1; + for (int i = -1; i <= 1; ++i) { + const double pphi = phi - phiRef + i*uf23::kTwoPi; + const double rr = rRef*exp(pphi * fTanPitch); + if (bestDist < 0 || std::abs(r-rr) < bestDist) { + bestDist = std::abs(r-rr); + iBest = i; + } + } + if (iBest == 0) { + const double phi0 = phi - log(r/rRef) / fTanPitch; + + // Eq. (16) + const double deltaPhi0 = uf23::DeltaPhi(phiRef, phi0); + const double delta = deltaPhi0 / fSpurWidth; + const double B = fDiskB1 * exp(-0.5*pow(delta, 2)); + + // Eq. (18) + const double wS = 5*uf23::degree; + const double phiC = fSpurCenter; + const double deltaPhiC = uf23::DeltaPhi(phiC, phi); + const double lC = fSpurLength; + const double gS = 1 - uf23::Sigmoid(std::abs(deltaPhiC), lC, wS); + + // Eq. (13) + const double hd = 1 - uf23::Sigmoid(std::abs(z), fDiskH, fDiskW); + + // Eq. (17) + const double bS = rRef/r * B * hd * gS; + const double bCyl[3] = {bS * fSinPitch, bS * fCosPitch, 0}; + const double cosPhi = x / r; + const double sinPhi = y / r; + return uf23::CylToCart(bCyl, cosPhi, sinPhi); + } + else + return Vector3d(0, 0, 0); + +} + +Vector3d +UF23Field::getSpiralField(const double x, const double y, const double z) + const +{ + // reference radius + const double rRef = 5*uf23::kpc; + // inner boundary of spiral field + const double rInner = 5*uf23::kpc; + const double wInner = 0.5*uf23::kpc; + // outer boundary of spiral field + const double rOuter = 20*uf23::kpc; + const double wOuter = 0.5*uf23::kpc; + + // cylindrical coordinates + const double r2 = x*x + y*y; + if (r2 == 0) + return Vector3d(0, 0, 0); + const double r = sqrt(r2); + const double phi = atan2(y, x); + + // Eq.(13) + const double hdz = 1 - uf23::Sigmoid(std::abs(z), fDiskH, fDiskW); + + // Eq.(14) times rRef divided by r + const double rFacI = uf23::Sigmoid(r, rInner, wInner); + const double rFacO = 1 - uf23::Sigmoid(r, rOuter, wOuter); + // (using lim r--> 0 (1-exp(-r^2))/r --> r - r^3/2 + ...) + const double rFac = r > 1e-5*uf23::pc ? (1-exp(-r*r)) / r : r * (1 - r2/2); + const double gdrTimesRrefByR = rRef * rFac * rFacO * rFacI; + + // Eq. (12) + const double phi0 = phi - log(r/rRef) / fTanPitch; + + // Eq. (10) + const double b = + fDiskB1 * cos(1 * (phi0 - fDiskPhase1)) + + fDiskB2 * cos(2 * (phi0 - fDiskPhase2)) + + fDiskB3 * cos(3 * (phi0 - fDiskPhase3)); + + // Eq. (11) + const double fac = hdz * gdrTimesRrefByR; + const double bCyl[3] = + { b * fac * fSinPitch, + b * fac * fCosPitch, + 0}; + + const double cosPhi = x / r; + const double sinPhi = y / r; + return uf23::CylToCart(bCyl, cosPhi, sinPhi); +} +} diff --git a/test/testMagneticField.cpp b/test/testMagneticField.cpp index 8528c2593..da735edc6 100644 --- a/test/testMagneticField.cpp +++ b/test/testMagneticField.cpp @@ -4,6 +4,7 @@ #include "crpropa/magneticField/CMZField.h" #include "crpropa/magneticField/PolarizedSingleModeMagneticField.h" #include "crpropa/magneticField/GalacticMagneticField.h" +#include "crpropa/magneticField/UF23Field.h" #include "crpropa/Grid.h" #include "crpropa/Units.h" #include "crpropa/Common.h" @@ -105,7 +106,7 @@ TEST(testPeriodicMagneticField, Exceptions) { TEST(testCMZMagneticField, SimpleTest) { ref_ptr field = new CMZField(); - + // check use-Values EXPECT_FALSE(field->getUseMCField()); EXPECT_TRUE(field->getUseICField()); @@ -125,7 +126,7 @@ TEST(testCMZMagneticField, SimpleTest) { TEST(testCMZMagneticField, TestICComponent) { ref_ptr field = new CMZField(); - Vector3d pos(10*pc,15*pc,-5*pc); + Vector3d pos(10*pc,15*pc,-5*pc); // check IC field at given position Vector3d bVec = field->getField(pos); @@ -136,8 +137,8 @@ TEST(testCMZMagneticField, TestICComponent) { } TEST(testCMZMagneticField, TestNTFField){ ref_ptr field = new CMZField(); - Vector3d pos(10*pc,15*pc,-5*pc); - + Vector3d pos(10*pc,15*pc,-5*pc); + // check NFTField at given position Vector3d bVec = field->getNTFField(pos); EXPECT_NEAR(bVec.getR(),1.692*muG, 1e-3*muG); @@ -147,7 +148,7 @@ TEST(testCMZMagneticField, TestNTFField){ } TEST(testCMZMagneticField, TestRaidoArcField){ ref_ptr field = new CMZField(); - Vector3d pos(10*pc,15*pc,-5*pc); + Vector3d pos(10*pc,15*pc,-5*pc); // check RadioArcField at given position Vector3d bVec = field->getRadioArcField(pos); @@ -160,7 +161,7 @@ TEST(testCMZMagneticField, TestRaidoArcField){ TEST(testCMZMagneticField, TestAzimutalComponent){ ref_ptr field = new CMZField(); Vector3d mid(12*pc, 9*pc, 20*pc); - Vector3d pos(9*pc, 10*pc, 25*pc); + Vector3d pos(9*pc, 10*pc, 25*pc); // simple Test for inner part Vector3d bVec = field->BAz(pos, mid, 100, 0.2, 60*pc); @@ -202,6 +203,138 @@ TEST(testLogarithmicSpiralField, SimpleTest) { EXPECT_NEAR(b.z, 0, 1E-10); } + +// UF23 reference values +std::vector> +getUF23ReferenceValues() +{ + using V = Vector3d; + std::vector> retVal; + retVal.push_back({ + V{-1.46860417e+00, 1.64555489e+00, 8.35702311e-01}, + V{0.00000000e+00, -4.22433497e-01, 1.83232985e-01}, + V{-3.00239177e-01, -2.97045767e-01, 1.83232985e-01}, + V{8.58382023e-04, 7.71891049e-03, 9.71705527e-01}, + V{-1.17276875e+00, -2.33013590e-01, 4.10798494e-01}, + V{2.63883569e-01, -8.79631081e-01, 5.54229961e-06}, + V{-1.71971907e-02, -2.00291358e-02, 4.16024875e-02}, + V{6.18960813e-01, -2.14437026e+00, 8.35702315e-01}, + V{-1.50883911e+00, 2.63874909e+00, 1.16578928e-02} + }); + retVal.push_back({ + V{-1.39447625e+00, 1.57135437e+00, 8.65163870e-01}, + V{0.00000000e+00, -4.44318027e-01, 1.54430718e-01}, + V{-3.15695875e-01, -3.12652066e-01, 1.54430718e-01}, + V{2.09678202e-03, 2.63832916e-03, 9.81077691e-01}, + V{-1.08603376e+00, -1.86359136e-01, 3.91730436e-01}, + V{2.11095801e-01, -7.04082667e-01, 8.86591950e-05}, + V{-1.38087843e-02, -1.52011002e-02, 3.06842086e-02}, + V{5.76543516e-01, -2.03544911e+00, 8.65163870e-01}, + V{-3.77778745e-01, 6.89644787e-01, 5.85347897e-02}, + }); + retVal.push_back({ + V{-1.04170166e+00, 1.93212739e+00, 2.95362499e+00}, + V{0.00000000e+00, -5.87995638e-01, 4.66925506e-01}, + V{-4.20802701e-01, -4.10291633e-01, 4.59557272e-01}, + V{7.79982996e-03, 1.50482351e-02, 5.50337531e+00}, + V{-1.34639246e+00, 6.80103301e-02, 8.60002649e-01}, + V{1.37242004e-01, -8.67085225e-01, 1.25106474e-01}, + V{-1.69325979e-02, -2.95454738e-02, 4.72616404e-02}, + V{4.53873303e-01, -2.03292501e+00, 9.95205454e-01}, + V{-1.30599234e+00, 2.25151057e+00, 2.56704152e-01}, + }); + retVal.push_back({ + V{-1.45840221e+00, 1.64227817e+00, 8.41586777e-01}, + V{0.00000000e+00, -6.98341816e-01, 1.84323895e-01}, + V{-4.95342500e-01, -4.92154113e-01, 1.84323895e-01}, + V{-5.27533653e-03, 1.48965596e-02, 9.86107885e-01}, + V{-1.47472798e+00, -3.33905274e-01, 4.11265918e-01}, + V{2.31202234e-01, -7.70722755e-01, 1.12129423e-05}, + V{-1.54200071e-02, -2.22012150e-02, 4.22733003e-02}, + V{5.92201401e-01, -2.10378256e+00, 8.41586782e-01}, + V{4.05057122e-03, -9.76673905e-03, 8.24321504e-03}, + }); + retVal.push_back({ + V{-1.50595160e+00, 1.67869437e+00, 8.29806168e-01}, + V{0.00000000e+00, -2.27289811e-01, 1.86869632e-01}, + V{-1.62291024e-01, -1.59076795e-01, 1.86869632e-01}, + V{6.71536463e-04, 7.88990042e-03, 9.63291786e-01}, + V{-9.35561196e-01, -1.55840707e-01, 4.13609050e-01}, + V{2.68120686e-01, -8.93743434e-01, 3.48927149e-06}, + V{-1.88291055e-02, -1.94245715e-02, 4.29970299e-02}, + V{6.50804557e-01, -2.18817590e+00, 8.29806172e-01}, + V{-1.58450595e+00, 2.77817014e+00, 1.10336683e-02}, + }); + retVal.push_back({ + V{-1.33098823e+00, 1.49556466e+00, 6.88642986e-01}, + V{0.00000000e+00, -4.99342453e-01, 1.16143813e-01}, + V{-3.54225500e-01, -3.51947350e-01, 1.16143813e-01}, + V{2.14352986e-03, 3.57971370e-03, 8.05219599e-01}, + V{-1.15136627e+00, -2.48536546e-01, 2.95713475e-01}, + V{1.18057505e-01, -3.98308301e-01, 9.11299352e-04}, + V{-1.06899006e-02, -1.12335952e-02, 2.33358311e-02}, + V{5.61264936e-01, -1.94601963e+00, 6.88642987e-01}, + V{-1.18200258e+00, 2.01179003e+00, 8.02453884e-02}, + }); + retVal.push_back({ + V{-3.65508814e-01, 4.74456797e-01, 5.76165841e-01}, + V{0.00000000e+00, 0.00000000e+00, 6.36668108e-02}, + V{-3.68784336e-03, -2.20689056e-03, 6.36668108e-02}, + V{-5.03208784e-02, 5.09887480e-02, 6.27665021e-01}, + V{-4.04821314e-01, -1.01426396e-02, 2.06022291e-01}, + V{-2.34287239e-02, -9.72536117e-02, 2.80624948e-02}, + V{-4.42698169e-03, -6.23429869e-03, 1.07555956e-02}, + V{3.25022555e-01, -1.18950549e+00, 5.76163734e-01}, + V{-8.65217092e-01, 1.85270104e+00, 3.73913202e-01}, + }); + retVal.push_back({ + V{-1.92580231e+00, 2.17134243e+00, 1.13723929e+00}, + V{0.00000000e+00, -4.73261099e-01, 2.65974345e-01}, + V{-3.36780079e-01, -3.32368900e-01, 2.65974345e-01}, + V{-5.24558924e-04, 1.51886238e-02, 1.33757258e+00}, + V{-1.47213976e+00, -2.82203258e-01, 5.71920142e-01}, + V{3.54206333e-01, -1.18108962e+00, 1.00077423e-04}, + V{-2.67268842e-02, -2.88588605e-02, 6.38364561e-02}, + V{7.90570955e-01, -2.84039611e+00, 1.13723931e+00}, + V{-1.97316856e+00, 3.44153600e+00, 3.20266742e-02}, + }); + return retVal; +} + + +TEST(testUF23Field, SimpleTest) { + const std::vector models = + { + UF23Field::base, UF23Field::neCL, UF23Field::expX, UF23Field::spur, + UF23Field::cre10, UF23Field::synCG, UF23Field::twistX, UF23Field::nebCor + }; + + using V = Vector3d; + const std::vector testPositions = + { + V{1, 1, 1}, V{0, 0, -8}, V{0.1, -0.1, -8}, V{0.1, 0.1, 0.1}, V{-1, 3, 4}, + V{-10, -3, 2}, V{-10, -10, 20}, V{-4, -2, 1}, V{6, 5, -0.1} + }; + + const std::vector> referenceValues = + getUF23ReferenceValues(); + + const double precision = 1e-6; + + for (unsigned int i = 0; i < models.size(); ++i) { + const auto& model = models[i]; + const UF23Field uf23Field(model); + for (unsigned int j = 0; j < testPositions.size(); ++j) { + const auto& position = testPositions[j]; + const auto val = uf23Field.getField(position*kpc); + const auto& refVal = referenceValues[i][j] * microgauss; + EXPECT_NEAR(val.x, refVal.x, precision); + EXPECT_NEAR(val.y, refVal.y, precision); + EXPECT_NEAR(val.z, refVal.z, precision); + } + } +} + TEST(testPolarizedSingleModeMagneticField, SimpleTest) { PolarizedSingleModeMagneticField B(2, 4, 0.5, Vector3d(1,1,1), Vector3d(0,1,0), Vector3d(1,0,0), "amplitude", "polarization", "elliptical"); Vector3d b = B.getField(Vector3d(1,1,2));