ChargeInConstantMagneticField.cc

Go to the documentation of this file.

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

#include <utl/RK5ODEIntegrator.h>

#include <RelativisticChargeInMagneticFieldODE.h>

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


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


class ConstantMagneticField {
public:
  ConstantMagneticField(const Vector3& b0)
  { fB0[0] = b0[0]; fB0[1] = b0[1]; fB0[2] = b0[2]; }

  void
  operator()(const Vector6& /*position*/, Vector3& b)
  {
    // ignore position
    b[0] = fB0[0];
    b[1] = fB0[1];
    b[2] = fB0[2];
  }

private:
  Vector3 fB0;
};


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

  typedef ConstantMagneticField MFType;
  const Vector3 b0 = { 1, -1, 2 };
  MFType b(b0);

  typedef RelativisticChargeInMagneticFieldODE<MFType> ODE;
  ODE chargeMotionODE(1, 1, b);

  RK5ODEIntegrator<ODE> helix(chargeMotionODE);

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

  for (int i = 0; i < 30; ++i) {

    // 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, position and velocity
    cout << it.GetX();
    for (int j = 0; j < 6; ++j)
      cout << ' ' << xu[j];
    cout << '\n';

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

  return EXIT_SUCCESS;
}