GMF_HMR.h

Go to the documentation of this file.

#ifndef _gal_GMF_HMR_ASS_A_h
#define _gal_GMF_HMR_ASS_A_h

#include <cmath>

#include <GalacticUnits.h>


namespace galactic {

  namespace gmf {


    namespace HMRConstants {

      const double kRange = 20 * kiloparsec;
      const double kL = 8.5 * kiloparsec;
      const double kRho1 = 2 * kiloparsec;
      const double kD = -0.5 * kiloparsec;
      const double kZ1 = 0.3 * kiloparsec;
      const double kZ2 = 4 * kiloparsec;
      const double kZ3 = 20 * parsec;
      const double kPitch = -170 * degree;
      const double kPsi = 10.55 * kiloparsec;
      const double kB0 = 3 * microgauss;

    }



    class HMR_ASS_A {
    public:
      template<typename PositionVector, typename FieldVector>
      void
      operator()(const PositionVector& position, FieldVector& b)
      {
        using namespace std;
        using namespace HMRConstants;

        const double& x = position[0];
        const double& y = position[1];
        const double& z = position[2];

        const double rho2 = x*x + y*y;
        const double r = sqrt(rho2 + z*z);

        b[2] = 0;

        // return with zero field if out of range
        if (r > kRange || !rho2) {
          b[0] = b[1] = 0;
          return;
        }

        const double rho = sqrt(rho2);

        const double b_rho =
          kB0 * kL *
          pow(tanh(rho / kRho1), 3) / rho2 *
          pow(cos(atan2(y, x) - log(rho / kPsi)/tan(kPitch) - M_PI), 2) *
          (0.5 / cosh(z / kZ1) + 0.5 / cosh(z / kZ2)) * tanh(z / kZ3);

        const double sinp = sin(kPitch);
        const double cosp = cos(kPitch);
        b[0] = b_rho * (sinp * x - cosp * y);
        b[1] = b_rho * (sinp * y + cosp * x);

        /*cout << "GMF " << b[0] << ' ' << b[1] << ' ' << b[2] << ' '
             << r/kiloparsec << ' ' << rho/kiloparsec << ' '
             << b_rho << endl;*/
      }

    };


    class HMR_BSS_S {
    public:
      template<typename PositionVector, typename FieldVector>
      void
      operator()(const PositionVector& position, FieldVector& b)
      {
        using namespace std;
        using namespace HMRConstants;

        const double& x = position[0];
        const double& y = position[1];
        const double& z = position[2];

        const double rho2 = x*x + y*y;
        const double r = sqrt(rho2 + z*z);

        b[2] = 0;

        // return with zero field if out of range
        if (r > kRange || !rho2) {
          b[0] = b[1] = 0;
          return;
        }

        const double rho = sqrt(rho2);

        const double b_rho =
          kB0 * kL *
          pow(tanh(rho / kRho1), 3) / rho2 *
          cos(atan2(y, x) - log(rho / kPsi)/tan(kPitch) - M_PI) *
          (0.5 / cosh(z / kZ1) + 0.5 / cosh(z / kZ2));

        const double sinp = sin(kPitch);
        const double cosp = cos(kPitch);
        b[0] = b_rho * (sinp * x - cosp * y);
        b[1] = b_rho * (sinp * y + cosp * x);
      }

    };

  }

}


#endif