// tsid.cc
// Eric Aubourg         CEA/DAPNIA/SPP   aout 1999
// adapted from TSid.java, E. Aubourg 1998.

#include "tsid.h"

const double TSid::SIDRATE  = 0.9972695677;
const double TSid::MJDRef   = 1545;      // 1/1/2000, 12h UT1

TSid::TSid(double lg) 
: longitude(lg)
{
  InitDUT();
}

void TSid::InitDUT() {
  dut1[1270.5] = 0.60846; // 1999,3,2 (april, 2 1999)
  dut1[1299.5] = 0.57975; // 1999,4,1
  dut1[1330.5] = 0.55305; // 1999,5,1
  dut1[1360.5] = 0.54229; // 1999,6,1
  dut1[1391.5] = 0.53916; // 1999,7,1
  dut1[1422.5] = 0.52822; // 1999,8,1
  dut1[1452.5] = 0.50279; // 1999,9,1
  dut1[1483.5] = 0.46517; // 1999,10,1
  dut1[1513.5] = 0.42700; // 1999,11,1
  dut1[1544.5] = 0.39116; // 2000,0,1
  dut1[1575.5] = 0.35876; // 2000,1,1
  dut1[1604.5] = 0.30979; // 2000,2,1
  dut1[1635.5] = 0.25614; // 2000,3,1
}

double TSid::getDUT1(double mjd) {
  map<double, double>::iterator i = dut1.upper_bound(mjd);
  if (i==dut1.end() || i==dut1.begin()) return -1;
  map<double, double>::iterator j = i; j--;
  double x = (mjd - (*j).first) / ((*i).first - (*j).first);
  return (1-x)*(*j).second + x * (*i).second;
}

void TSid::setLongitude(double lg) {
  longitude = lg;
}

double TSid::getLST(double mjd) {
  double dut1 = getDUT1(mjd);
  double T = ((mjd-MJDRef) + dut1/86400.)/36525.;
    // T = siecles julien depuis 1/1/2000, 12h UT1
  double GMST = 24110.54841 + 8640184.812866 * T +
     0.093104 * T * T - 0.0000062 * T * T * T; // a 0h UT1
  GMST += ((mjd-.5 - (int)(mjd-.5))*86400. + dut1);
  double LST = GMST + (longitude/15. * 3600.);
  LST -= 86400 * int(LST/86400);
  if (LST<0) LST += 86400.;
  return LST;
}