NcTimestamp.cxx

Go to the documentation of this file.

//    by the introduction of Leap Seconds when abs(UT1-UTC) exceeds 0.9 sec.
//    This time synchronisation is coordinated by the International Earth Rotation and
//    Reference Systems Service (IERS) via a daily monitoring of the Earth Orientation
//    Parameters (EOP). 
//    Depending on dUT=UT1-UTC, these Leap Seconds can be positive or negative.
//    The introduction of Leap Seconds into UTC started at 01-jan-1972 00:00:00 UT.
//
//    An overview of the history of introduced Leap seconds (TAI-UTC) is online available at :
//                 https://hpiers.obspm.fr/iers/bul/bulc/Leap_Second.dat
//                     or  http://maia.usno.navy.mil/ser7/tai-utc.dat
//
//    The time difference dUT=UT1-UTC is monitored on a daily basis and the data are available at :
//                 https://hpiers.obspm.fr/iers/series/opa/eopc04
//                     or  http://maia.usno.navy.mil/ser7/ser7.dat
//
//    The accuracy of dUT=UT1-UTC is about 10 microseconds, so in case of accurate astronomical
//    timing the user is advised to specify dUT in the various Set() facilities or use SetUT().
//    An automatic setting of dUT is provided based on the IERS data files. 
//    Please refer to the member function LoadUTCparameterFiles() for further details.
//
// 4) In some cases Unix Time (also called POSIX Time or UNIX Epoch Time) is used.
//    Unix Time is closely releated to UTC and represents the (fractional) elapsed second count
//    since the start of the Unix Epoch, in which every Unix Time day contains exactly 86400 seconds.
//    This implies that for Unix Time the UTC leap seconds should be ignored, as explained below.
//    The UNIX Time EPOCH starts at 01-jan-1970 00:00:00 UTC which corresponds to JD=2440587.5
//    (i.e. the start of MJD=40587, as explained in the discussion on Julian Date below).
//    Synchronization with UTC is obtained by continuing the second count when a UTC leap second
//    occurs and then the Unix Time count jumps up or down 1 second at the beginning of the new
//    day after the UTC leap second has expired.
//    In case of a negative UTC leap second (which has not occurred so far) this would result
//    in a jump forward of 1 second in Unix Time, introducing a "gap" of 1 second in the timing
//    at the start of the new day.
//    In case of a positive UTC leap second this results in a jump back of 1 second in Unix Time
//    at the start of a new day.
//    In this case there exist 2 ambiguous Unix Times over a 1 second time interval,
//    namely at the beginning of the introduction of the UTC leap second and at the moment that
//    the UTC leap second expires, i.e. the start of the new day.
//    This NcTimestamp facility provides setting and retrieval of Unix Time, but for accurate
//    timing the user is advised to use one of the time scales mentioned below. 
// 
// Supported time scales :
// ----------------------
// UT  : Universal Time.
//       1 day is the time span between two successive solar meridian transitions.
//       This implies that UT is based on the actual rotation of the Earth and as such suited
//       for (synchronization of) astronomical observations at various locations on Earth.
//       The reference time is defined by transitions over the meridian at Greenwich.
// ST  : Siderial Time (see also the details below).
//       1 day is the time span between two successive stellar meridian transitions.
//       This implies that ST is based on the actual rotation of the Earth and its rotation
//       around the Sun. Siderial Time is very well suited for astronomical observations.
//       The reference time is defined by transitions over the meridian at Greenwich.
//       Because of the stability of the rotation of the Earth around the Sun, the Siderial Time
//       is derived from the Universal Time and as such not treated as a separate time scale.
// TAI : International Atomic Time (Temps Atomique International).
//       1 day is the time span of 86400 standard atomic (SI) seconds at average sea level.
//       The standard atomic (SI) second is derived from a set of Cs atomic clocks.
//       This implies that TAI is not related to the actual rotation of the Earth and as such
//       will not run at the same pace as UT.
//       The start epoch of TAI is 01-jan-1958 00:00:00 UT, at which time UT1-UTC was about 0.
// GPS : Global Positioning System time.
//       This satellite based timing system is based on TAI and broadcast via a satellite network.
//       The start epoch of GPS is 06-jan-1980 00:00:00 UTC, at which the number of Leap Seconds was 19.
//       This implies that at that time TAI=UTC+19 sec. and consequently we always have TAI=GPS+19 sec.
//       GPS is broadcast in two formats : (w,sow) and (w,dow,sod) with
//       w   = week number after the start epoch.   1 week=7 days or 1 week=604800 seconds.
//       sow = (fractional) second count in the current week
//       dow = day number in the current week (sunday=0, monday=1 etc.)
//       sod = (fractional) second count in the current day
//       Early implementations used to reset the week count after a cycle of 1024 weeks and
//       provided the corresponding cycle count. Please refer to the SetGPS() memberfunctions
//       for further details.
// TT  : Terrestrial Time.
//       This timing system is based on TAI and provides the date/time at average sea level.
//       It has been introduced for observations from the surface of the Earth, to be consistent
//       with General Relativity. Since planetary orbits are very stable in time and not related
//       to the rotation of the Earth, TT is mainly used for observations of the Solar system.
//       TT was synchronized with TAI at 01-jan-1977 00:00:00 TAI, to indicate 00:00:32.184 in
//       order to provide a continuation of its (now obsolete) predecessor Ephemeris Time (ET).
//       This implies that always TT=TAI+32.184 sec.
//
// All Gregorian dates/times (e.g. 12-aug-1982 13:20:12) are treated on basis of
// the time scales mentioned above. So, there is no need for leap second treatment
// like this is the case for Coordinated Universal Time (UTC).  
//
// In order to enable a precise measurement of elapsed time over (very) long periods,
// a continuous day counting system has been introduced, called the Julian Date.
// The Julian Date (JD) indicates the number of days since noon (12:00:00) on
// 01 jan -4712 (i.e. noon 01 jan 4713 BC), being day 0 of the Julian calendar.
// It is custom to couple the Julian Date to UT so that it serves astronomical observations.
// However, this NcTimestamp facility allows to use a continuous day counting system
// with all the supported time scales mentioned above.
//
// The Modified Julian Date (MJD) indicates the number of days since midnight (00:00:00)
// on 17-nov-1858, which corresponds to 2400000.5 days after day 0 of the
// Julian calendar.
//
// The Truncated Julian Date (TJD) corresponds to 2440000.5 days after day 0
// of the Julian calendar and consequently TJD=MJD-40000.
// This TJD date indication was used by the Vela and CGRO satellite missions in
// view of Gamma Ray Burst investigations.
//
// The Julian Epoch (JE) indicates the fractional elapsed Julian year count
// since the start of the Gregorian year count.
// A Julian year is defined to be 365.25 days and starts at 01-jan 12:00:00.
// As such, the integer part of JE corresponds to the usual Gregorian year count,
// apart from 01-jan before 12:00:00.
// So, 01-jan-1965 12:00:00 UT corresponds to JE=1965.0
//
// The Besselian Epoch (BE) indicates the fractional elapsed Besselian year count
// since the start of the Gregorian year count.
// A Besselian (or tropical) year is defined to be 365.242198781 days.
// The date 31-dec-1949 22:09:46.862 UT corresponds to BE=1950.0
//
// The Besselian and Julian epochs are used in astronomical catalogs
// to denote values of time varying observables like e.g. right ascension.
//
// Because of the fact that the Julian date indicators are all w.r.t. UT or TAI derived
// time scales, they provide an absolute time scale irrespective of timezone or daylight
// saving time (DST).
//
// In view of astronomical observations and positioning it is convenient
// to have also a UT equivalent related to stellar meridian transitions.
// This is achieved by the Greenwich Sidereal Time (GST).
// The GST is defined as the right ascension of the objects passing
// the meridian over Greenwich.
// Due to the rotation of the Earth around the Sun, a sidereal day
// lasts 86164.09 seconds (23h 56m 04.09s) compared to the mean solar
// day of 86400 seconds (24h).
// Furthermore, precession of the earth's spin axis results in the fact
// that the zero point of right ascension (vernal equinox) gradually 
// moves along the celestial equator.
// In addition, tidal friction and ocean and atmospheric effects will
// induce seasonal variations in the earth's spin rate and polar motion
// of the earth's spin axis.
// Taking the above effects into account leads to what is called
// the Greenwich Mean Sidereal Time (GMST).
// In case also the nutation of the earth's spin axis is taken into
// account we speak of the Greenwich Apparent Sidereal Time (GAST).
//
// This NcTimestamp facility allows for picosecond precision, in view
// of time of flight analyses for particle physics experiments.
// For normal date/time indication the standard nanosecond precision
// will in general be sufficient.
// Picosecond precision can be obtained by invokation of GetPs() or GetDifference().
// Note that when the fractional JD, MJD and TJD counts are used instead
// of the integer (days,sec,ns) specification, the nanosecond precision
// may be lost due to computer accuracy w.r.t. floating point operations.
//
// The TTimeStamp EPOCH starts at 01-jan-1970 00:00:00
// which corresponds to JD=2440587.5 or the start of MJD=40587 or TJD=587.
// Using the corresponding MJD of this EPOCH allows construction of
// the yy-mm-dd hh:mm:ss:ns TTimeStamp from a given input (M/T)JD and time.
// Obviously this TTimeStamp implementation would prevent usage of values
// smaller than JD=2440587.5 or MJD=40587 or TJD=587.
// Furthermore, due to a limitation on the "seconds since the EPOCH start" count
// in TTimeStamp, the latest accessible date/time is until 19-jan-2038 02:14:07.
// However, this NcTimestamp facility provides support for the full range
// of (M/T)JD values, but the setting of the corresponding TTimeStamp parameters
// is restricted to the values allowed by the TTimeStamp implementation.
// For these earlier/later (M/T)JD values, the standard TTimeStamp parameters will
// be set corresponding to the start of the TTimeStamp EPOCH.
// This implies that for these earlier/later (M/T)JD values the TTimeStamp parameters
// do not match the Julian parameters of NcTimestamp.
// As such the standard TTimeStamp parameters do not appear on the print output
// when invoking the Date() memberfunction for these earlier/later (M/T)JD values.  
//
// Examples :
// ==========
//
// Note : All TTimeStamp functionality is available as well.
//
// NcTimestamp t;
//
// t.Date();
//
// // Set a specific Date/Time in Universal Time (UT1)
// // without recording UTC nor the corresponding International Atomic Time (TAI)
// t.SetUT("22-08-2016","14:03:29.7358",0,"U");
//
// // Set a specific Date/Time in Universal Time (UTC)
// // without recording UT1 nor the corresponding International Atomic Time (TAI)
// t.SetUT("22-08-2016","14:03:29.7358",0,"N");
//
// // Set a specific Date/Time in Universal Time (UTC)
// // and also record UT1 and the corresponding International Atomic Time (TAI)
// // using manual "leap" and "dut" settings
// Int_t leap=34;
// Double_t dut=-0.14;
// t.SetUT("17-05-2011","05:43:18.31468",0,"M",leap,dut);
//
// // Set a specific Date/Time in Universal Time (UTC)
// // and also record UT1 and the corresponding International Atomic Time (TAI)
// // using automatic "leap" and "dut" settings from the IERS data files
// t.LoadUTCparameterFiles("leap.txt","dut.txt");
// t.SetUT("17-05-2011","05:43:18.31468",0,"A");
//
// // Set a specific Date/Time in Global Positioning System time (GPS)
// // and also record the corresponding Universal Times (UTC and UT1) and TAI
// // using automatic "leap" and "dut" settings from the IERS data files
// t.LoadUTCparameterFiles("leap.txt","dut.txt");
// t.SetTAI("GPS","17-05-2011","05:43:18.31468",0,"A",0,0);
// 
// // Retrieve Julian Date to ns precision
// Int_t jd,jsec,jns;
// t.GetJD(jd,jsec,jns);
// // Get the remaining ps precision
// Int_t ps=GetPs();
//
// // Retrieve fractional Truncated Julian Date
// Double_t tjd=t.GetTJD();
//
// // Retrieve fractional Julian Epoch
// Double_t je=t.GetJE();
//
// // Set to a specific Modified Julian Date representing UTC
// // without recording the corresponding International Atomic Time (TAI)
// Int_t mjd=50537;
// Int_t mjsec=1528;
// Int_t mjns=185643;
// Int_t mjps=35;
// t.SetMJD(mjd,mjsec,mjns,mjps,"N");
//
// t.Date();
//
// // Set to a specific Modified Julian Date representing UTC
// // and also record the corresponding International Atomic Time (TAI)
// // using automatic "leap" and "dut" settings from the IERS data files
// Int_t mjd=58457;
// Int_t mjsec=1528;
// Int_t mjns=185643;
// Int_t mjps=35;
// t.LoadUTCparameterFiles("leap.txt","dut.txt");
// t.SetMJD(mjd,mjsec,mjns,mjps,"A");
//
// t.Date();
//
// // Set to a specific Modified Julian Date representing UT1
// // and also record the corresponding International Atomic Time (TAI)
// // using automatic "leap" and "dut" settings from the IERS data files
// Int_t mjd=58457;
// Int_t mjsec=1528;
// Int_t mjns=185643;
// Int_t mjps=35;
// t.LoadUTCparameterFiles("leap.txt","dut.txt");
// t.SetMJD(mjd,mjsec,mjns,mjps,"U");
//
// t.Date();
//
// // Time intervals for e.g. Trigger or Time Of Flight analysis
// NcEvent evt;
// NcTrack* tx=evt.GetTrack(5);
// NcTimestamp* timex=tx->GetTimestamp();
// Double_t dt=evt.GetDifference(timex,"ps");
// NcTimestamp trig((NcTimestamp)evt);
// trig.Add(0,0,2,173);
// NcSignal* sx=evt.GetHit(23);
// NcTimestamp* timex=sx->GetTimestamp();
// Double_t dt=trig.GetDifference(timex,"ps");
// Int_t d,s,ns,ps;
// trig.GetDifference(timex,d,s,ns,ps);
//
// // Some practical conversion facilities
// // Note : They don't influence the actual date/time settings
// //        and as such can also be invoked as NcTimestamp::Convert(...) etc...
// Int_t y=1921;
// Int_t m=7;
// Int_t d=21;
// Int_t hh=15;
// Int_t mm=23;
// Int_t ss=47;
// Int_t ns=811743;
// Double_t jdate=t.GetJD(y,m,d,hh,mm,ss,ns);
//
// Int_t days,secs,nsecs;
// Double_t date=421.1949327;
// t.Convert(date,days,secs,nsecs);
//
// days=875;
// secs=23;
// nsecs=9118483;
// date=t.Convert(days,secs,nsecs);
//
// Double_t mjdate=40563.823744;
// Double_t epoch=t.GetJE(mjdate,"mjd");
//
//--- Author: Nick van Eijndhoven 28-jan-2005 Utrecht University
//- Modified: Nick van Eijndhoven, IIHE-VUB Brussel, UTC June 28, 2026  11:10
~~~
**/
 
#include "NcTimestamp.h"
#include "Riostream.h"
 
ClassImp(NcTimestamp); // Class implementation to enable ROOT I/O
 
NcTimestamp::NcTimestamp() : TTimeStamp()
{
 
 FillJulian();
 fJps=0;
 fUtc=0;
 fLeap=0;
 fDut=0;
 fTmjd=0;
 fTsec=0;
 fTns=0;
 fTps=0;
 fUTCdata=0;
}

NcTimestamp::NcTimestamp(TTimeStamp& t) : TTimeStamp(t)
{
 
 FillJulian();
 fJps=0;
 fUtc=0;
 fLeap=0;
 fDut=0;
 fTmjd=0;
 fTsec=0;
 fTns=0;
 fTps=0;
 fUTCdata=0;
}

NcTimestamp::~NcTimestamp()
{
 
 if (fUTCdata)
 {
  delete fUTCdata;
  fUTCdata=0;
 }
}

NcTimestamp::NcTimestamp(const NcTimestamp& t) : TTimeStamp(t)
{
 
 fMJD=t.fMJD;
 fJsec=t.fJsec;
 fJns=t.fJns;
 fJps=t.fJps;
 fCalcs=t.fCalcs;
 fCalcns=t.fCalcns;
 fUtc=t.fUtc;
 fLeap=t.fLeap;
 fDut=t.fDut;
 fTmjd=t.fTmjd;
 fTsec=t.fTsec;
 fTns=t.fTns;
 fTps=t.fTps;
 fUTCdata=0;
 TTree* tx=t.fUTCdata;
 if (tx) fUTCdata=(TTree*)tx->Clone();
}

NcTimestamp& NcTimestamp::operator=(const NcTimestamp& t)
{
 
 NcTimestamp ts(t);
 
 fUTCdata=0;
 TTree* tx=ts.fUTCdata;
 if (tx) fUTCdata=(TTree*)tx->Clone();
 
 Int_t mjd=0;
 Int_t sec=0;
 Int_t ns=0;
 Int_t ps=0;
 Int_t leap=0;
 Double_t dut=0;
 
 ts.GetMJD(mjd,sec,ns);
 ps=ts.fJps;
 
 Int_t flag=ts.GetUTCparameters(leap,dut);
 
 TString utc="N";
 if (flag==1) utc="M";
 if (flag==-1) utc="A";
 if (flag==-2) utc="U";
 if (flag==-3) utc="U";
 
 SetMJD(mjd,sec,ns,ps,utc,leap,dut);
 
 // Correct for double setting of dut
 if (utc=="A" || utc=="M") AddSec(-dut);
 
 return *this;
}

void NcTimestamp::Date(Int_t mode,Double_t offset)
{
 
 Int_t mjd,mjsec,mjns;
 GetMJD(mjd,mjsec,mjns);
 
 TString month[12]={"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};
 TString day[7]={"Mon","Tue","Wed","Thu","Fri","Sat","Sun"};
 UInt_t y,m,d,wd;
 Int_t hh,mm,ss,ns,ps;
 Double_t gat,gast;
 Bool_t date=kFALSE;
 
 // The UT date and time
 if (abs(mode)==1 || abs(mode)==3)
 {
  if (mjd>=40587 && (mjd<65442 || (mjd==65442 && mjsec<8047)))
  {
   GetDate(kTRUE,0,&y,&m,&d);
   wd=GetDayOfWeek(kTRUE,0);
   printf(" %-s, %02d %-s %-d ",day[wd-1].Data(),d,month[m-1].Data(),y);
   date=kTRUE;
  }
  else
  {
   cout << " Time ";
   date=kFALSE;
  }
  GetUT(hh,mm,ss,ns,ps);
  printf("%02d:%02d:%02d.%09d%03d",hh,mm,ss,ns,ps);
  if (!fUtc)
  {
   cout << " (UTC) ";
  }
  else
  {
   cout << " (UT1) ";
  }
 
  // The GMST time information
  GetGMST(hh,mm,ss,ns,ps);
  printf("%02d:%02d:%02d.%09d%03d (GMST)\n",hh,mm,ss,ns,ps);
 
  // Equation of Time and Equation of Equinoxes
  Double_t eot;
  Double_t eox=Almanac(0,0,0,0,"Sun",0,0,0,&eot,10);
  // Convert to fractional hours
  eox/=3600.;
  eot/=3600.;
  if (date)
  {
   cout << "                 ";
  }
  else
  {
   cout << "     ";
  }
  if (eot>=0) cout << " ";
  PrintTime(eot,3); cout << " (LAT-LMT)";
  cout << "     ";
  if (eox>=0) cout << " ";
  PrintTime(eox,3); cout << " (LAST-LMST)"; cout << endl;
 
  // Greenwich apparent time information
  if (mode<0)
  {
   gat=GetLT(eot); // Obtain GAT via the Equation of Time as offset
   gast=GetGAST();
   if (date)
   {
    cout << "                  ";
   }
   else
   {
    cout << "      ";
   }
   PrintTime(gat,12); cout << " (GAT) ";
   PrintTime(gast,12); cout << " (GAST)" << endl;
  }
 
  // Local time information
  if (offset)
  {
   // Determine the new date by including the offset
   NcTimestamp t2(*this);
   t2.Add(offset);
   Int_t mjd2,mjsec2,mjns2;
   t2.GetMJD(mjd2,mjsec2,mjns2);
   if (mjd2>=40587 && (mjd2<65442 || (mjd2==65442 && mjsec2<8047)))
   {
    t2.GetDate(kTRUE,0,&y,&m,&d);
    wd=t2.GetDayOfWeek(kTRUE,0);
    printf(" %-s, %02d %-s %-d ",day[wd-1].Data(),d,month[m-1].Data(),y);
    date=kTRUE;
   }
   else
   {
    cout << " Time ";
    date=kFALSE;
   }
   // Determine the local time by including the offset w.r.t. the original timestamp
   Double_t hlt=0;
   Double_t hlst=0;
   if (mode>0)
   {
    hlt=GetLT(offset);
    hlst=GetLMST(offset);
    PrintTime(hlt,12);  cout << " (LMT) ";
    PrintTime(hlst,12); cout << " (LMST)" << endl;
   }
   else
   {
    hlt=GetLT(offset+eot); // Obtain LAT via the Equation of Time as extra offset
    hlst=GetLAST(offset);
    PrintTime(hlt,12);  cout << " (LAT) ";
    PrintTime(hlst,12); cout << " (LAST)" << endl;
   }
  }
 }
 
 // Julian parameter information
 if (abs(mode)==2 || abs(mode)==3)
 {
  Int_t jd,jsec,jns;
  GetJD(jd,jsec,jns);
  Int_t tjd,tjsec,tjns;
  GetTJD(tjd,tjsec,tjns);
  printf(" Julian Epoch : %-.20f Besselian Epoch : %-.20f\n",GetJE(),GetBE());
  printf(" JD  : %7d sec : %5d ns : %9d ps : %3d Fractional : %25.17f\n",jd,jsec,jns,fJps,GetJD());
  printf(" MJD : %7d sec : %5d ns : %9d ps : %3d Fractional : %25.17f\n",mjd,mjsec,mjns,fJps,GetMJD());
  printf(" TJD : %7d sec : %5d ns : %9d ps : %3d Fractional : %25.17f\n",tjd,tjsec,tjns,fJps,GetTJD());
  if (fUtc && fUtc!=-3)
  {
   printf(" TAI : %7d sec : %5d ns : %9d ps : %3d Fractional : %25.17f\n",fTmjd,fTsec,fTns,fTps,GetTAI());
  }
 }
 
 // TAI related information
 if (mode==4 && fUtc && fUtc!=-3)
 {
  printf( " Cumulated (TAI-UTC) leap seconds: %-3d UT1-UTC : %-.6f sec.",fLeap,fDut);
  if (fUtc<0) cout << " (IERS database)" << endl;
  if (fUtc>0) cout << " (Manual setting)" << endl;
 
  // A dummy timestamp is used to obtain the TAI corresponding date indicator
  NcTimestamp tx;
  tx.SetMJD(fTmjd,fTsec,fTns,fTps);
 
  if (fTmjd>=40587 && (fTmjd<65442 || (fTmjd==65442 && fTsec<8047)))
  {
   tx.GetDate(kTRUE,0,&y,&m,&d);
   wd=tx.GetDayOfWeek(kTRUE,0);
   printf(" %-s, %02d %-s %-d ",day[wd-1].Data(),d,month[m-1].Data(),y);
   date=kTRUE;
  }
  else
  {
   cout << " Time ";
   date=kFALSE;
  }
 
  // Determine the TAI derived times
  GetTAI(hh,mm,ss,ns,ps,"TAI");
  printf("%02d:%02d:%02d.%09d%03d (TAI) ",hh,mm,ss,ns,ps);
 
  GetTAI(hh,mm,ss,ns,ps,"UTC");
  printf("%02d:%02d:%02d.%09d%03d (UTC)\n",hh,mm,ss,ns,ps);
 
  GetTAI(hh,mm,ss,ns,ps,"GPS");
  if (!date)
  {
    cout << "      ";
  }
  else
  {
   cout << "                  ";
  }
  printf("%02d:%02d:%02d.%09d%03d (GPS) ",hh,mm,ss,ns,ps);
 
  GetTAI(hh,mm,ss,ns,ps,"TT");
  printf("%02d:%02d:%02d.%09d%03d (TT)\n",hh,mm,ss,ns,ps);
 }
}

Double_t NcTimestamp::GetJD(Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Int_t ss,Int_t ns) const
{
 
 if (y<0 || m<1 || m>12 || d<1 || d>31) return -1;
 if (hh<0 || hh>23 || mm<0 || mm>59 || ss<0 || ss>59 || ns<0 || ns>1e9) return -1;
 
 // The UT daytime in fractional hours
 Double_t ut=double(hh)+double(mm)/60.+(double(ss)+double(ns)*1.e-9)/3600.;
 
 Double_t JD=0;
 
 JD=367*y-int(7*(y+int((m+9)/12))/4)
    -int(3*(int((y+(m-9)/7)/100)+1)/4)
    +int(275*m/9)+d+1721028.5+ut/24.;
 
 return JD;
}

Double_t NcTimestamp::GetMJD(Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Int_t ss,Int_t ns) const
{
 
 Double_t JD=GetJD(y,m,d,hh,mm,ss,ns);
 
 if (JD<0) return JD;
 
 Double_t MJD=JD-2400000.5;
 
 return MJD; 
} 

Double_t NcTimestamp::GetTJD(Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Int_t ss,Int_t ns) const
{
 
 Double_t JD=GetJD(y,m,d,hh,mm,ss,ns);
 
 if (JD<0) return JD;
 
 Double_t TJD=JD-2440000.5;
 
 return TJD; 
} 

Double_t NcTimestamp::GetJE(Double_t date,TString mode) const
{
 
 if ((mode != "jd") && (mode != "mjd") && (mode != "tjd")) return -99999;
 
 Double_t jd=date;
 if (mode=="mjd") jd=date+2400000.5;
 if (mode=="tjd") jd=date+2440000.5;
 
 Double_t je=2000.+(jd-2451545.)/365.25;
 
 return je;
}

Double_t NcTimestamp::GetBE(Double_t date,TString mode) const
{
 
 if ((mode != "jd") && (mode != "mjd") && (mode != "tjd")) return -99999;
 
 Double_t jd=date;
 if (mode=="mjd") jd=date+2400000.5;
 if (mode=="tjd") jd=date+2440000.5;
 
 Double_t be=1900.+(jd-2415020.31352)/365.242198781;
 
 return be;
}

void NcTimestamp::Convert(Double_t date,Int_t& days,Int_t& secs,Int_t& ns) const
{
 
 days=int(date);
 date=date-double(days);
 Int_t daysecs=24*3600;
 date=date*double(daysecs);
 secs=int(date);
 date=date-double(secs);
 ns=int(date*1.e9);
}

Double_t NcTimestamp::Convert(Int_t days,Int_t secs,Int_t ns) const
{
 
 Double_t frac=double(secs)+double(ns)*1.e-9;
 Int_t daysecs=24*3600;
 frac=frac/double(daysecs);
 Double_t date=double(days)+frac;
 return date;
}

void NcTimestamp::Convert(Double_t h,Int_t& hh,Int_t& mm,Int_t& ss,Int_t& ns,Int_t& ps) const
{
 
 // Neglect sign of h
 h=fabs(h);
 
 hh=int(h);
 h=h-double(hh);
 h=h*60.;
 mm=int(h);
 h=h-double(mm);
 h=h*60.;
 ss=int(h);
 h=h-double(ss);
 h=h*1.e9;
 ns=int(h);
 h=h-double(ns);
 h=h*1000.;
 ps=int(h);
}

void NcTimestamp::Convert(Double_t h,Int_t& hh,Int_t& mm,Double_t& ss) const
{
 
 // Neglect sign of h
 h=fabs(h);
 
 hh=int(h);
 h=h-double(hh);
 h=h*60.;
 mm=int(h);
 h=h-double(mm);
 ss=h*60.;
}

Double_t NcTimestamp::Convert(Int_t hh,Int_t mm,Int_t ss,Int_t ns,Int_t ps) const
{
 
 // Neglect the sign of the input values
 hh=abs(hh);
 mm=abs(mm);
 ss=abs(ss);
 ns=abs(ns);
 ps=abs(ps);
 
 Double_t h=hh;
 h+=double(mm)/60.+(double(ss)+double(ns)*1.e-9+double(ps)*1.e-12)/3600.;
 
 return h;
}

Double_t NcTimestamp::Convert(Int_t hh,Int_t mm,Double_t ss) const
{
 
 // Neglect the sign of the input values
 hh=abs(hh);
 mm=abs(mm);
 ss=fabs(ss);
 
 Double_t h=hh;
 h+=double(mm)/60.+ss/3600.;
 
 return h;
}

void NcTimestamp::PrintTime(Double_t h,Int_t ndig) const
{
 
 Int_t hh,mm,ss;
 ULong64_t sfrac;
 Double_t s;
 
 while (h<-24)
 {
  h+=24.;
 }
 while (h>24)
 {
  h-=24.;
 }
   
 Convert(h,hh,mm,s);
 ss=Int_t(s);
 s-=Double_t(ss);
 s*=pow(10.,ndig);
 sfrac=ULong64_t(s);
 
 if (h<0) printf("%-s","-");
 printf("%02d:%02d:%02d.%0*llu",hh,mm,ss,ndig,sfrac);
}

void NcTimestamp::FillJulian()
{
 
 UInt_t y,m,d,hh,mm,ss;
 
 GetDate(kTRUE,0,&y,&m,&d);
 GetTime(kTRUE,0,&hh,&mm,&ss);
 Int_t ns=GetNanoSec();
 
 Double_t mjd=GetMJD(y,m,d,hh,mm,ss,ns);
 
 fMJD=int(mjd);
 fJsec=GetSec()%(24*3600); // Daytime in elapsed seconds
 fJns=ns;                  // Remaining fractional elapsed second in nanoseconds
 
 // Store the TTimeStamp seconds and nanoseconds values
 // for which this Julian calculation was performed.
 fCalcs=GetSec();
 fCalcns=GetNanoSec();
}

void NcTimestamp::GetMJD(Int_t& mjd,Int_t& sec,Int_t& ns)
{
 
 if (fCalcs != GetSec() || fCalcns != GetNanoSec())
 {
  FillJulian();
  SetUTCparameters("A",0,0);
 }
 
 mjd=fMJD;
 sec=fJsec;
 ns=fJns;
}

Double_t NcTimestamp::GetMJD()
{
 
 Int_t mjd=0;
 Int_t sec=0;
 Int_t ns=0;
 GetMJD(mjd,sec,ns);
 
 Double_t date=Convert(mjd,sec,ns);
 
 return date;
}

void NcTimestamp::GetTJD(Int_t& tjd,Int_t& sec,Int_t& ns)
{
 
 Int_t mjd=0;
 GetMJD(mjd,sec,ns);
 
 tjd=mjd-40000;
}

Double_t NcTimestamp::GetTJD()
{
 
 Int_t tjd=0;
 Int_t sec=0;
 Int_t ns=0;
 GetTJD(tjd,sec,ns);
 
 Double_t date=Convert(tjd,sec,ns);
 
 return date;
}

Int_t NcTimestamp::GetTAI(Int_t& d,Int_t& sec,Int_t& ns,Int_t& ps,Bool_t tmjd)
{
 
 // Make sure to have the updated parameters
 GetMJD(d,sec,ns);
 FillTAI();
 
 d=0;
 sec=0;
 ns=0;
 ps=0;
 
 if (!fUtc) return 0;
 
 d=fTmjd;
 sec=fTsec;
 ns=fTns;
 ps=fTps;
 
 if (!tmjd) d-=36204;
 
 return fUtc;
}

Double_t NcTimestamp::GetTAI(Bool_t tmjd)
{
 
 if (!fUtc) return 0;
 
 Int_t d=0;
 Int_t s=0;
 Int_t ns=0;
 Int_t ps=0;
 GetTAI(d,s,ns,ps,tmjd);
 
 Double_t days=Convert(d,s,ns);
 
 return days;
}

Int_t NcTimestamp::GetTAI(Int_t& hh,Int_t& mm,Int_t& ss,Int_t& ns,Int_t& ps,TString type)
{
 
 hh=0;
 mm=0;
 ss=0;
 ns=0;
 ps=0;
 
 if (type!="TAI" && type!="UTC" && type!="GPS" && type!="TT") return 0;
 
 Int_t d,sec,nsec,psec;
 
 // Use a dummy timestamp to easily correct for the various offsets
 NcTimestamp tx=(*this);
 if (type=="UTC") tx.Add(0,-fLeap,0,0);
 if (type=="GPS") tx.Add(0,-19,0,0);
 if (type=="TT") tx.Add(0,32,184000000,0);
 
 // Use the UTC parameters of the actual (unmodified) timestamp
 tx.fLeap=fLeap;
 tx.fDut=fDut;
 
 tx.GetTAI(d,sec,nsec,psec);
 
 hh=sec/3600;
 sec=sec%3600;
 mm=sec/60;
 ss=sec%60;
 ns=nsec;
 ps=psec;
 
 return fUtc;
}

Double_t NcTimestamp::GetUnixTime()
{
 
 Double_t t=0;
 
 Int_t days=0;
 Int_t secs=0;
 Int_t ns=0;
 Int_t mjd=0;
 GetMJD(mjd,secs,ns);
 Int_t ps=GetPs();
 
 // Get UTC from the stored UT1 via UTC=UT1-dUT
 if (fUtc)
 {
  NcTimestamp tx;
  tx.SetMJD(mjd,secs,ns,ps);
  tx.AddSec(-fDut);
  tx.GetMJD(mjd,secs,ns);
  ps=tx.GetPs();
 }
 
 days=mjd-40587;
 
 t=(days*86400)+secs;
 t+=double(ns)*1.e-9+double(ps)*1.e-12;
 
 return t;
}

void NcTimestamp::GetJD(Int_t& jd,Int_t& sec,Int_t& ns)
{
 
 Int_t mjd=0;
 GetMJD(mjd,sec,ns);
 
 jd=mjd+2400000;
 sec+=12*3600;
 if (sec >= 24*3600)
 {
  sec-=24*3600;
  jd+=1;
 }
}

Double_t NcTimestamp::GetJD()
{
 
 Int_t jd=0;
 Int_t sec=0;
 Int_t ns=0;
 GetJD(jd,sec,ns);
 
 Double_t date=Convert(jd,sec,ns);
 
 return date;
}

Double_t NcTimestamp::GetJE()
{
 
 Double_t jd=GetJD();
 Double_t je=GetJE(jd);
 return je;
}

Double_t NcTimestamp::GetBE()
{
 
 Double_t jd=GetJD();
 Double_t be=GetBE(jd);
 return be;
}

void NcTimestamp::SetMJD(Int_t mjd,Int_t sec,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut)
{
 
 if (sec<0 || sec>=24*3600 || ns<0 || ns>=1e9 || ps<0 || ps>=1000)
 {
  cout << " *NcTimestamp::SetMJD* Invalid input."
       << " sec : " << sec << " ns : " << ns << " ps : " << ps << endl; 
  return;
 }
 
 fMJD=mjd;
 fJsec=sec;
 fJns=ns;
 fJps=ps;
 
 Int_t epoch=40587; // MJD of the start of the epoch
 Int_t limit=65442; // MJD of the latest possible TTimeStamp date/time
 
 Int_t date,time;
 if (mjd<epoch || mjd>limit || (mjd==limit && sec>=8047))
 {
  Set(0,kFALSE,0,kFALSE);
  date=GetDate();
  time=GetTime();
  Set(date,time,0,kTRUE,0);
 }
 else
 {
  // The elapsed time since start of EPOCH
  Int_t days=mjd-epoch;
  UInt_t secs=days*24*3600;
  secs+=sec;
  Set(secs,kFALSE,0,kFALSE);
  date=GetDate();
  time=GetTime();
  Set(date,time,ns,kTRUE,0);
 }
 
 // Denote that the Julian and TTimeStamp parameters are synchronised,
 // even in the case the MJD falls outside the TTimeStamp validity range.
 // The latter still allows retrieval of Julian parameters for these
 // earlier times.
 fCalcs=GetSec();
 fCalcns=GetNanoSec();
 
 // Update the UTC parameters and corresonding TAI time recording
 SetUTCparameters(utc,leap,dut);
 
 if (utc=="A" || utc=="M") AddSec(fDut);
}

void NcTimestamp::SetMJD(Double_t mjd,TString utc,Int_t leap,Double_t dut)
{
 
 Int_t days=0;
 Int_t secs=0;
 Int_t ns=0;
 Convert(mjd,days,secs,ns);
 SetMJD(days,secs,ns,0,utc,leap,dut);
}

void NcTimestamp::SetJD(Int_t jd,Int_t sec,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut)
{
 
 Int_t mjd=jd-2400000;
 sec-=12*3600;
 if (sec<0)
 {
  sec+=24*3600;
  mjd-=1;
 }
 
 SetMJD(mjd,sec,ns,ps,utc,leap,dut);
}

void NcTimestamp::SetJD(Double_t jd,TString utc,Int_t leap,Double_t dut)
{
 
 Int_t days=0;
 Int_t secs=0;
 Int_t ns=0;
 Convert(jd,days,secs,ns);
 
 SetJD(days,secs,ns,0,utc,leap,dut);
}

void NcTimestamp::SetTJD(Int_t tjd,Int_t sec,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut)
{
 
 Int_t mjd=tjd+40000;
 
 SetMJD(mjd,sec,ns,ps,utc,leap,dut);
}

void NcTimestamp::SetTJD(Double_t tjd,TString utc,Int_t leap,Double_t dut)
{
 
 Int_t days=0;
 Int_t secs=0;
 Int_t ns=0;
 Convert(tjd,days,secs,ns);
 
 SetTJD(days,secs,ns,0,utc,leap,dut);
}

void NcTimestamp::FillTAI()
{
 
 if (!fUtc || fUtc==-3)
 {
  fTmjd=0;
  fTsec=0;
  fTns=0;
  fTps=0;
  return;
 }
 
 // Use memberfunction to ensure most recent values for the stored UT1 
 GetMJD(fTmjd,fTsec,fTns);
 fTps=GetPs();
 
 // Dummy timestamp to easily obtain TAI based day etc. counts
 // It is essential not to use UTC parameters here in order to prevent an infinite loop
 NcTimestamp tx;
 tx.SetMJD(fTmjd,fTsec,fTns,fTps,"N"); // This stores UT1 as if it were UTC
 
 tx.AddSecCalc(-fDut,kFALSE);    // Convert the stored UT1 into UTC via UTC=UT1-dUT
 tx.AddCalc(0,fLeap,0,0,kFALSE); // Account for the leap seconds
 
 // Retrieve the corresponding TAI day etc. count
 tx.GetMJD(fTmjd,fTsec,fTns);
 fTps=tx.GetPs();
}

Int_t NcTimestamp::SetTAI(TString type,TString date,TString time,Int_t mode,TString utc,Int_t leap,Double_t dut)
{
 
 Int_t ibad=0;
 
 if (type!="UTC" && type!="GPS" && type!="TAI" && type!="TT") ibad=1;
 
 if (utc!="M" && utc!="A") ibad=1;
 
 if (utc=="M" && fabs(dut)>0.9) ibad=1;
 
 if (utc=="A" && !fUTCdata) ibad=1;
 
 // In case utc="A" check whether the corresponding IERS database info is available
 if (utc=="A")
 {
  NcTimestamp tx; // Dummy timestamp for easy MJD retrieval
  tx.SetUT(date,time,mode,utc);
  Int_t ien=GetUTCparameters(tx.fMJD,leap,dut);
  if (ien<0) ibad=1;
 }
 
 if (ibad)
 {
  SetJD(0,"N");
  return fUtc;
 }
 
 SetUT(date,time,mode,utc,leap,dut);
 if (type != "UTC") Add(0,-fLeap,0,0);      // Account for the leap seconds
 if (type == "GPS") Add(0,19,0,0);          // Account for TAI-GPS=19 sec.
 if (type == "TT") Add(0,-32,-184000000,0); // Account for TAI-TT=-32.184 sec.
 
 return fUtc;
}

Int_t NcTimestamp::SetTAI(Int_t d,Int_t sec,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut,Bool_t tmjd)
{
 
 Int_t ibad=0;
 
 if (sec<0 || sec>86400 || ns<0 || ns>999999999 || ps<0 || ps>999) ibad=1;
 
 if (utc!="M" && utc!="A") ibad=1;
 
 if (utc=="M" && fabs(dut)>0.9) ibad=1;
 
 if (utc=="A" && !fUTCdata) ibad=1;
 
 // Set the corresponding MJD
 Int_t mjd=d;
 if (!tmjd) mjd+=36204;
 
 // In case utc="A" check whether the corresponding IERS database info is available
 if (utc=="A")
 {
  Int_t ien=GetUTCparameters(mjd,leap,dut);
  if (ien<0) ibad=1;
 }
 
 if (ibad)
 {
  SetJD(0,"N");
  return fUtc;
 }
 
 SetMJD(mjd,sec,ns,ps);
 SetUTCparameters(utc,leap,dut); // Update the UTC parameters and corresonding TAI time recording
 Add(0,-fLeap,0,0); // Account for the leap seconds
 AddSec(fDut);      // Account for dUT=UT1-UTC
 
 // Set the corresponding TAI day count etc.
 FillTAI();
 
 return fUtc;
}

Int_t NcTimestamp::SetTAI(Double_t tai,TString utc,Int_t leap,Double_t dut,Bool_t tmjd)
{
 
 Int_t days=0;
 Int_t secs=0;
 Int_t ns=0;
 Convert(tai,days,secs,ns);
 
 SetTAI(days,secs,ns,0,utc,leap,dut,tmjd);
 
 return fUtc;
}

Int_t NcTimestamp::SetGPS(Int_t w,Int_t sow,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut,Int_t icycle)
{
 
 if (w<0 || sow<0 || sow>604800 || ns<0 || ns>999999999 || ps<0 || ps>999 || icycle<0)
 {
  SetJD(0,"N");
  return fUtc;
 }
 
 // Correct the week count for the cycle number if needed
 if (icycle) w+=icycle*1024;
 
 Int_t days=8040+w*7;
 sow+=19;
 Int_t daysecs=24*3600;
 Int_t days2=sow/daysecs;
 days+=days2;
 Int_t secs=sow%daysecs;
 
 SetTAI(days,secs,ns,ps,utc,leap,dut);
 
 return fUtc;
}

Int_t NcTimestamp::SetGPS(Int_t w,Int_t dow,Int_t sod,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut,Int_t icycle)
{
 
 if (w<0 || dow<0 || dow>7 || sod<0 || sod>86400 || ns<0 || ns>999999999 || ps<0 || ps>999 || icycle<0)
 {
  SetJD(0,"N");
  return fUtc;
 }
 
 // Correct the week count for the cycle number if needed
 if (icycle) w+=icycle*1024;
 
 Int_t days=8040+w*7+dow;
 sod+=19;
 
 SetTAI(days,sod,ns,ps,utc,leap,dut);
 
 return fUtc;
}

Int_t NcTimestamp::SetUnixTime(Double_t sec,TString utc,Int_t leap,Double_t dut)
{
 
 // Determine the fractional day count since the start of the Unix Epoch
 Double_t tday=sec/86400.;
 
 Int_t days=0;
 Int_t s=0;
 Int_t ns=0;
 Convert(tday,days,s,ns);
 
 // Determine the remaining elapsed picoseconds
 Int_t iword=int(sec);
 sec=sec-double(iword);
 sec*=1e9;
 iword=int(sec);
 sec=sec-double(iword);
 Int_t ps=int(sec*1000.);
 
 // Unix time is related to UTC and not to UT1.
 if (utc=="U") utc="A";
 
 SetMJD(40587,0,0,0); // Start of the Unix Epoch
 SetUTCparameters(utc,leap,dut); // Update the UTC parameters and corresonding TAI time recording
 Add(days,s,ns,ps);      // Add the elapsed time
 if (fUtc) AddSec(fDut); // Correct for dUT=UT1-UTC
 
 return fUtc;
}

void NcTimestamp::SetNs(Int_t ns)
{
 
 if (ns>=0 && ns<=999999999) fJns=ns; 
}

Int_t NcTimestamp::GetNs() const
{
 
 return fJns; 
}

void NcTimestamp::SetPs(Int_t ps)
{
 
 if (ps>=0 && ps<=999) fJps=ps; 
}

Int_t NcTimestamp::GetPs() const
{
 
 return fJps; 
}

Int_t NcTimestamp::GetUTCparameters(Int_t& leap,Double_t& dut) const
{
 
 leap=fLeap;
 dut=fDut;
 
 return fUtc; 
}

Int_t NcTimestamp::GetUTCparameters(Int_t mjd,Int_t& leap,Double_t& dut) const
{
 
 leap=0;
 dut=0;
 
 if (!fUTCdata) return -1;
 
 Int_t nen=fUTCdata->GetEntries();
 
 if (!nen) return -1;
 
 Int_t dbmjd=0;
 Int_t dbleap=0;
 Double_t dbdut=0;
 
 fUTCdata->SetBranchAddress("mjd",&dbmjd);
 fUTCdata->SetBranchAddress("lsec",&dbleap);
 fUTCdata->SetBranchAddress("dut",&dbdut);
 
 // Data of the first entry
 fUTCdata->GetEntry(0);
 Int_t ien=mjd-dbmjd;
 
 if (ien<0 || ien>=nen) return -1; // Specified mjd not in range of database
 
 fUTCdata->GetEntry(ien);
 if (dbmjd==mjd) // Specified mjd is found in database
 {
  leap=dbleap;
  dut=dbdut;
 }
 else
 {
  ien=-1;
  leap=0;
  dut=0;
 }
 
 return ien;
}

Int_t NcTimestamp::SetUTCparameters(TString utc,Int_t leap,Double_t dut)
{
 
 fUtc=0;
 if (utc=="U") fUtc=-3;
 fLeap=0;
 fDut=0;
 
 Int_t ibad=0;
 
 if (utc!="N" && utc!="M" && utc!="A" && utc!="U") ibad=1;
 
 if (utc=="N" || ((utc=="A" || utc=="U") && !fUTCdata)) ibad=1;
 
 if (utc=="M" && fabs(dut)>0.9) ibad=1;
 
 if (ibad)
 {
  FillTAI();
  return fUtc;
 }
 
 // From here only utc="M", utc="A" or utc="U"
 
 if (utc=="M")
 {
  fUtc=1;
  fLeap=leap;
  fDut=dut;
 
  FillTAI();
  return fUtc;
 }
 
 // Automatic setting of the UTC parameters from the loaded data files
 Int_t nen=fUTCdata->GetEntries();
 
 if (nen<=0) // No entries in the IERS data TTree
 {
  FillTAI();
  return fUtc;
 }
 
 Int_t mjd=0;
 
 fUTCdata->SetBranchAddress("mjd",&mjd);
 fUTCdata->SetBranchAddress("lsec",&leap);
 fUTCdata->SetBranchAddress("dut",&dut);
 
 // Get the starting mjd of the IERS daily data
 // and determine the entry for the current MJD info
 fUTCdata->GetEntry(0);
 Int_t ien=fMJD-mjd;
 if (ien>=0 && ien<nen)
 {
  fUTCdata->GetEntry(ien);
  if (mjd==fMJD)
  {
   fUtc=-1;
   if (utc=="U") fUtc=-2;
   fLeap=leap;
   fDut=dut;
  }
 }
 
 FillTAI();
 
 return fUtc;
}

TTree* NcTimestamp::LoadUTCparameterFiles(TString leapfile,TString dutfile)
{
 
 // Expand the input file pathnames 
 leapfile=gSystem->ExpandPathName(leapfile.Data());
 dutfile=gSystem->ExpandPathName(dutfile.Data());
  
 if (fUTCdata)
 {
  delete fUTCdata;
  fUTCdata=0;
 }
 
 // The Leap Second input data file 
 ifstream fleap;
 fleap.clear();
 fleap.open(leapfile.Data());
 if (!fleap.good())
 {
  cout << " *NcTimestamp::LoadUTCparameterFiles* Data file for Leap Seconds not found ***" << endl;
  cout << " File name provided was : " << leapfile.Data() << endl;
  return 0;
 }
 
 // The dUT=UT1-UTC input data file 
 ifstream fdut;
 fdut.clear();
 fdut.open(dutfile.Data());
 if (!fdut.good())
 {
  cout << " *NcTimestamp::LoadUTCparameterFiles* Data file for dUT=UT-UTC not found ***" << endl;
  cout << " File name provided was : " << dutfile.Data() << endl;
  return 0;
 }
 
 // Determine the number of characters in the Leap Second input file
 // to reserve sufficient array storage of all Leap Second entries
 fleap.seekg(0,fleap.end); // Position at end of file
 Int_t ndim=fleap.tellg();
 
 // The storage arrays for the Leap second data
 Int_t* lmjd=new Int_t[ndim];
 Int_t* leap=new Int_t[ndim];

 // Read the Leap Second data //
 
 fleap.seekg(0); // Position at begin of file
 
 // Read title lines until the first data line is found
 string line;
 Int_t i=0;
 while (getline(fleap,line))
 {
  if (line.find("1972")!=line.npos) break;
  i++;
 }
 
 // Go to the beginning of the file and skip the title lines preceding the data lines
 fleap.seekg(0);
 for (Int_t j=0; j<i; j++)
 {
  getline(fleap,line);
 }
 
 // Read the data
 Float_t rmjd=0;
 Int_t lsec=0;
 Float_t x; // Dummy variable for skipping non-requested data columns
 i=0;
 while (fleap >> rmjd >> x >> x >> x >> lsec)
 {
  lmjd[i]=int(rmjd);
  leap[i]=lsec;
  i++;
 }
 
 // The number of actual Leap Second entries
 Int_t nleap=i;

 // Read the dUT=UT1-UTC //
 
 fdut.seekg(0); // Position at begin of file
 
 // Read title lines until the first data line is found
 i=0;
 while (getline(fdut,line))
 {
  if (line.find("1962")!=line.npos) break;
  i++;
 }
 
 // Go to the beginning of the file and skip the title lines preceding the data lines
 fdut.seekg(0);
 for (Int_t j=0; j<i; j++)
 {
  getline(fdut,line);
 }
 
 // Read the dUT daily data and fill the TTree structure
 Int_t mjd=0;
 Double_t dut=0;
 
 // The produced output structure
 fUTCdata=new TTree("UTCdata","Daily UTC leap second and dUT=UT-UTC parameter data");
 
 // Prevent this Tree to be automatically coupled to a user defined output file
 fUTCdata->SetDirectory(0);
 
 // The output variables for the Tree
 fUTCdata->Branch("mjd",&mjd,"mjd/I");
 fUTCdata->Branch("lsec",&lsec,"lsec/I");
 fUTCdata->Branch("dut",&dut,"dut/D");
 
 while (fdut >> x >> x >> x >> mjd >> x >> x >> dut >> x >> x >> x >> x >> x >> x >> x >> x >> x)
 {
  lsec=0;
  // Retrieve the corresponding Leap Second info
  for (Int_t j=nleap-1; j>=0; j--)
  {
   if (mjd>=lmjd[j])
   {
    lsec=leap[j];
    break;
   }
  }
  fUTCdata->Fill();
 }
 
 delete[] lmjd;
 delete[] leap; 
 
 // Retrieve the UTC parameters for the current timestamp
 if (!fUtc) // Reference time is UTC
 {
  SetUTCparameters("A",0,0);
  AddSec(fDut);
 }
 else // Reference time is UT1
 {
  SetUTCparameters("U",0,0);
 }
 
 return fUTCdata;
}

TTree* NcTimestamp::GetIERSdatabase() const
{
 
 return fUTCdata;
}

void NcTimestamp::Add(Int_t d,Int_t s,Int_t ns,Int_t ps)
{
 
 AddCalc(d,s,ns,ps);
}

void NcTimestamp::Add(Double_t hours)
{
 
 AddCalc(hours);
}

void NcTimestamp::AddSec(Double_t seconds)
{
 
 AddSecCalc(seconds);
}

void NcTimestamp::AddCalc(Int_t d,Int_t s,Int_t ns,Int_t ps,Bool_t utcpar)
{
 
 Int_t days=0;
 Int_t secs=0;
 Int_t nsec=0;
 // Use Get functions to ensure updated Julian parameters. 
 GetMJD(days,secs,nsec);
 Int_t psec=GetPs();
 
 psec+=ps%1000;
 nsec+=ps/1000;
 while (psec<0)
 {
  nsec-=1;
  psec+=1000;
 }
 while (psec>999)
 {
  nsec+=1;
  psec-=1000;
 }
 
 nsec+=ns%1000000000;
 secs+=ns/1000000000;
 while (nsec<0)
 {
  secs-=1;
  nsec+=1000000000;
 }
 while (nsec>999999999)
 {
  secs+=1;
  nsec-=1000000000;
 }
 
 secs+=s%(24*3600);
 days+=s/(24*3600);
 while (secs<0)
 {
  days-=1;
  secs+=24*3600;
 }
 while (secs>=24*3600)
 {
  days+=1;
  secs-=24*3600;
 }
 
 days+=d;
 
 TString utc="N";
 if (fUtc==1) utc="M";
 if (fUtc==-1) utc="A";
 if (fUtc==-2) utc="U";
 if (fUtc==-3) utc="U";
 Int_t leap=fLeap;
 Double_t dut=fDut;
 SetMJD(days,secs,nsec,psec);
 
 if (utcpar) SetUTCparameters(utc,leap,dut);
}

void NcTimestamp::AddCalc(Double_t hours,Bool_t utcpar)
{
 
 Int_t d,s,ns,ps;
 Double_t h=fabs(hours);
 d=int(h/24.);
 h-=double(d)*24.;
 h*=3600.;
 s=int(h);
 h-=double(s);
 h*=1.e9;
 ns=int(h);
 h-=double(ns);
 ps=int(h*1000.);
 if (hours>0) AddCalc(d,s,ns,ps,utcpar);
 if (hours<0) AddCalc(-d,-s,-ns,-ps,utcpar);
}

void NcTimestamp::AddSecCalc(Double_t seconds,Bool_t utcpar)
{
 
 Int_t s,ns,ps;
 Double_t a=fabs(seconds);
 s=int(a);
 a-=double(s);
 a*=1.e9;
 ns=int(a);
 a-=double(ns);
 ps=int(a*1000.);
 if (seconds>0) AddCalc(0,s,ns,ps,utcpar);
 if (seconds<0) AddCalc(0,-s,-ns,-ps,utcpar);
}

Int_t NcTimestamp::GetDifference(NcTimestamp* t,Int_t& d,Int_t& s,Int_t& ns,Int_t& ps,TString type)
{
 
 d=0;
 s=0;
 ns=0;
 ps=0;
 
 if (!t || (type!="UT" && type!="TAI")) return 0;
 
 Int_t tUtc=t->fUtc;
 if (type=="TAI" && (!fUtc || !tUtc || fUtc==-3 || tUtc==-3)) return 0;
 
 Int_t d1=0;
 Int_t s1=0;
 Int_t ns1=0;
 Int_t ps1=0;
 
 Int_t d2=0;
 Int_t s2=0;
 Int_t ns2=0;
 Int_t ps2=0;
 
 NcTimestamp ts1(*t);
 NcTimestamp ts2(*this);
 
 // Check for a mixed use of UT1 and UTC
 if (type=="UT")
 {
  Double_t dut=0;
  if (!fUtc && tUtc) // This timestamp is in UTC and the other in UT1
  {
   dut=t->fDut;
   ts1.AddSec(-dut);
  }
  if (fUtc && !tUtc) // This timestamp is in UT1 and the other in UTC
  {
   ts2.AddSec(-fDut);
  }
 }
 
 // Use Get functions to ensure updated Julian and TAI parameters. 
 if (type=="UT")
 {
  ts1.GetMJD(d1,s1,ns1);
  ps1=ts1.GetPs();
  ts2.GetMJD(d2,s2,ns2);
  ps2=ts2.GetPs();
 }
 if (type=="TAI")
 {
  t->GetTAI(d1,s1,ns1,ps1);
  GetTAI(d2,s2,ns2,ps2);
 }
 
 d=d1-d2;
 s=s1-s2;
 ns=ns1-ns2;
 ps=ps1-ps2;
 
 if (!d && !s && !ns && !ps) return 0;
 
 Int_t sign=0;
 
 if (d>0) sign=1;
 if (d<0) sign=-1;
 
 if (!sign && s>0) sign=1;
 if (!sign && s<0) sign=-1;
 
 if (!sign && ns>0) sign=1; 
 if (!sign && ns<0) sign=-1;
 
 if (!sign && ps>0) sign=1; 
 if (!sign && ps<0) sign=-1;
 
 // In case the input stamp was earlier, take the reverse difference
 // to simplify the algebra.
 if (sign<0)
 {
  d=-d;
  s=-s;
  ns=-ns;
  ps=-ps;
 }
 
 // Here we always have a positive time difference
 // and can now unambiguously correct for other negative values.
 if (ps<0)
 {
  ns-=1;
  ps+=1000;
 }
 
 if (ns<0)
 {
  s-=1;
  ns+=1000000000;
 }
 
 if (s<0)
 {
  d-=1;
  s+=24*3600;
 }
 
 return sign;
}

Int_t NcTimestamp::GetDifference(NcTimestamp& t,Int_t& d,Int_t& s,Int_t& ns,Int_t& ps,TString type)
{
 
 return GetDifference(&t,d,s,ns,ps,type);
}

Double_t NcTimestamp::GetDifference(NcTimestamp* t,TString u,Int_t mode,TString type)
{
 
 if (!t || mode<1 || mode>3 || (type!="UT" && type!="TAI")) return 0;
 
 if (u!="d" && u!="s" && u!="ns" && u!="ps") return 0;
 
 
 Int_t tUtc=t->fUtc;
 if (type=="TAI" && (!fUtc || !tUtc || fUtc==-3 || tUtc==-3)) return 0;
 
 Double_t dt=0;
 
 Int_t d1=0;
 Int_t s1=0;
 Int_t ns1=0;
 Int_t ps1=0;
 
 Int_t d2=0;
 Int_t s2=0;
 Int_t ns2=0;
 Int_t ps2=0;
 
 NcTimestamp ts1(*t);
 NcTimestamp ts2(*this);
 
 // Check for a mixed use of UT1 and UTC
 if (type=="UT")
 {
  Double_t dut=0;
  if (!fUtc && tUtc) // This timestamp is in UTC and the other in UT1
  {
   dut=t->fDut;
   ts1.AddSec(-dut);
  }
  if (fUtc && !tUtc) // This timestamp is in UT1 and the other in UTC
  {
   ts2.AddSec(-fDut);
  }
 }
 
 // Use Get functions to ensure updated Julian and TAI parameters. 
 if (type=="UT")
 {
  ts1.GetMJD(d1,s1,ns1);
  ps1=ts1.GetPs();
  ts2.GetMJD(d2,s2,ns2);
  ps2=ts2.GetPs();
 }
 if (type=="TAI")
 {
  t->GetTAI(d1,s1,ns1,ps1);
  GetTAI(d2,s2,ns2,ps2);
 }
 
 Long_t dd=d1-d2;
 Long_t ds=s1-s2;
 Long_t dns=ns1-ns2;
 Long_t dps=ps1-ps2;
 
 // Time difference for the specified units only
 if (mode==3)
 {
  if (u=="d") dt=dd;
  if (u=="s") dt=ds;
  if (u=="ns") dt=dns;
  if (u=="ps") dt=dps;
  return dt;
 }
 
 // Suppress elapsed time for the larger units than specified
 if (mode==2)
 {
  if (u=="s") dd=0;
  if (u=="ns")
  {
   dd=0;
   ds=0;
  }
  if (u=="ps")
  {
   dd=0;
   ds=0;
   dns=0;
  }
 }
 
 // Compute the time difference as requested 
 if (u=="s" || u=="d")
 {
  // The time difference in (fractional) seconds
  dt=double(dd*24*3600+ds)+(double(dns)*1e-9)+(double(dps)*1e-12);
  if (u=="d") dt=dt/double(24*3600);
 }
 if (u=="ns") dt=(double(dd*24*3600+ds)*1e9)+double(dns)+(double(dps)*1e-3);
 if (u=="ps") dt=(double(dd*24*3600+ds)*1e12)+(double(dns)*1e3)+double(dps);
 
 return dt;
}

Double_t NcTimestamp::GetDifference(NcTimestamp& t,TString u,Int_t mode,TString type)
{
 
 return GetDifference(&t,u,mode,type);
}

void NcTimestamp::SetUT(Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Int_t ss,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut)
{
 
 if (d<1 || d>31 || m<1 || m>12 || hh<0 || hh>23 || mm<0 || mm>59 || ss<0 || ss>59 || ns<0 || ns>999999999 || ps<0 || ps>999)
 {
  cout << " *NcTimestamp::SetUT* Incompatible argument(s) Day=" << d << " Month=" << m << " Year=" << y
       << " hour=" << hh << " min=" << mm << " sec=" << ss << " ns=" << ns << " ps=" << ps << endl;
  cout << " ==> TAI related time recording is disabled and JD=0 has been set." << endl;
  SetJD(0,"N");
  return;
 }
 
 Int_t day=GetDayOfYear(d,m,y);
 Int_t secs=hh*3600+mm*60+ss;
 SetUT(y,day-1,secs,ns,ps,utc,leap,dut);
}

void NcTimestamp::SetUT(Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Double_t s,TString utc,Int_t leap,Double_t dut)
{
 
 Int_t ss=int(s);
 s-=double(ss);
 Int_t ns=s*1.e9;
 s-=double(ns)*1.e-9;
 Int_t ps=s*1.e12;
 SetUT(y,m,d,hh,mm,ss,ns,ps,utc,leap,dut);
}

void NcTimestamp::SetUT(Int_t y,Int_t m,Int_t d,TString time,TString utc,Int_t leap,Double_t dut)
{
 
 Long64_t iword=0;
 Int_t hh=0;
 Int_t mm=0;
 Int_t ss=0;
 Int_t ns=0;
 Int_t ps=0;
 time.ReplaceAll(":","");
 
 // Unpack the hhmmss integer part
 iword=time.Atoll();
 hh=iword/10000;
 iword=iword%10000;
 mm=iword/100;
 iword=iword%100;
 ss=iword;
 
 // Unpack the fractional part in ns and ps
 time.Remove(0,7); // Remove the integer part, including the decimal "."
 Int_t length=time.Length();
 time.Append('0',12-length); // Increase the string to represent the integer number of ns
 iword=time.Atoll();
 ns=iword/1000;
 ps=iword%1000;
 SetUT(y,m,d,hh,mm,ss,ns,ps,utc,leap,dut);
}

void NcTimestamp::SetUT(TString date,TString time,Int_t mode,TString utc,Int_t leap,Double_t dut)
{
 
 Int_t iword=0;
 Int_t utdate=0;
 Int_t year=0;
 Int_t month=0;
 Int_t day=0;
 TString datex=date;
 datex.ReplaceAll("-","");
 datex.ReplaceAll("/","");
 utdate=datex.Atoi();
 iword=utdate;
 if (mode==0)
 {
  day=iword/1000000;
  iword=iword%1000000;
  month=iword/10000;
  iword=iword%10000;
  year=iword;
 }
 if (mode==1)
 {
  year=iword/10000;
  iword=iword%10000;
  month=iword/100;
  iword=iword%100;
  day=iword;
 }
 if (mode==2)
 {
  month=iword/1000000;
  iword=iword%1000000;
  day=iword/10000;
  iword=iword%10000;
  year=iword;
 }
 if (mode==3)
 {
  year=iword/10000;
  iword=iword%10000;
  day=iword/100;
  iword=iword%100;
  month=iword;
 }
 
 if (day<1 || day>31 || month<1 || month>12)
 {
  cout << " *NcTimestamp::SetUT* Incompatible argument(s) Date: " << date << " Time: " << time << " mode: " << mode << endl;
  cout << " ==> TAI related time recording is disabled and JD=0 has been set." << endl;
  SetJD(0,"N");
 }
 else
 {
  SetUT(year,month,day,time,utc,leap,dut);
 }
}

void NcTimestamp::SetUT(Int_t y,Int_t d,Int_t s,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut)
{
 
 Double_t jd=GetJD(y,1,1,0,0,0,0);
 SetJD(jd);
 
 Int_t mjd,sec,nsec;
 GetMJD(mjd,sec,nsec);
 SetMJD(mjd,0,0,0);
 SetUTCparameters(utc,leap,dut); // Update the UTC parameters and corresonding TAI time recording
 Add(d,s,ns,ps);
 
 // Determine UT1 using dut=UT1-UTC if UTC was provided as input
 if (utc!="U") AddSec(fDut);
}

void NcTimestamp::GetUT(Int_t& hh,Int_t& mm,Int_t& ss,Int_t& ns,Int_t& ps)
{
 
 Int_t mjd,sec,nsec,psec;
 
 GetMJD(mjd,sec,nsec);
 psec=GetPs();
 
 hh=sec/3600;
 sec=sec%3600;
 mm=sec/60;
 ss=sec%60;
 ns=nsec;
 ps=psec;
}

Double_t NcTimestamp::GetUT()
{
 
 Int_t hh,mm,ss,ns,ps;
 
 GetUT(hh,mm,ss,ns,ps);
 
 Double_t ut=Convert(hh,mm,ss,ns,ps);
 
 return ut;
}

void NcTimestamp::GetGMST(Int_t& hh,Int_t& mm,Int_t& ss,Int_t& ns,Int_t& ps)
{
 
 Int_t mjd,sec,nsec,psec;
 
 // The current UT based timestamp data
 GetMJD(mjd,sec,nsec);
 psec=fJps;
 
 // The basis for the daily corrections in units of Julian centuries w.r.t. J2000.
 // Note : Epoch J2000 starts at 01-jan-2000 12:00:00 UT.
 Double_t tau=(GetJD()-2451545.)/36525.;
 
 // Syncronise sidereal time with current timestamp
 NcTimestamp sid;
 sid.SetMJD(mjd,sec,nsec,psec);
 
 // Add offset for GMST start value defined as 06:41:50.54841 at 01-jan 00:00:00 UT
 sec=6*3600+41*60+50;
 nsec=548410000;
 psec=0;
 sid.Add(0,sec,nsec,psec);
 
 // Daily correction for precession and polar motion 
 Double_t addsec=8640184.812866*tau+0.093104*pow(tau,2)-6.2e-6*pow(tau,3);
 sec=int(addsec);
 addsec-=double(sec);
 nsec=int(addsec*1.e9);
 addsec-=double(nsec)*1.e-9;
 psec=int(addsec*1.e12);
 sid.Add(0,sec,nsec,psec);
 
 sid.GetMJD(mjd,sec,nsec);
 psec=sid.GetPs();
 
 hh=sec/3600;
 sec=sec%3600;
 mm=sec/60;
 ss=sec%60;
 ns=nsec;
 ps=psec;
}

Double_t NcTimestamp::GetGMST()
{
 
 Int_t hh,mm,ss,ns,ps;
 
 GetGMST(hh,mm,ss,ns,ps);
 
 Double_t gst=Convert(hh,mm,ss,ns,ps);
 
 return gst;
}

Double_t NcTimestamp::GetGAST()
{
 
 Double_t da=Almanac();
 
 // Convert to fractional hours
 da/=3600.;
 
 Double_t gast=GetGMST()+da;
 
 while (gast<0)
 {
  gast+=24.;
 }
 while (gast>24.)
 {
  gast-=24.;
 }
 
 return gast;
}

Double_t NcTimestamp::GetLT(Double_t offset)
{
 
 // Current UT time in fractional hours
 Double_t h=GetUT();
 
 h+=offset;
 
 while (h<0)
 {
  h+=24.;
 }
 while (h>24)
 {
  h-=24.;
 }
 
 return h;
}

Double_t NcTimestamp::GetLAT(Double_t offset)
{
 
 // Current UT time in fractional hours
 Double_t h=GetUT();
 
 // Equation of Time
 Double_t eot; // LAT-LMT
 Almanac(0,0,0,0,"Sun",0,0,0,&eot,10);
 eot/=3600.; // Convert to fractional hours
 
 h+=offset+eot;
 
 while (h<0)
 {
  h+=24.;
 }
 while (h>24)
 {
  h-=24.;
 }
 
 return h;
}

Double_t NcTimestamp::GetLMST(Double_t offset)
{
 
 // Current GMST time in fractional hours
 Double_t h=GetGMST();
 
 h+=offset;
 
 while (h<0)
 {
  h+=24.;
 }
 while (h>24)
 {
  h-=24.;
 }
 
 return h;
}

Double_t NcTimestamp::GetLAST(Double_t offset)
{
 
 // Current GAST time in fractional hours
 Double_t h=GetGAST();
 
 h+=offset;
 
 while (h<0)
 {
  h+=24.;
 }
 while (h>24)
 {
  h-=24.;
 }
 
 return h;
}

void NcTimestamp::SetLT(Double_t dt,Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Int_t ss,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut)
{
 
 SetUT(y,m,d,hh,mm,ss,ns,ps,utc,leap,dut);
 Add(-dt);
}

void NcTimestamp::SetLT(Double_t dt,Int_t y,Int_t m,Int_t d,Int_t hh,Int_t mm,Double_t s,TString utc,Int_t leap,Double_t dut)
{
 
 SetUT(y,m,d,hh,mm,s,utc,leap,dut);
 Add(-dt);
}

void NcTimestamp::SetLT(Double_t dt,Int_t y,Int_t m,Int_t d,TString time,TString utc,Int_t leap,Double_t dut)
{
 
 SetUT(y,m,d,time,utc,leap,dut);
 Add(-dt);
}

void NcTimestamp::SetLT(Double_t dt,TString date,TString time,Int_t mode,TString utc,Int_t leap,Double_t dut)
{
 
 SetUT(date,time,mode,utc,leap,dut);
 Add(-dt);
}

void NcTimestamp::SetLT(Double_t dt,Int_t y,Int_t d,Int_t s,Int_t ns,Int_t ps,TString utc,Int_t leap,Double_t dut)
{
 
 SetUT(y,d,s,ns,ps,utc,leap,dut);
 Add(-dt);
}


Double_t NcTimestamp::GetJD(Double_t e,TString mode) const
{
 
 Double_t jd=0;
 
 if (mode=="J" || mode=="j") jd=(e-2000.0)*365.25+2451545.0;
 
 if (mode=="B" || mode=="b") jd=(e-1900.0)*365.242198781+2415020.31352;
 
 return jd;
}

Double_t NcTimestamp::GetMJD(Double_t e,TString mode) const
{
 
 Double_t mjd=GetJD(e,mode)-2400000.5;
 
 return mjd;
}

Double_t NcTimestamp::GetTJD(Double_t e,TString mode) const
{
 
 Double_t tjd=GetJD(e,mode)-2440000.5;
 
 return tjd;
}

Double_t NcTimestamp::Almanac(Double_t* dpsi,Double_t* deps,Double_t* eps,Double_t* dl,TString name,Double_t* el,Double_t* eb,Double_t* dr,Double_t* value,Int_t j)
{
 
 Double_t pi=acos(-1.);
 
 Double_t td;      // Time difference in fractional Julian days w.r.t. the start of J2000.
 Double_t tc;      // Time difference in fractional Julian centuries w.r.t. the start of J2000.
 Double_t tm;      // Time difference in fractional Julian millennia w.r.t. the start of J2000.
 const Int_t nvals=11;
 Double_t val[nvals]; // Array to hold the additional (orbital) parameters
 
 // Initialize the solar system body related values 
 if (el) *el=720;
 if (eb) *eb=720;
 if (dr) *dr=-1;
 for (Int_t i=0; i<nvals; i++)
 {
  val[i]=720;
 }
 val[0]=-1; // Semi major axis of the orbit
 val[1]=-1; // Eccentricity of the orbit
 val[10]=0; // Equation of Time, this will always be calculated correctly if requested
 
 // Initialize the requested additional observable value
 if (value)
 {
  *value=0;
  if (j>=0 && j<nvals) *value=val[j];
 }
 
 td=GetJD()-2451545.0;
 tc=td/36525.;
 tm=tc/10.;
 
 // Fundamental solar system variables (in arcseconds) w.r.t. the J2000.0 equinox.
 // The expressions are taken from the USNO circular 179.
 // epsilon : Mean obliquity of the ecliptic
 // l       : Mean anomaly of the Moon
 // lp      : Mean anomaly of the Sun
 // f       : Mean argument of latitude of the moon
 // d       : Mean elongation of the Moon from the Sun
 // om      : Mean longitude of the Moon's mean ascending node
 Double_t epsilon=84381.406-46.836769*tc-0.0001831*pow(tc,2)+0.00200340*pow(tc,3)-0.000000576*pow(tc,4)-0.0000000434*pow(tc,5);
 Double_t l=485868.249036+1717915923.2178*tc+31.8792*pow(tc,2)+0.051635*pow(tc,3)-0.00024470*pow(tc,4);
 Double_t lp=1287104.79305+129596581.0481*tc-0.5532*pow(tc,2)+0.000136*pow(tc,3)-0.00001149*pow(tc,4);
 Double_t f=335779.526232+1739527262.8478*tc-12.7512*pow(tc,2)-0.001037*pow(tc,3)+0.00000417*pow(tc,4);
 Double_t d=1072260.70369+1602961601.2090*tc-6.3706*pow(tc,2)+0.006593*pow(tc,3)-0.00003169*pow(tc,4);
 Double_t om=450160.398036-6962890.5431*tc+7.4722*pow(tc,2)+0.007702*pow(tc,3)-0.00005939*pow(tc,4);
 
 // General precession in longitude (in arcseconds) w.r.t. J2000.0
 // according to Jean Meeus Ch.21 (page 136)
 Double_t prec=5029.0966*tc+1.11113*pow(tc,2)-0.000006*pow(tc,3);
 
 if (eps) *eps=epsilon;
 if (dl) *dl=prec;
 
 // Convert to radians for use with goniometric functions
 Double_t fac=pi/(180.*3600.);
 epsilon*=fac;
 l*=fac;
 lp*=fac;
 f*=fac;
 d*=fac;
 om*=fac;
 
 //The IAU 2000A nutation series expansion.
 Double_t phi[28]={om,2.*(f-d+om),2.*(f+om),2.*om,lp,lp+2.*(f-d+om),l,
                   2.*f+om,l+2.*(f+om),2.*(f-d+om)-lp,2.*(f-d)+om,2.*(f+om)-l,2.*d-l,l+om,
                   om-l,2.*(f+d+om)-l,l+2.*f+om,2.*(f-l)+om,2.*d,2.*(f+d+om),2.*(f-d+om-lp),
                   2.*(d-l),2.*(l+d+om),l+2.*(f-d+om),2.*f+om-l,2.*l,2.*f,lp+om};
 Double_t s[28]={-17.2064161,-1.3170907,-0.2276413, 0.2074554, 0.1475877,-0.0516821, 0.0711159,
                  -0.0387298,-0.0301461, 0.0215829, 0.0128227, 0.0123457, 0.0156994, 0.0063110,
                  -0.0057976,-0.0059641,-0.0051613, 0.0045893, 0.0063384,-0.0038571, 0.0032481,
                  -0.0047722,-0.0031046, 0.0028593, 0.0020441, 0.0029243, 0.0025887,-0.0014053};
 Double_t sd[28]={-0.0174666,-0.0001675,-0.0000234, 0.0000207,-0.0003633, 0.0001226, 0.0000073,
                  -0.0000367,-0.0000036,-0.0000494, 0.0000137, 0.0000011, 0.0000010, 0.0000063,
                  -0.0000063,-0.0000011,-0.0000042, 0.0000050, 0.0000011,-0.0000001, 0.0000000,
                   0.0000000,-0.0000001, 0.0000000, 0.0000021, 0.0000000, 0.0000000,-0.0000025};
 Double_t cp[28]={ 0.0033386,-0.0013696, 0.0002796,-0.0000698, 0.0011817,-0.0000524,-0.0000872,
                   0.0000380, 0.0000816, 0.0000111, 0.0000181, 0.0000019,-0.0000168, 0.0000027,
                  -0.0000189, 0.0000149, 0.0000129, 0.0000031,-0.0000150, 0.0000158, 0.0000000,
                  -0.0000018, 0.0000131,-0.0000001, 0.0000010,-0.0000074,-0.0000066, 0.0000079};
 Double_t c[28]= { 9.2052331, 0.5730336, 0.0978459,-0.0897492, 0.0073871, 0.0224386,-0.0006750,
                   0.0200728, 0.0129025,-0.0095929,-0.0068982,-0.0053311,-0.0001235,-0.0033228,
                   0.0031429, 0.0025543, 0.0026366,-0.0024236,-0.0001220, 0.0016452,-0.0013870,
                   0.0000477, 0.0013238,-0.0012338,-0.0010758,-0.0000609,-0.0000550, 0.0008551};
 Double_t cd[28]={ 0.0009086,-0.0003015,-0.0000485, 0.0000470,-0.0000184,-0.0000677, 0.0000000,
                   0.0000018,-0.0000063, 0.0000299,-0.0000009, 0.0000032, 0.0000000, 0.0000000,
                   0.0000000,-0.0000011, 0.0000000,-0.0000010, 0.0000000,-0.0000011, 0.0000000,
                   0.0000000,-0.0000011, 0.0000010, 0.0000000, 0.0000000, 0.0000000,-0.0000002};
 Double_t sp[28]={ 0.0015377,-0.0004587, 0.0001374,-0.0000291,-0.0001924,-0.0000174, 0.0000358,
                   0.0000318, 0.0000367, 0.0000132, 0.0000039,-0.0000004, 0.0000082,-0.0000009,
                  -0.0000075, 0.0000066, 0.0000078, 0.0000020, 0.0000029, 0.0000068, 0.0000000,
                  -0.0000025, 0.0000059,-0.0000003,-0.0000003, 0.0000013, 0.0000011,-0.0000045};
 
 Double_t dp=0,de=0,da=0;
 for (Int_t i=0; i<28; i++)
 {
  dp+=(s[i]+sd[i]*tc)*sin(phi[i])+cp[i]*cos(phi[i]);
  de+=(c[i]+cd[i]*tc)*cos(phi[i])+sp[i]*sin(phi[i]);
 }
 
 da=dp*cos(epsilon)+0.00264096*sin(om)+0.00006352*sin(2.*om)
     +0.00001175*sin(2.*f-2.*d+3.*om)+0.00001121*sin(2.*f-2.*d+om)
     -0.00000455*sin(2.*f-2.*d+2.*om)+0.00000202*sin(2.*f+3.*om)+0.00000198*sin(2.*f+om)
     -0.00000172*sin(3.*om)-0.00000087*tc*sin(om);
 
 if (dpsi) *dpsi=dp;
 if (deps) *deps=de;
 
 // Convert to seconds
 da/=15.;

 // Determination of the Equation of Time via a recursive invokation //
 if (j==10)
 {
  Double_t xdpsi,xdeps,xeps,xlambda,xbeta;
  Almanac(&xdpsi,&xdeps,&xeps,0,"Sun",&xlambda,&xbeta);
 
  // Convert from arcsec to degrees
  xdpsi/=3600.;
  xdeps/=3600.;
  xeps/=3600.;
 
  // Correct for nutation to get the true values
  xeps+=xdeps;
  xlambda+=xdpsi;
 
  while (xlambda<0) { xlambda+=360.; }
  while (xlambda>360) { xlambda-=360.; }
 
  // Convert to radians for gonio
  Double_t epsr=xeps*pi/180.;
  Double_t lambdar=xlambda*pi/180.;
  Double_t betar=xbeta*pi/180.;
 
  // True Right Ascension of the Sun (see Ch.13 p.93 of J. Meeus) 
  Double_t x=sin(lambdar)*cos(epsr)-tan(betar)*sin(epsr);
  Double_t y=cos(lambdar);
  Double_t alpha=0;
  if (x || y) alpha=atan2(x,y)*180./pi;
 
  while (alpha<0) { alpha+=360.; }
  while (alpha>360) { alpha-=360.; }
 
  // Mean longitude of the Sun (see Ch.28 p.183 of J. Meeus)
  Double_t L0=280.4664567+360007.6982779*tm+0.03032028*pow(tm,2)+pow(tm,3)/49931.-pow(tm,4)/15300.-pow(tm,5)/2000000.;
 
  while (L0<0) { L0+=360.; }
  while (L0>360) { L0-=360.; }
 
  // Determine the Equation of Time in degrees (see Ch.28 p.183 of J. Meeus)
  Double_t eot=L0-0.0057183-alpha+xdpsi*cos(xeps*pi/180.);
 
  // The Equation of Time never exceeds 20 minutes (= 5 degrees)
  while (eot<-5) { eot+=360.; }
  while (eot>5) { eot-=360.; }
 
  // Convert the Equation of Time from degrees to seconds
  eot*=240.;
 
  val[10]=eot;
  if (value) *value=eot;
 }

 // Determination of the mean orbital elements and true ecliptic coordinates //
 // of a requested solar system body, for the mean equinox of the date.      //
 
 // The definitions and expressions are the ones used in the book of Jean Meeus
 // "Astronomical Algorithms" (2nd edition of August 2009), esp. chapters 31-33.
 
 // The various observables are :
 // a      : Semi major axis (in AU) of the orbit
 // e      : Eccentricity of the orbit
 // inc    : Inclination (in degrees) of the orbit with the ecliptic
 // omega  : Mean ecliptic longitude (in degrees) of the ascending node
 // lp     : Mean orbital longitude (in degrees) of the perihelion
 // l      : Mean orbital longitude (in degrees) of the body
 // omega2 : Orbital argument (in degrees) of the perihelion
 // m      : Mean orbital anomaly (in degrees) of the body
 // ec     : Equation of the center (in degrees)
 // nu     : True anomaly (in degrees) of the body
 // ltrue  : Heliocentric longitude of the body
 // btrue  : Heliocentric latitude of the body
 // r      : Distance between the body and the sun
 // lambda : Geocentric ecliptic longitude
 // beta   : Geocentric ecliptic latitude
 
 Double_t a=0;
 Double_t e=0;
 Double_t inc=0;
 Double_t omega=0;
 Double_t omega2=0;
 Double_t m=0;
 Double_t ec=0;
 Double_t nu=0;
 Double_t ltrue=0;
 Double_t btrue=0;
 Double_t r=0;
 Double_t lambda=0;
 Double_t beta=0;
 
 // Polynomial coefficients for a (in AU)of the 8 major planets 
 Double_t aa0[8]={0.387098310,0.723329820,1.000001018,1.523679342,5.202603209,9.554909192,19.218446062,30.110386869};
 Double_t aa1[8]={0,0,0,0,0.0000001913,-0.0000021390,-0.0000000372,-0.0000001663};
 Double_t aa2[8]={0,0,0,0,0,0.000000004,0.00000000098,0.00000000069};
 Double_t aa3[8]={0,0,0,0,0,0,0,0};
 // Polynomial coefficients for e of the 8 major planets 
 Double_t ea0[8]={0.20563175,0.00677192,0.01670863,0.09340065,0.04849793,0.05554814,0.04638122,0.00945575};
 Double_t ea1[8]={0.000020407,-0.000047765,-0.000042037,0.000090484,0.000163225,-0.000346641,-0.000027293,0.000006033};
 Double_t ea2[8]={-0.0000000283,0.0000000981,-0.0000001267,-0.0000000806,-0.0000004714,-0.0000006436,0.0000000789,0};
 Double_t ea3[8]={-0.00000000018,0.00000000046,0.00000000014,-0.00000000025,-0.00000000201,0.00000000340,0.00000000024,-0.00000000005};
 // Polynomial coefficients for inc (in degrees) of the 8 major planets 
 Double_t ia0[8]={7.004986,3.394662,0,1.849726,1.303267,2.488879,0.773197,1.769953};
 Double_t ia1[8]={0.0018215,0.0010037,0,-0.0006011,-0.0054965,-0.0037362,0.0007744,-0.0093082};
 Double_t ia2[8]={-0.00001810,-0.00000088,0,0.00001276,0.00000466,-0.00001519,0.00003749,-0.00000708};
 Double_t ia3[8]={0.000000056,-0.000000007,0,-0.000000007,-0.000000002,0.000000087,-0.000000092,0.000000027};
 // Polynomial coefficients for omega (in degrees) of the 8 major planets
 Double_t oa0[8]={48.330893,76.679920,0,49.558093,100.464407,113.665503,74.005957,131.784057};
 Double_t oa1[8]={1.1861883,0.9011206,0,0.7720959,1.0209774,0.8770880,0.5211278,1.1022039};
 Double_t oa2[8]={0.00017542,0.00040618,0,0.00001557,0.00040315,-0.00012176,0.00133947,0.00025952};
 Double_t oa3[8]={0.000000215,-0.000000093,0,0.000002267,0.000000404,-0.000002249,0.000018484,-0.000000637};
 // Polynomial coefficients for lp (in degrees) of the 8 major planets
 Double_t pa0[8]={77.456119,131.563703,102.937348,336.060234,14.331207,93.057237,173.005291,48.120276};
 Double_t pa1[8]={1.5564776,1.4022288,1.7195366,1.8410449,1.6126352,1.9637613,1.4863790,1.4262957};
 Double_t pa2[8]={0.00029544,-0.00107618,0.00045688,0.00013477,0.00103042,0.00083753,0.00021406,0.00038434};
 Double_t pa3[8]={0.000000009,-0.000005678,-0.000000018,0.000000536,-0.000004464,0.000004928,0.000000434,0.000000020};
 // Polynomial coefficients for l (in degrees) of the 8 major planets
 Double_t la0[8]={252.250906,181.979801,100.466457,355.433000,34.351519,50.077444,314.055005,304.348665};
 Double_t la1[8]={149474.0722491,58519.2130302,36000.7698278,19141.6964471,3036.3027748,1223.5110686,429.8640561,219.8833092};
 Double_t la2[8]={0.00030350,0.00031014,0.00030322,0.00031052,0.00022330,0.00051908,0.00030390,0.00030882};
 Double_t la3[8]={0.000000018,0.000000015,0.000000020,0.000000016,0.000000037,-0.000000030,0.000000026,0.000000018};
 
 TString names[10]={"Mercury","Venus","Earth","Mars","Jupiter","Saturn","Uranus","Neptune","Sun","Moon"};
 
 Int_t k=-1;
 Int_t geo=1;
 for (Int_t jbody=0; jbody<10; jbody++)
 {
  if (name.Contains("*")) geo=0;
  if (name.Contains(names[jbody].Data()))
  {
   k=jbody;
   break;
  }
 }
 
 if (k<0) return da; // Non-supported solar system body
 
 if (!geo && k==8) return da; // Request for heliocentric data of the Sun itself
 
 if (geo && k==2) return da; // Request for geocentric data of the Earth itself
 
 // In case geocentric data for the Sun are requested, the heliocentric data
 // of the Earth are used to construct the corresponding data for the Sun.
 Int_t sun=0;
 if (k==8)
 {
  k=2;
  sun=1;
 }
 
 // In case heliocentric data for the Moon are requested, the corresponding data
 // for the Earth are provided in view of negligible differences within the
 // accuracy of the algorithms used here.
 Int_t moon=0;
 if (k==9) moon=1;
 if (!geo && k==9) k=2;  
 
 lambda=0;
 beta=0;

 // Determination of the geocentric data for the Moon //
 
 if (geo && k==9)
 {
  // Low-precision geocentric ecliptic coordinates (in degrees) of the Moon.
  // Source : Astronomical Alamanac 2012 page D22.
  // Maximal errors : 0.3 degr. in lambda, 0.2 degr. in beta, 0.003 degr. in plax and 0.2 R_Earth in r
  // A more accurate method is the series expansion given in the book of Jean Meeus.
  lambda=218.32+481267.881*tc
         +6.29*sin((135.+477198.87*tc)*pi/180.)-1.27*sin((259.3-413335.36*tc)*pi/180.)
         +0.66*sin((235.7+890534.22*tc)*pi/180.)+0.21*sin((269.9+954397.74*tc)*pi/180.)
         -0.19*sin((357.5+35999.05*tc)*pi/180.)-0.11*sin((186.5+966404.03*tc)*pi/180.);
  beta=5.13*sin((93.3+483202.02*tc)*pi/180.)+0.28*sin((228.2+960400.89*tc)*pi/180.)
       -0.28*sin((318.3+6003.15*tc)*pi/180.)-0.17*sin((217.6-407332.21*tc)*pi/180.);
  Double_t plax=0.9508
                +0.0518*cos((135.+477198.87*tc)*pi/180.)+0.0095*cos((259.3-413335.36*tc)*pi/180.)
                +0.0078*cos((235.7+890534.22*tc)*pi/180.)+0.0028*cos((269.9+954397.74*tc)*pi/180.);
  r=1./sin(plax*pi/180.);
 
  // Convert r into km using an average Earth radius of 6367.45 km
  r*=6367.45;
 
  while (lambda<0) { lambda+=360.; }
  while (lambda>360) { lambda-=360.; }
 
  if (el) *el=lambda;
  if (eb) *eb=beta;
  if (dr) *dr=r;
 
  return da;
 }

 // Determination of the heliocentric data for the requested solar system body //
 
 a=0;
 e=0;
 inc=0;
 omega=0;
 l=0;
 lp=0;
 
 a=aa0[k]+aa1[k]*tc+aa2[k]*pow(tc,2)+aa3[k]*pow(tc,3);
 e=ea0[k]+ea1[k]*tc+ea2[k]*pow(tc,2)+ea3[k]*pow(tc,3);
 inc=ia0[k]+ia1[k]*tc+ia2[k]*pow(tc,2)+ia3[k]*pow(tc,3);
 omega=oa0[k]+oa1[k]*tc+oa2[k]*pow(tc,2)+oa3[k]*pow(tc,3);
 lp=pa0[k]+pa1[k]*tc+pa2[k]*pow(tc,2)+pa3[k]*pow(tc,3);
 l=la0[k]+la1[k]*tc+la2[k]*pow(tc,2)+la3[k]*pow(tc,3);
 
 while (omega<0) { omega+=360.; }
 while (omega>360) { omega-=360.; }
 
 m=l-lp;
 
 while (m<0) { m+=360.; }
 while (m>360) { m-=360.; }
 
 ec=(2.*e-(pow(e,3)/4.)+(5.*pow(e,5)/96.))*sin(m*pi/180.)+((5.*pow(e,2)/4.)-(11.*pow(e,4)/24.))*sin(2.*m*pi/180.)
    +((13.*pow(e,3)/12.)-(43.*pow(e,5)/64.))*sin(3.*m*pi/180.)+(103.*pow(e,4)*sin(4.*m*pi/180.)/96.)
    +(1097.*pow(e,5)*sin(5.*m*pi/180.)/960.);
 ec*=180./pi;
 
 omega2=lp-omega;
 
 while (omega2<0) { omega2+=360.; }
 while (omega2>360) { omega2-=360.; }
 
 while (lp<0) { lp+=360.; }
 while (lp>360) { lp-=360.; }
 
 while (l<0) { l+=360.; }
 while (l>360) { l-=360.; }
 
 nu=m+ec;
 
 while (nu<0) { nu+=360.; }
 while (nu>360) { nu-=360.; }
 
 // Store the orbital parameters in the additional values array
 if (!sun && !moon)
 {
  val[0]=a;
  val[1]=e;
  val[2]=inc;
  val[3]=omega;
  val[4]=lp;
  val[5]=l;
  val[6]=omega2;
  val[7]=m;
  val[8]=ec;
  val[9]=nu;
 }
 
 // Make requested value available
 if (value)
 {
  *value=0;
  if (j>=0 && j<nvals) *value=val[j];
 }
 
 r=a*(1.-e*e)/(1.+e*cos(nu*pi/180.));
 
 // Use sine rule to obtain the latitude in radians
 Double_t sinb=sin(inc*pi/180.)*sin((l-omega+ec)*pi/180.);
 btrue=asin(sinb);
 
 ltrue=omega;
 // Use Neper's rule to obtain the extra term of the longitude
 Double_t arg=l-omega+ec;
 while (arg<0) { arg+=360.; }
 while (arg>360) { arg-=360.; }
 Double_t extra=0;
 if (cos(btrue))
 {
  Double_t cosl=cos(arg*pi/180.)/cos(btrue);
  extra=acos(cosl)*180./pi;
  if (arg>180) extra=-extra;
 }
 
 btrue*=180./pi;
 ltrue+=extra;
 
 // Convert heliocentric Earth data into geocentric Sun data if requested
 if (sun)
 {
  btrue=-btrue;
  ltrue+=180.;
 }
 
 while (ltrue<0) { ltrue+=360.; }
 while (ltrue>360) { ltrue-=360.; } 
 
 if (el) *el=ltrue;
 if (eb) *eb=btrue;
 if (dr) *dr=r;
 
 if (!geo || sun) return da; // Heliocentric (or geocentric Sun) coordinates were requested

 // Convert into geocentric ecliptic coordinates //
 
 // The algorithm used here is the one outlined in Ch. 33 of the book of Jean Meeus.
 // In view of the accuracy of the current algorithm, the effects of light-time
 // and aberration are not taken into account here.
 
 // Determine the heliocentric coordinates of the Earth
 // via recursive invokation of this memberfunction
 Double_t l0,b0,r0;
 Almanac(0,0,0,0,"Earth*",&l0,&b0,&r0);
 
 Double_t x=r*cos(btrue*pi/180.)*cos(ltrue*pi/180.)-r0*cos(b0*pi/180.)*cos(l0*pi/180.);
 Double_t y=r*cos(btrue*pi/180.)*sin(ltrue*pi/180.)-r0*cos(b0*pi/180.)*sin(l0*pi/180.);
 Double_t z=r*sin(btrue*pi/180.)-r0*sin(b0*pi/180.);
 
 lambda=atan2(y,x)*180./pi;
 beta=atan2(z,sqrt(x*x+y*y))*180./pi;
 r=sqrt(x*x+y*y+z*z);
 
 while (lambda<0) { lambda+=360.; }
 while (lambda>360) { lambda-=360.; }
 
 if (el) *el=lambda;
 if (eb) *eb=beta;
 if (dr) *dr=r;
 
 return da;
}

void NcTimestamp::SetEpoch(Double_t e,TString mode,TString utc,Int_t leap,Double_t dut)
{
 
 Double_t jd=GetJD(e,mode);
 SetJD(jd,utc,leap,dut);
}

Double_t NcTimestamp::GetEpoch(TString mode)
{
 
 Double_t e=0;
 if (mode=="B" || mode=="b") e=GetBE();
 if (mode=="J" || mode=="j") e=GetJE();
 return e;
}

TString NcTimestamp::GetDayTimeString(TString mode,Int_t ndig,Double_t offset,TString* date,TString* time,Bool_t full)
{
 
 Bool_t set=kFALSE;
 
 TString sdate=mode;
 sdate+=" information unavailable";
 
 TString stime=mode;
 stime+=" information unavailable";
 
 TString daytime=mode;
 daytime+=" information unavailable";
 
 TString month[12]={"Jan","Feb","Mar","Apr","May","Jun","Jul","Aug","Sep","Oct","Nov","Dec"};
 TString day[7]={"Mon","Tue","Wed","Thu","Fri","Sat","Sun"};
 
 UInt_t y=0;
 UInt_t m=0;
 UInt_t d=0;
 UInt_t wd=0;
 Int_t hh,mm,ss,ns,ps;
 Double_t s=0;
 ULong64_t sfrac=0;
 Double_t gat=0;
 Double_t gast=0;
 Bool_t bdate=kFALSE;
 
 Int_t mjd,mjsec,mjns;
 GetMJD(mjd,mjsec,mjns);
 
 // Check whether UT1 is the main reference
 if (mode=="UT1" && fUtc) mode="UT";
 
 // Check whether UTC is the main reference
 if (mode=="UTC" && !fUtc) mode="UT";
 
 // UT related information :  UT, GAT, GMST and GAST date and time
 if (mode=="UT" || mode=="GAT" || mode=="GMST" || mode=="GAST")
 {
  if (mjd>=40587 && (mjd<65442 || (mjd==65442 && mjsec<8047)))
  {
   GetDate(kTRUE,0,&y,&m,&d);
   wd=GetDayOfWeek(kTRUE,0);
   bdate=kTRUE;
  }
  if (mode=="UT") GetUT(hh,mm,ss,ns,ps);
  if (mode=="GMST") GetGMST(hh,mm,ss,ns,ps);
  if (mode=="GAT")
  {
   // Obtain GAT as LAT without offset
   gat=GetLAT(0);
   Convert(gat,hh,mm,ss,ns,ps);
  }
  if (mode=="GAST")
  {
   gast=GetGAST();
   Convert(gast,hh,mm,ss,ns,ps);
  }
 
  set=kTRUE;
 }
 
 // TAI related information : TAI, UTC, GPS and TT date and time
 if (mode=="TAI" || mode=="UTC" || mode=="GPS" || mode=="TT")
 {
  if (fUtc && fUtc!=-3)
  {
   // A dummy timestamp is used to obtain the TAI corresponding date indicator
   NcTimestamp tx;
   tx.SetMJD(fTmjd,fTsec,fTns,fTps);
 
   // Convert to the corresponding UTC, GPS or TT timestamps
   if (mode=="UTC") tx.AddSec(-fLeap);
   if (mode=="GPS") tx.AddSec(-19);
   if (mode=="TT") tx.AddSec(32.184);
 
   Int_t tmjd,tsec,tns;
   tx.GetMJD(tmjd,tsec,tns); 
 
   if (tmjd>=40587 && (tmjd<65442 || (tmjd==65442 && tsec<8047)))
   {
    tx.GetDate(kTRUE,0,&y,&m,&d);
    wd=tx.GetDayOfWeek(kTRUE,0);
    bdate=kTRUE;
   }
   GetTAI(hh,mm,ss,ns,ps,mode);
 
   set=kTRUE;
  }
 }
 
 // Local time information
 if (mode=="LMT" || mode=="LAT" || mode=="LMST" || mode=="LAST")
 {
  // Determine the new date by including the offset
  NcTimestamp t2(*this);
  t2.Add(offset);
  Int_t mjd2,mjsec2,mjns2;
  t2.GetMJD(mjd2,mjsec2,mjns2);
  if (mjd2>=40587 && (mjd2<65442 || (mjd2==65442 && mjsec2<8047)))
  {
   t2.GetDate(kTRUE,0,&y,&m,&d);
   wd=t2.GetDayOfWeek(kTRUE,0);
   bdate=kTRUE;
  }
  // Determine the local time by including the offset w.r.t. the original timestamp
  Double_t hlt=0;
  if (mode=="LMT") hlt=GetLT(offset);
  if (mode=="LAT") hlt=GetLAT(offset);
  if (mode=="LMST") hlt=GetLMST(offset);
  if (mode=="LAST") hlt=GetLAST(offset);
  Convert(hlt,hh,mm,ss,ns,ps);
 
  set=kTRUE;
 }
 
 // No match found with requested date/time system
 if (!set)
 {
  if (date) *date=sdate;
  if (time) *time=stime;
  return daytime;
 }
 
 // Create the date and time info string
 sdate="";
 stime="";
 
 // The date string
 if (bdate)
 {
  if (full)
  {
   sdate.Form("%-s %02i-%-s-%-i",day[wd-1].Data(),d,month[m-1].Data(),y);
  }
  else
  {
   sdate.Form("%02i-%02i-%-i",d,m,y);
  }
 }
 
 // The time string
 if (hh<10) stime+="0";
 stime+=hh;
 stime+=":";
 if (mm<10) stime+="0";
 stime+=mm;
 stime+=":";
 if (ss<10) stime+="0";
 stime+=ss;
 if (ndig>0)
 {
  stime+=".";
  s=double(ns)*1e-9+double(ps)*1e-12;
  s*=pow(10.,ndig);
  sfrac=ULong64_t(s);
  ULong64_t isfrac=pow(10,ndig-1);
  while (sfrac<isfrac)
  {
   stime+="0";
   isfrac/=10;
  }
  if (sfrac) stime+=sfrac;
 }
 
 // Construct the combined date/time string
 daytime=sdate;
 daytime+=" ";
 daytime+=stime;
 
 // Add the time system indicator
 if (mode=="UT")
 {
  mode="UT1";
  if (!fUtc) mode="UTC";
 }
 if (sdate=="") sdate="---";
 if (full)
 {
  sdate+=" ";
  sdate+=mode;
  stime+=" ";
  stime+=mode;
  daytime+=" ";
  daytime+=mode;
 }
 
 if (date) *date=sdate;
 if (time) *time=stime;
 
 return daytime;
}

void NcTimestamp::SetSystemTime()
{
 
 // Retrieve the current system clock time
 Set();
 
 // Store this time as UTC
 FillJulian();
 Int_t ret=SetUTCparameters("A",0,0);
 
 // Use the UTC parameters (if any) to determine UT1
 if (ret && ret!=-3)
 {
  AddSec(fDut); // Correct UTC for dUT=UT1-UTC to obtain UT1 as main reference
 }
 else
 {
  SetUTCparameters("N",0,0); // Keep the system clock time as main (UTC) reference
 }
 FillTAI();
}

Bool_t NcTimestamp::IsUT1() const
{
 
 Bool_t ut1=kTRUE;
 if (!fUtc) ut1=kFALSE;
 
 return ut1;
}