MagneticFieldModel.cc

Go to the documentation of this file.

#include "MagneticFieldModel.h"

#include <utl/Vector.h>
#include <utl/Point.h>
#include <utl/CoordinateSystemPtr.h>
#include <utl/AugerUnits.h>
#include <utl/UTCDate.h>

#include <utl/AugerException.h>
#include <utl/Reader.h>

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

#include <string>

fwk::MagneticFieldModel::MagneticFieldModel()
{
  Init();
}

fwk::MagneticFieldModel& fwk::MagneticFieldModel::instance()
{
  static MagneticFieldModel f;
  return f;
}

void
fwk::MagneticFieldModel::Init()
{
  using namespace utl;
  using namespace std;
  using namespace fwk;
  Branch top =
    CentralConfig::GetInstance()->GetTopBranch("MagneticFieldModel");

  Branch models = top.GetChild("BaseModelList");

  for (Branch model = models.GetFirstChild(); model; model = model.GetNextSibling())
    {
      double year;
      model.GetChild("BaseYear").GetData(year);
      size_t order;
      model.GetChild("Order").GetData(order);
      vector<double> g;
      vector<double> h;
      model.GetChild("G").GetData(g);
      model.GetChild("H").GetData(h);

      fModels.push_back(ParametricGeoMagneticField(year,order,&g[0],&h[0]));
    }

  Branch changeModel =  top.GetChild("ChangeModel");
  {
    size_t order;
    changeModel.GetChild("Order").GetData(order);
    vector<double> g;
    vector<double> h;
    changeModel.GetChild("G").GetData(g);
    changeModel.GetChild("H").GetData(h);

    fChangeModel = new ParametricGeoMagChange(order,&g[0],&h[0]);
  }
}


double fwk::MagneticFieldModel::GetDeclination(const utl::Point &position,
                                               const utl::TimeStamp &time)
{
  utl::CoordinateSystemPtr cs = fwk::LocalCoordinateSystem::Create(position);

  return (90*utl::degree - GetMagneticFieldVector(position, time).GetPhi(cs));
}

double fwk::MagneticFieldModel::GetInclination(const utl::Point &position,
                                               const utl::TimeStamp &time)
{
  utl::CoordinateSystemPtr cs = fwk::LocalCoordinateSystem::Create(position);

  return (GetMagneticFieldVector(position, time).GetTheta(cs) - 90*utl::degree);
}


utl::Vector fwk::MagneticFieldModel::GetMagneticFieldVector(const utl::Point &position,
                                                            const utl::TimeStamp &time)
{
  return instance().ModelForTime(time).Get(position);
}


double
fwk::MagneticFieldModel::decimalYear(int year, int month, int day)
{
  int days[12] = { 0, 31, 59, 90, 120, 151, 181, 212, 243, 273, 304, 334};

  int leap_year = (((year % 4) == 0) &&  (((year % 100) != 0) || ((year % 400) == 0)));

  double day_in_year = (days[month - 1] + day + (month > 2 ? leap_year : 0));

  //did you know a day in a leap year is a little shorter ...
  return ((double) year + (day_in_year / (365.0 + leap_year)));
}


fwk::ParametricGeoMagneticField 
fwk::MagneticFieldModel::ModelForTime(const utl::UTCDate &date)
{
  int year = date.GetYear();
  int month = date.GetMonth();
  int day = date.GetDay();

  return ModelForYear(decimalYear(year, month, day));
}

fwk::ParametricGeoMagneticField 
fwk::MagneticFieldModel::ModelForYear(double year)
{
  double firstYear = fModels.front().GetBaseYear();
  if (year < firstYear)
    {
      std::ostringstream msg;
      msg << "No magnet field information before year " << year << "\n";
      ERROR(msg.str());
      abort();
      // There may be some sense in returning models.front(), but rather let user provide
      // older magnet field data or let her explicitely ask for the first year.
    }

  double lastYear = fModels.back().GetBaseYear();
  if (year > lastYear)
    {
      return fChangeModel -> Extrapolate (year - lastYear, fModels.back());
    }

  std::vector<ParametricGeoMagneticField>::iterator modelit = fModels.begin();
  do
    {
      ParametricGeoMagneticField &model = *modelit;
      ParametricGeoMagneticField &nextModel = *(++modelit);

      if (nextModel.GetBaseYear() >= year)
        {
          return model.Interpolate (year, nextModel);
        }
    }
  while (modelit != fModels.end());

  ERROR("This point should never be reached."); abort();
}