ParametricGeoMagneticField.cc

Go to the documentation of this file.

#include "ParametricGeoMagneticField.h"

#include <utl/AugerUnits.h>
#include <utl/SVector.h>
#include <utl/Vector.h>
#include <utl/Point.h>

#include <fwk/LocalCoordinateSystem.h>
#include <fwk/CoordinateSystemRegistry.h>

#include <iostream>


namespace fwk {

  ParametricGeoMagneticField::ParametricGeoMagneticField(const double baseyear, const size_t order,
                                                         const double g[], const double h[]) :
    fBaseYear(baseyear),
    fOrder(order),
    fG(std::valarray<double>(g, order*(order+1)/2)),
    fH(std::valarray<double>(h, order*(order+1)/2))
  { }


  utl::SVector<3, double>
  ParametricGeoMagneticField::GetXYZ(const double latitude, const double longitude, const double elevation)
    const
  {
    // Implementation follows the reference implementation
    // of the [IGRF](http://www.ngdc.noaa.gov/IAGA/vmod/igrf.html)

    using namespace utl;
    double sl[fOrder];
    double cl[fOrder];
    {
      const double sin_longitude = sl[0] = sin(longitude);
      const double cos_longitude = cl[0] = cos(longitude);
      for (size_t i = 1; i < fOrder; ++i) {
        sl[i] = sl[i-1] * cos_longitude + cl[i-1] * sin_longitude;
        cl[i] = cl[i-1] * cos_longitude - sl[i-1] * sin_longitude;
      }
    }

    const size_t npq = (fOrder * (fOrder + 3)) / 2;
    double p[npq];
    double q[npq];

    const double sin_latitude = sin(latitude);
    const double cos_latitude = cos(latitude);
    p[0] = 2 * sin_latitude;
    p[1] = 2 * cos_latitude;
    p[2] = 4.5 * sin_latitude * sin_latitude - 1.5;
    p[3] = 3 * sqrt(3) * cos_latitude * sin_latitude;
    q[0] = -cos_latitude;
    q[1] = sin_latitude;
    q[2] = -3 * cos_latitude * sin_latitude;
    q[3] = sqrt(3) * (sin_latitude*sin_latitude - cos_latitude * cos_latitude);

    for (size_t k = 4, n = 2, m = 2; k < npq; ++k, ++m) {
      if (n < m) {
        m = 0;
        ++n;
      }
      if (m == n) {
        const double aa = sqrt(1 - 0.5/n);
        const size_t j = k - n - 1;
        p[k] =  aa * (1 + 1./n) * cos_latitude * p[j];
        q[k] =  aa * (cos_latitude * q[j] + sin_latitude/n * p[j]);
      } else {
        const double fn = n;
        const double fm = m;
        const double aa = sqrt (fn*fn - fm*fm);
        const double bb = sqrt ((fn - 1)*(fn - 1) - (fm*fm)) / aa;
        const double cc = (2*fn - 1) / aa;
        const size_t i = k - n;
        const size_t j = k - 2*n + 1;
        p[k] = (fn + 1) * (cc * sin_latitude / fn * p[i] - bb/(fn - 1) * p[j]);
        q[k] = cc * (sin_latitude * q[i] - cos_latitude / fn * p[i]) - bb * q[j];
      }
    }

    double x = 0;
    double y = 0;
    double z = 0;

    const double earth_radius = 6371.2 * utl::kilometer;
    const double ratio = earth_radius / elevation;

    for (size_t n = 1, k = 0; n < fOrder; ++n) {
      const double rr = pow(ratio, n + 2);
      const double aa = rr * fG[k];

      x += aa * q[k];
      z -= aa * p[k];

      ++k;
      for (size_t m = 0; m < n; ++m, ++k) {
        const double aa = rr * fG[k];
        const double bb = rr * fH[k];
        const double cc = aa * cl[m] + bb * sl[m];

        x += cc * q[k];
        z -= cc * p[k];

        if (cos_latitude > 0) {
          const double fm = m + 1;
          const double fn = n;
          y +=  (aa * sl[m] - bb * cl[m]) * fm * p[k] / ((fn + 1) * cos_latitude);
        } else
          y +=  (aa * sl[m] - bb * cl[m]) * q[k] * sin_latitude;
      }
    }

    const double r[] = { x*nano*tesla, y*nano*tesla, z*nano*tesla };
    const utl::SVector<3, double> result(r);
    return result;
  }


  double 
  ParametricGeoMagneticField::GetDeclination(const utl::Point& where)
    const
  {
    const utl::CoordinateSystemPtr cs = fwk::LocalCoordinateSystem::Create(where);
    return Get(where).GetPhi(cs) - 90*utl::degree; // magnetic field declination in WEST direction
  }


  utl::Vector
  ParametricGeoMagneticField::Get(const utl::Point& where)
    const
  {
    const utl::CoordinateSystemPtr earthCS = fwk::CoordinateSystemRegistry::Get("ECEF");

    // alas ...
    // double longitude = where.GetTheta(earthCS);
    // double latitude = where.GetPhi(earthCS);
    const double x = where.GetX(earthCS);
    const double y = where.GetY(earthCS);
    const double z = where.GetZ(earthCS);

    const double longitude = atan2(y, x);
    const double rho = where.GetRho(earthCS); // sqrt(x*x + y*y);
    const double latitude = atan2(z, rho);
    const double height = where.GetR(earthCS); // sqrt(x*x + y*y + z*z);

    const utl::SVector<3, double> xyz = GetXYZ(latitude, longitude, height);

    const utl::CoordinateSystemPtr localCS = LocalCoordinateSystem::Create(where);
    return utl::Vector(xyz[1], xyz[0], -xyz[2], localCS);
  }


  void
  ParametricGeoMagneticField::Grow(const size_t order)
  {
    if (order < fOrder)
      WARNING("Forgetting coefficients by shrinking model. You asked for it.");

    fOrder = order;
    fG.resize(order*(order+1) / 2, 0.);
    fH.resize(order*(order+1) / 2, 0.);
  }


  ParametricGeoMagneticField
  ParametricGeoMagneticField::Interpolate(const double year,
                                          const ParametricGeoMagneticField& later)
    const
  {
    assert(year > fBaseYear);
    assert(year < later.GetBaseYear());

    const double factor = (year - fBaseYear) / (later.GetBaseYear() - fBaseYear);

    ParametricGeoMagneticField model = *this; // copy

    if (later.fOrder > model.fOrder)
      model.Grow(later.fOrder);

    if (model.fOrder == later.fOrder) {
      model.fG += factor * (later.fG - model.fG);
      model.fH += factor * (later.fH - model.fH);
      return model;
    }

    // unlikely and maybe senseless case
    ERROR("Order adaption of shrinking models unimplemented (laziness failure)");
    abort();
  }


  /* Meant for testing
     Parameters should rather be read from file/db.
  */
  ParametricGeoMagneticField
  ParametricGeoMagneticField::IGRF2010()
  {
    double g[] = {
      -29496.50,
      -1585.90,
      -2396.60,
      3026.00,
      1668.60,
      1339.70,
      -2326.30,
      1231.70,
      634.20,
      912.60,
      809.00,
      166.60,
      -357.10,
      89.70,
      -231.10,
      357.20,
      200.30,
      -141.20,
      -163.10,
      -7.70,
      72.80,
      68.60,
      76.00,
      -141.40,
      -22.90,
      13.10,
      -77.90,
      80.40,
      -75.00,
      -4.70,
      45.30,
      14.00,
      10.40,
      1.60,
      4.90,
      24.30,
      8.20,
      -14.50,
      -5.70,
      -19.30,
      11.60,
      10.90,
      -14.10,
      -3.70,
      5.40,
      9.40,
      3.40,
      -5.30,
      3.10,
      -12.40,
      -0.80,
      8.40,
      -8.40,
      -10.10,
      -2.00,
      -6.30,
      0.90,
      -1.10,
      -0.20,
      2.50,
      -0.30,
      2.20,
      3.10,
      -1.00,
      -2.80,
      3.00,
      -1.50,
      -2.10,
      1.60,
      -0.50,
      0.50,
      -0.80,
      0.40,
      1.80,
      0.20,
      0.80,
      3.80,
      -2.10,
      -0.20,
      0.30,
      1.00,
      -0.70,
      0.90,
      -0.10,
      0.50,
      -0.40,
      -0.40,
      0.20,
      -0.80,
      0.00,
      -0.20,
      -0.90,
      0.30,
      0.40,
      -0.40,
      1.10,
      -0.30,
      0.80,
      -0.20,
      0.40,
      0.00,
      0.40,
      -0.30,
      -0.30
    };

    double h[] = {
      0,
      4945.10,
      0.00,
      -2707.70,
      -575.40,
      0.00,
      -160.50,
      251.70,
      -536.80,
      0.00,
      286.40,
      -211.20,
      164.40,
      -309.20,
      0.00,
      44.70,
      188.90,
      -118.10,
      0.10,
      100.90,
      0.00,
      -20.80,
      44.20,
      61.50,
      -66.30,
      3.10,
      54.90,
      0.00,
      -57.80,
      -21.20,
      6.60,
      24.90,
      7.00,
      -27.70,
      -3.40,
      0.00,
      10.90,
      -20.00,
      11.90,
      -17.40,
      16.70,
      7.10,
      -10.80,
      1.70,
      0.00,
      -20.50,
      11.60,
      12.80,
      -7.20,
      -7.40,
      8.00,
      2.20,
      -6.10,
      7.00,
      0.00,
      2.80,
      -0.10,
      4.70,
      4.40,
      -7.20,
      -1.00,
      -4.00,
      -2.00,
      -2.00,
      -8.30,
      0.00,
      0.10,
      1.70,
      -0.60,
      -1.80,
      0.90,
      -0.40,
      -2.50,
      -1.30,
      -2.10,
      -1.90,
      -1.80,
      0.00,
      -0.80,
      0.30,
      2.20,
      -2.50,
      0.50,
      0.60,
      0.00,
      0.10,
      0.30,
      -0.90,
      -0.20,
      0.80,
      0.00,
      -0.80,
      0.30,
      1.70,
      -0.60,
      -1.20,
      -0.10,
      0.50,
      0.10,
      0.50,
      0.40,
      -0.20,
      -0.50,
      -0.80
    };

    return ParametricGeoMagneticField(2010, 14, g, h);
  }

}