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