ChargeInHMR_ASS_A.cc

Go to the documentation of this file.

#include <iostream>
#include <cmath>
#include <cstdlib>

#include <utl/RK5ODEIntegrator.h>

#include <GalacticUnits.h>
#include <RelativisticChargeInMagneticFieldODE.h>
using namespace std;
#include <GMF_HMR.h>

using namespace std;
using namespace utl;
using namespace galactic;


typedef double Vector3[3];
typedef double Vector6[6];


template<typename Vector>
inline
double
RadialDistance(const Vector& v)
{
  const double r2 = v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
  return sqrt(r2);
}


int
main()
{
  typedef double Vector6[6];

  typedef gmf::HMR_BSS_S GMFType;
  //typedef gmf::HMR_ASS_A GMFType;
  GMFType hmr;

  typedef RelativisticChargeInMagneticFieldODE<GMFType> ODE;

  for (double energy = 0.01*EeV; energy <= 100*EeV; energy *= 1.3) {

    ODE chargeMotionODE(1*eplus, energy, hmr);

    RK5ODEIntegrator<ODE> integrator(chargeMotionODE);

    const double x0 = 0;
    const double dxFirst = 0.1*lightyear;
    // first 3 elements are position, last 3 are (unit) velocity
    const double sunGCDistance = 8.5*kiloparsec;
    const Vector6 x0u0 = { -sunGCDistance, 0, 50*parsec, 1, 0, 0 };
    const double accuracy = 1e-5;
    AdaptiveRK5Iterator<ODE, Vector6> it =
      integrator.AdaptiveBegin(x0, dxFirst, x0u0, accuracy);

    while (RadialDistance(it.GetY()) < 40*kiloparsec) {

      // normalize velocity vector to prevent numerical dissipation
      // this is optional but should be used for low accuracy runs
      Vector6& xu = it.GetY();
      /*const double uNorm =
        sqrt(xu[3]*xu[3] + xu[4]*xu[4] + xu[5]*xu[5]);
      xu[3] /= uNorm;
      xu[4] /= uNorm;
      xu[5] /= uNorm;*/

      // output current time (in units of ct = distance)
      cout << it.GetX() / kiloparsec;
      // output position (in kiloparsec)
      for (int j = 0; j < 3; ++j)
        cout << ' ' << xu[j] / kiloparsec;
      // output unit velocity
      for (int j = 3; j < 6; ++j)
        cout << ' ' << xu[j];
      cout << '\n';

      // advance with optimal step for this accuracy
      ++it;
    }

    cout << endl;

  }

  return EXIT_SUCCESS;
}