Line data    Source code

       1              : #ifndef _UF23Field_h_
       2              : #define _UF23Field_h_
       3              : 
       4              : #include <vector>
       5              : #include "crpropa/magneticField/MagneticField.h"
       6              : 
       7              : namespace crpropa {
       8              : /**
       9              :  * \addtogroup MagneticFields
      10              :  * @{
      11              :  */
      12              : 
      13              : 
      14              : /**
      15              :  @class UF23Field
      16              :  @brief UF23Field Galactic magnetic field model
      17              : 
      18              :  Implements the eight coherent magnetic field models of UF23
      19              :  See: M. Unger and G.R. Farrar, Astrophys.J. 970 (2024) 95, arXiv:2311.12120
      20              : 
      21              :  Assumes a galactocentric coordinate system with the Galactic center
      22              :  at the origin, the x-axis pointing in the opposite direction of the
      23              :  Sun, and the z-axis pointing towards Galactic North.
      24              : 
      25              :  */
      26              : 
      27            8 : class UF23Field : public MagneticField  {
      28              : public:
      29              :   /// model variations (see Tab.2 of UF23 paper)
      30              :   enum ModelType {
      31              :     base,
      32              :     neCL,
      33              :     expX,
      34              :     spur,
      35              :     cre10,
      36              :     synCG,
      37              :     twistX,
      38              :     nebCor
      39              :   };
      40              : 
      41              : 
      42              : public:
      43              :   /**
      44              :      @brief constructor
      45              :      @param mt model type  (see Tab.2 of UF23 paper)
      46              :      @param maxRadiusInKpc maximum radius of field in kpc
      47              :   */
      48              :   UF23Field(const ModelType mt);
      49              :   /// no default constructor
      50              :   UF23Field() = delete;
      51              : 
      52              :   Vector3d getField(const Vector3d& pos) const;
      53              : 
      54              : private:
      55              : 
      56              :   /**
      57              :      @brief calculate coherent magnetic field at a given position
      58              :      @param posInKpc position with components given in kpc
      59              :      @return coherent field in microgauss
      60              :   */
      61              :   Vector3d operator()(const Vector3d& posInKpc) const;
      62              : 
      63              :   /// model parameters, see Table 3 of UF23 paper
      64              :   enum EPar {
      65              :     eDiskB1 = 0,
      66              :     eDiskB2,
      67              :     eDiskB3,
      68              :     eDiskH,
      69              :     eDiskPhase1,
      70              :     eDiskPhase2,
      71              :     eDiskPhase3,
      72              :     eDiskPitch,
      73              :     eDiskW,
      74              :     ePoloidalA,
      75              :     ePoloidalB,
      76              :     ePoloidalP,
      77              :     ePoloidalR,
      78              :     ePoloidalW,
      79              :     ePoloidalZ,
      80              :     ePoloidalXi,
      81              :     eSpurCenter,
      82              :     eSpurLength,
      83              :     eSpurWidth,
      84              :     eStriation,
      85              :     eToroidalBN,
      86              :     eToroidalBS,
      87              :     eToroidalR,
      88              :     eToroidalW,
      89              :     eToroidalZ,
      90              :     eTwistingTime,
      91              :     eNpar
      92              :   };
      93              : 
      94              :   /// model type given in constructor
      95              :   const ModelType fModelType;
      96              :   /// maximum galacto-centric radius beyond which B=0
      97              :   const double fMaxRadiusSquared;
      98              : 
      99              :   // parameters are stored in array
     100              :   double fParameters[eNpar] = { 0 };
     101              :   // references to parameters for convience
     102              :   double& fDiskB1       = fParameters[eDiskB1];
     103              :   double& fDiskB2       = fParameters[eDiskB2];
     104              :   double& fDiskB3       = fParameters[eDiskB3];
     105              :   double& fDiskH        = fParameters[eDiskH];
     106              :   double& fDiskPhase1   = fParameters[eDiskPhase1];
     107              :   double& fDiskPhase2   = fParameters[eDiskPhase2];
     108              :   double& fDiskPhase3   = fParameters[eDiskPhase3];
     109              :   double& fDiskPitch    = fParameters[eDiskPitch];
     110              :   double& fDiskW        = fParameters[eDiskW];
     111              :   double& fPoloidalA    = fParameters[ePoloidalA];
     112              :   double& fPoloidalB    = fParameters[ePoloidalB];
     113              :   double& fPoloidalP    = fParameters[ePoloidalP];
     114              :   double& fPoloidalR    = fParameters[ePoloidalR];
     115              :   double& fPoloidalW    = fParameters[ePoloidalW];
     116              :   double& fPoloidalZ    = fParameters[ePoloidalZ];
     117              :   double& fPoloidalXi   = fParameters[ePoloidalXi];
     118              :   double& fSpurCenter   = fParameters[eSpurCenter];
     119              :   double& fSpurLength   = fParameters[eSpurLength];
     120              :   double& fSpurWidth    = fParameters[eSpurWidth];
     121              :   double& fStriation    = fParameters[eStriation];
     122              :   double& fToroidalBN   = fParameters[eToroidalBN];
     123              :   double& fToroidalBS   = fParameters[eToroidalBS];
     124              :   double& fToroidalR    = fParameters[eToroidalR];
     125              :   double& fToroidalW    = fParameters[eToroidalW];
     126              :   double& fToroidalZ    = fParameters[eToroidalZ];
     127              :   double& fTwistingTime = fParameters[eTwistingTime];
     128              : 
     129              :   // some pre-calculated derived parameter values
     130              :   double fSinPitch  = 0;
     131              :   double fCosPitch  = 0;
     132              :   double fTanPitch  = 0;
     133              : 
     134              :   /// major field components
     135              :   Vector3d getDiskField(const Vector3d& pos) const;
     136              :   Vector3d getHaloField(const Vector3d& pos) const;
     137              : 
     138              :   /// sub-components depending on model type
     139              :   /// -- Sec. 5.2.2
     140              :   Vector3d getSpiralField(const double x, const double y, const double z) const;
     141              :   /// -- Sec. 5.2.3
     142              :   Vector3d getSpurField(const double x, const double y, const double z) const;
     143              :   /// -- Sec. 5.3.1
     144              :   Vector3d getToroidalHaloField(const double x, const double y, const double z) const;
     145              :   /// -- Sec. 5.3.2
     146              :   Vector3d getPoloidalHaloField(const double x, const double y, const double z) const;
     147              :   /// -- Sec. 5.3.3
     148              :   Vector3d getTwistedHaloField(const double x, const double y, const double z) const;
     149              : 
     150              : };
     151              : /** @} */
     152              : } // namespace crpropa
     153              : #endif