ConversionUtil.cc

Go to the documentation of this file.

#include <utl/config.h>
#include "ConversionUtil.h"
#include "ErrorPropagation.h"

#include <utl/AugerException.h>
#include <utl/TimeStamp.h>
#include <utl/UTCDateTime.h>
#include <utl/ModifiedJulianDate.h>
#include <utl/NumericalErrorPropagation.h>
#include <utl/CorrelationMatrix.h>

#include <evt/Event.h>

#include <TVector3.h>

#include <RecEvent.h>
#include <RecShower.h>

#include <iostream>
#include <string>
#include <cmath>
#include <vector>

using namespace std;
using namespace otoa;
using namespace utl;


namespace otoa {

  unsigned int
  TimeStamp2YYMMDD(const utl::TimeStamp& timest)
  {
    const UTCDateTime t = timest;
    return (t.GetYear() % 100)*10000 + t.GetMonth()*100 + t.GetDay();
  }


  unsigned int
  TimeStamp2HHMMSS(const utl::TimeStamp& timest)
  {
    const UTCDateTime t = timest;
    return t.GetHour()*10000 + t.GetMinute()*100 + t.GetSecond();
  }


  double
  TimeStamp2MoonCycle(const utl::TimeStamp& timest)
  {
    const double kFirstFullMoon = 1795.157467;  // moon cycle at 2004/01/07
    const double kMoonCycle = 29.53;  // duration of one moon cycle in days
    return ModifiedJulianDate(timest) / kMoonCycle - kFirstFullMoon;
  }


  double
  TimeStamp2GMST(const utl::TimeStamp& ts)
  {
    // a la CDAS
    const UTCDateTime p = ts;
    const TimeStamp time0 = UTCDateTime(p.GetYear(), p.GetMonth(), p.GetDay()).GetTimeStamp();
    const double julianDay0 = ModifiedJulianDate(time0) + 2400000.5;
    const double tt = (julianDay0 - 2415020.0) / 36525.0;
    const double dgmst0 = (0.000001075*tt + 100.0021359)*tt + 0.276919398;
    const double gmst0 = (dgmst0 - int(dgmst0)) * 24.0;
    return gmst0 + (p.GetHour() + p.GetMinute()/60. + p.GetSecond()/3600.) * 366.25 / 365.25;
  }


  double
  CalculateNutationCorrection(const double julianDay)
  {
    // a la CDAS
    const  double tt = (julianDay - 2415020.0) / 36525.;  // time (Julius century) from 1900

    const double llt = (0.000303*tt + 36000.7689)*tt + 279.6967;
    const double ldt = (-0.001133*tt + 481267.8831)*tt + 270.4342;
    const double mmt = (-0.000150*tt + 35999.0498)*tt + 358.4758;
    const double mdt = (0.009192*tt + 477198.8491)*tt + 296.1046;
    const double omegat = (0.002078*tt - 1934.1420)*tt + 259.1833;

    const double ll = fmod(llt, 360.);  // mean ecliptic longtitude of the sun (deg)
    const double ld = fmod(ldt, 360.);  // mean ecliptic longtitude of the moon(deg)
    const double mm = fmod(mmt, 360.);  // mean anomaly of the sun  (deg)
    const double md = fmod(mdt, 360.);  // mean anomaly of the moon (deg)

    double omega = fmod(omegat, 360.);  // ecliptic lg. of ascending node of moon
    if (omega < 0)
      omega += 360;
    /* not used in the correction
       const double d_psi_t = -(17.2327 + 0.01737 * tt) * sin(deg*(omega))
       - (1.2729 + 0.00013 * tt) * sin(deg*(2.0 * ll))
       + 0.2088 * sin(deg*(2.0 * omega))
       - 0.2037 * sin(deg*(2.0 * ld))
       + (0.1261 - 0.00031 * tt) * sin(deg*(mm))
       + 0.0675 * sin(deg*(md))
       - (0.0497 - 0.00012 * tt) * sin(deg*(2.0 * ll + mm))
       - 0.0342 * sin(deg*(2.0 * ld - omega))
       - 0.0261 * sin(deg*(2.0 * ld + mm))
       + 0.0214 * sin(deg*(2.0 * ll - mm))
       - 0.0149 * sin(deg*(2.0 * (ll - ld) + md))
       + 0.0124 * sin(deg*(2.0 * ll - omega))
       + 0.0114 * sin(deg*(2.0 * ld - md));
    */
    const double depsit =
        (9.2100 + 0.00091*tt) * cos(omega*deg)
      + (0.5522 - 0.00029*tt) * cos(2*ll*deg)
      -  0.0904               * cos(2*omega*deg)
      +  0.0884               * cos(2*ld*deg)
      +  0.0216               * cos((2*ll + mm)*deg)
      +  0.0183               * cos((2*ld - omega)*deg)
      +  0.0113               * cos((2*ld + md)*deg)
      -  0.0093               * cos((2*ll - mm)*deg)
      -  0.0006               * cos((2*ll - omega)*deg);

    //vector<double> result;
    //result.push_back(d_psi_t / 3600.0);    //sol.d_psi
    //result.push_back( d_epsi_t / 3600.0); //sol.d_epsi

    /* central equation of the sun, not used in the correction
       const double cc = (1.919460 - 0.004789 * tt - 0.000014 * tt * tt) * sin(deg*(mm))
       + (0.020094 - 0.000100 * tt) * sin(deg*(2.0 * mm))
       + 0.000293 * sin(deg*(3.0 * mm));*/

    //result.push_back(ll + cc);    //sol.lambda_s

    const double epsi =
      ((0.000000503*tt - 0.00000164)*tt - 0.0130125)*tt + 23.452294 +
      0.00256*cos(omega*deg);  //sol.epsi
    //result.push_back(epsi*deg);    //sol.epsi
    double nuthour = depsit / 3600 * cos(epsi*deg);
    return nuthour;
  }


  void
  LinearFit(const vector<double>& x,
            const vector<double>& y,
            const vector<double>& ey,
            double& a0,
            double& a1,
            double& chi2)
  {
    double s1 = 0;
    double sx = 0;
    double sxx = 0;
    double sy = 0;
    double sxy = 0;
    double syy = 0;

    for (unsigned int i = 0; i < x.size(); ++i) {
      if (ey[i] > 0) {
        const double wi = 1 / (ey[i]*ey[i]);
        const double wX = wi * x[i];
        const double wY = wi * y[i];
        s1  += wi;
        sx  += wX;
        sxx += wX * x[i];
        sy  += wY;
        sxy += wX * y[i];
        syy += wY * y[i];
      }
    }

    const double d = s1*sxx - sx*sx;

    a0 = (sxx*sy - sx * sxy) / d;
    a1 = (-sx*sy + s1 * sxy) / d;

    chi2 = syy - a0*sy - a1*sxy;
  }


  void
  FillCelestialCoordinates(RecShower& recShower)
  {
    const TVector3& coreUTM = recShower.GetCoreUTMCS();
    TimeStamp coreTime(recShower.GetCoreTimeSecond(), recShower.GetCoreTimeNanoSecond());

    const double julianDay = ModifiedJulianDate(coreTime) + 2400000.5;
    const double gmst = TimeStamp2GMST(coreTime);

    // error on julian day and gmst negligible, relative error << 1e-5
    //   1e-5 if our clocks are wrong by second (time conversion GPS <-> UTC can cause this but very unlikely)
    //   1e-12 if we are wrong by some 100 nanoseconds (order of time offset between SD/FD)
    vector<Parameter> locParameters;
    locParameters.push_back(Parameter(julianDay, 0.));
    locParameters.push_back(Parameter(gmst, 0.));
    locParameters.push_back(Parameter(coreUTM.X(), recShower.GetCoreEastingError()));
    locParameters.push_back(Parameter(coreUTM.Y(), recShower.GetCoreNorthingError()));
    locParameters.push_back(Parameter(coreUTM.Z(), 0.));  // zero by definition
    locParameters.push_back(Parameter(recShower.GetZenith(), recShower.GetZenithError()));
    locParameters.push_back(Parameter(recShower.GetAzimuth(), recShower.GetAzimuthError()));

    CorrelationMatrix locCorrelations(7);
    // we only have these, others are set to zero automatically
    locCorrelations.Set(2, 3, recShower.GetCoreNorthingEastingCorrelation());
    locCorrelations.Set(5, 6, recShower.GetZenithAzimuthCorrelation());

    err::CalculateEquatorialCoordinates calcEqu;
    NumericalErrorPropagation pEqu(calcEqu, locParameters, locCorrelations);

    const double ra  = pEqu.GetFunctionResult()[0];
    const double dec = pEqu.GetFunctionResult()[1];

    // safeguards for NaNs
    const double raErr = std::isnan(pEqu.GetPropagatedErrors()[0]) ?
      ra : pEqu.GetPropagatedErrors()[0];
    const double decErr = std::isnan(pEqu.GetPropagatedErrors()[1]) ?
      dec : pEqu.GetPropagatedErrors()[1];
    const double raDecCorr = std::isnan(pEqu.GetPropagatedCorrelations().Get(0, 1)) ?
      0 : pEqu.GetPropagatedCorrelations().Get(0, 1);

    recShower.SetRightAscension(ra);
    recShower.SetRightAscensionError(raErr);
    recShower.SetDeclination(dec);
    recShower.SetDeclinationError(decErr);
    recShower.SetRDCorrelation(raDecCorr);

    // -------- galactic coordinates
    err::CalculateGalacticCoordinates calcGal;
    NumericalErrorPropagation pGal(calcGal, locParameters, locCorrelations);

    const double galLong = pGal.GetFunctionResult()[0];
    const double galLat  = pGal.GetFunctionResult()[1];

    // saveguards for NaNs
    const double galLongErr = std::isnan(pGal.GetPropagatedErrors()[0]) ?
      galLong : pGal.GetPropagatedErrors()[0];
    const double galLatErr = std::isnan(pGal.GetPropagatedErrors()[1]) ?
      galLat : pGal.GetPropagatedErrors()[1];
    const double galLatLongCorr = std::isnan(pGal.GetPropagatedCorrelations().Get(0, 1)) ?
      0 : pGal.GetPropagatedCorrelations().Get(0, 1);

    recShower.SetGalacticLatitude(galLat);
    recShower.SetGalacticLongitude(galLong);
    recShower.SetGalLatError(galLatErr);
    recShower.SetGalLongError(galLongErr);
    recShower.SetLongLatCorrelation(galLatLongCorr);
  }

}