Line data    Source code

       1              : #include "crpropa/magneticField/PT11Field.h"
       2              : #include "crpropa/Units.h"
       3              : 
       4              : #include <algorithm>
       5              : 
       6              : namespace crpropa {
       7              : 
       8            0 : PT11Field::PT11Field() : useASS(true), useBSS(false), useHalo(true) {
       9              :         // disk parameters
      10            0 :         d = - 0.6 * kpc;
      11            0 :         R_sun = 8.5 * kpc;
      12            0 :         R_c = 5.0 * kpc;
      13            0 :         z0_D = 1.0 * kpc;
      14            0 :         B0_D = 2.0 * muG;
      15              : 
      16              :         // halo parameters
      17            0 :         z0_H = 1.3 * kpc;
      18            0 :         R0_H = 8.0 * kpc;
      19            0 :         B0_Hn = 4.0 * muG;
      20            0 :         B0_Hs = 4.0 * muG;
      21            0 :         z11_H = 0.25 * kpc;
      22            0 :         z12_H = 0.4 * kpc;
      23              : 
      24              :         // set ASS specific parameters
      25            0 :         setUseASS(true);
      26            0 : }
      27              : 
      28            0 : void PT11Field::SetParams() {
      29            0 :         cos_pitch = cos(pitch);
      30            0 :         sin_pitch = sin(pitch);
      31            0 :         PHI = cos_pitch / sin_pitch * log1p(d / R_sun) - M_PI / 2;
      32            0 :         cos_PHI = cos(PHI);
      33            0 : }
      34              : 
      35            0 : void PT11Field::setUseASS(bool use) {
      36            0 :         useASS = use;
      37            0 :         if (not(use))
      38              :                 return;
      39              : 
      40            0 :         if (useBSS) {
      41              :                 std::cout << "PT11Field: Disk field changed to ASS" << std::endl;
      42            0 :                 useBSS = false;
      43              :         }
      44              : 
      45            0 :         pitch = -5.0 * M_PI / 180;
      46            0 :         B0_Hs = 2.0 * muG;
      47            0 :         SetParams();
      48              : }
      49              : 
      50            0 : void PT11Field::setUseBSS(bool use) {
      51            0 :         useBSS = use;
      52            0 :         if (not(use))
      53              :                 return;
      54              : 
      55            0 :         if (useASS) {
      56              :                 std::cout << "PT11Field: Disk field changed to BSS" << std::endl;
      57            0 :                 useASS = false;
      58              :         }
      59              : 
      60            0 :         pitch = -6.0 * M_PI / 180;
      61            0 :         B0_Hs = 4.0 * muG;
      62            0 :         SetParams();
      63              : }
      64              : 
      65            0 : void PT11Field::setUseHalo(bool use) {
      66            0 :         useHalo = use;
      67            0 : }
      68              : 
      69            0 : bool PT11Field::isUsingASS() {
      70            0 :         return useASS;
      71              : }
      72              : 
      73            0 : bool PT11Field::isUsingBSS() {
      74            0 :         return useBSS;
      75              : }
      76              : 
      77            0 : bool PT11Field::isUsingHalo() {
      78            0 :         return useHalo;
      79              : }
      80              : 
      81            0 : Vector3d PT11Field::getField(const Vector3d& pos) const {
      82            0 :         double r = sqrt(pos.x * pos.x + pos.y * pos.y);  // in-plane radius
      83              : 
      84              :         Vector3d b(0.);
      85              : 
      86              :         // disk field
      87            0 :         if ((useASS) or (useBSS)) {
      88              :                 // PT11 paper has B_theta = B * cos(p) but this seems because they define azimuth clockwise, while we have anticlockwise.
      89              :                 // see Tinyakov 2002 APh 18,165: "local field points to l=90+p" so p=-5 deg gives l=85 and hence clockwise from above.
      90              :                 // so to get local B clockwise in our system, need minus (like Sun etal).
      91              :                 // Ps base their system on Han and Qiao 1994 A&A 288,759 which has a diagram with azimuth clockwise, hence confirmed.
      92              : 
      93              :                 // PT11 paper define Earth position at (+8.5, 0, 0) kpc; but usual convention is (-8.5, 0, 0)
      94              :                 // thus we have to rotate our position by 180 degree in azimuth
      95            0 :                 double theta = M_PI - pos.getPhi();  // azimuth angle theta: PT11 paper uses opposite convention for azimuth
      96              :                 // the following is equivalent to sin(pi - phi) and cos(pi - phi) which is computationally slower
      97            0 :                 double cos_theta = - pos.x / r;
      98            0 :                 double sin_theta = pos.y / r;
      99              : 
     100              :                 // After some geometry calculations (on whiteboard) one finds:
     101              :                 // Bx = +cos(theta) * B_r - sin(theta) * B_{theta}
     102              :                 // By = -sin(theta) * B_r - cos(theta) * B_{theta}
     103              :                 // Use from paper: B_theta = B * cos(pitch)     and B_r = B * sin(pitch)
     104            0 :                 b.x = sin_pitch * cos_theta - cos_pitch * sin_theta;
     105            0 :                 b.y = - sin_pitch * sin_theta - cos_pitch * cos_theta;
     106              :                 b *= -1;        // flip magnetic field direction, as B_{theta} and B_{phi} refering to 180 degree rotated field
     107              : 
     108            0 :                 double bMag = cos(theta - cos_pitch / sin_pitch * log(r / R_sun) + PHI);
     109            0 :                 if (useASS)
     110            0 :                         bMag = fabs(bMag);
     111            0 :                 bMag *= B0_D * R_sun / std::max(r, R_c) / cos_PHI * exp(-fabs(pos.z) / z0_D);
     112              :                 b *= bMag;
     113              :         }
     114              : 
     115              :         // halo field
     116            0 :         if (useHalo) {
     117            0 :                 double bMag = (pos.z > 0 ? B0_Hn : - B0_Hs);
     118            0 :                 double z1 = (fabs(pos.z) < z0_H ? z11_H : z12_H);
     119            0 :                 bMag *= r / R0_H * exp(1 - r / R0_H) / (1 + pow((fabs(pos.z) - z0_H) / z1, 2.));
     120              :                 // equation (8) in paper: theta uses now the conventional azimuth definition in contrast to equation (3)
     121              :                 // cos(phi) = pos.x / r (phi going counter-clockwise)
     122              :                 // sin(phi) = pos.y / r
     123              :                 // unitvector of phi in polar coordinates: (-sin(phi), cos(phi), 0)
     124            0 :                 b += bMag * Vector3d(-pos.y / r, pos.x / r, 0);
     125              :         }
     126              : 
     127            0 :         return b;
     128              : }
     129              : 
     130              : } // namespace crpropa