Line data    Source code

       1             : #ifndef CRPROPA_GALACTICMAGNETICFIELD_H
       2             : #define CRPROPA_GALACTICMAGNETICFIELD_H
       3             : 
       4             : #include "crpropa/magneticField/MagneticField.h"
       5             : #include <cmath>
       6             : 
       7             : namespace crpropa {
       8             : /**
       9             :  * \addtogroup MagneticFields
      10             :  * @{
      11             :  */
      12             : 
      13             : /**
      14             :  @class ToroidalHaloField
      15             :  @brief Galactic halo field model from Prouza & Smida 2003 and Sun et al. 2008
      16             :  */
      17             : class ToroidalHaloField: public MagneticField {
      18             :         double b0; // halo field strength
      19             :         double z0; // vertical position
      20             :         double z1; // vertical scale
      21             :         double r0; // radial scale
      22             : 
      23             : public:
      24             :         /**
      25             :          * Constructor
      26             :          * @param b0 halo field strength
      27             :          * @param z0 vertical position
      28             :          * @param z1 vertical scale
      29             :          * @param r0 radial scale
      30             :         */
      31           1 :         ToroidalHaloField(double b0 = 1., double z0 = 1., double z1 = 1., double r0 = 1.) {
      32             :                 setParameters(b0, z0, z1, r0);
      33             :         }
      34             : 
      35             :         void setParameters(double b0, double z0, double z1, double r0) {
      36           1 :                 this->b0 = b0;
      37           1 :                 this->z0 = z0;
      38           1 :                 this->z1 = z1;
      39           1 :                 this->r0 = r0;
      40             :         }
      41             : 
      42           2 :         Vector3d getField(Vector3d pos) {
      43           2 :                 double r = sqrt(pos.x * pos.x + pos.y * pos.y) / r0; // in-plane radius in units of the radial scale
      44           2 :                 double b = b0 / (1 + pow((std::fabs(pos.z) - z0) / z1, 2)) * r * exp(1 - r);
      45           2 :                 double phi = pos.getPhi(); // azimuth
      46           2 :                 return Vector3d(cos(phi), sin(phi), 0) * b;
      47             :         }
      48             : };
      49             : /** @} */
      50             : 
      51             : /**
      52             :  * \addtogroup MagneticFields
      53             :  * @{
      54             :  */
      55             : 
      56             : /**
      57             :  @class LogarithmicSpiralField
      58             :  @brief Galactic disk field model of axisymmetric (ASS) or bisymmetric (BSS) logarithmic spiral shape
      59             :  */
      60             : class LogarithmicSpiralField: public MagneticField {
      61             : private:
      62             :         bool isBSS;   // true for BSS, false for ASS
      63             :         double b0;    // field strength
      64             :         double pitch; // pitch angle [rad]
      65             :         double rsol;  // distance of sun to galactic center
      66             :         double rc;    // radius of central region with constant field strength
      67             :         double d;     // distance to the first field reversal
      68             :         double z0;    // vertical attenuation length
      69             : 
      70             :         double phase; // phase of the spiral arms
      71             :         double cosPhase;
      72             :         double sinPitch;
      73             :         double cosPitch;
      74             :         double tanPitch;
      75             : 
      76           1 :         void updatePitchAngle() {
      77           1 :                 sinPitch = sin(pitch);
      78           1 :                 cosPitch = cos(pitch);
      79           1 :                 tanPitch = tan(pitch);
      80           1 :         }
      81             : 
      82           1 :         void updatePhase() {
      83           1 :                 phase = log(1 + d / rsol) / tanPitch - M_PI / 2;
      84           1 :                 cosPhase = cos(phase);
      85           1 :         }
      86             : 
      87             : public:
      88             :         /**
      89             :          * Constructor
      90             :          * @param isBSS switch for the magnetic field model
      91             :          *                              true for BSS, false for ASS
      92             :          * @param b0    magnetic field strength
      93             :          * @param pitch pitch angle [rad]
      94             :          * @param rsol  distance of sun from Galactic center
      95             :          * @param rc    radius of central region with constant field
      96             :          * @param d             distance to first field reversal
      97             :          * @param z0    vertical attenuation length
      98             :         */
      99             :         LogarithmicSpiralField(bool isBSS = true, double b0 = 1., double pitch = M_1_PI/4.,
     100           1 :                 double rsol = 8.5*kpc, double rc = 3*kpc, double d = 5*kpc, double z0 = 3*kpc) {
     101             :                 setParameters(isBSS, b0, pitch, rsol, rc, d, z0);
     102             :         }
     103             : 
     104             :         void setParameters(bool isBSS, double b0, double pitch, double rsol,
     105             :                         double rc, double d, double z0) {
     106           1 :                 this->isBSS = isBSS;
     107           1 :                 this->b0 = b0;
     108           1 :                 this->pitch = pitch;
     109           1 :                 this->rsol = rsol;
     110           1 :                 this->rc = rc;
     111           1 :                 this->d = d;
     112           1 :                 this->z0 = z0;
     113           1 :                 updatePitchAngle();
     114           1 :                 updatePhase();
     115             :         }
     116             : 
     117           1 :         Vector3d getField(Vector3d pos) const {
     118           1 :                 double r = sqrt(pos.x * pos.x + pos.y * pos.y); // in-plane radius
     119           1 :                 double b = b0 / cosPhase * rsol / std::max(r, rc);
     120             : 
     121           1 :                 double phi = pos.getPhi();
     122           1 :                 double c = cos(phi - log(r / rsol) / tanPitch + phase);
     123           1 :                 if (isBSS)
     124           1 :                         c = std::fabs(c);
     125           1 :                 b *= c * exp(-std::fabs(pos.z) / z0);
     126             : 
     127           1 :                 return Vector3d(cosPitch * cos(phi), sinPitch * sin(phi), 0) * b;
     128             :         }
     129             : };
     130             : /** @} */
     131             : 
     132             : }// namespace crpropa
     133             : 
     134             : #endif // CRPROPA_GALACTICMAGNETICFIELD_H