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