Fix computation of Earth vector and get time difference in TT timesystem

Jürgen Knödlseder
1 parent 810e8ce5
Showing 3 changed files with 26 additions and 13 deletions
NEWS
src/obs/GEphemerides.cpp
src/obs/GPulsar.cpp
 New Features and Important Changes in GammaLib 2.0.0
  
-14 January 2022
+17 January 2022
  
  
 1. Introduction
@@ -53,6 +53,8 @@ The following methods have been added:
 - GObservation::model_grad(GModel&, GModelPar&)
 - GObservation::grad_step_size()
 - GGti::reduce(GGti&)
+- GTime::leap_seconds()
+- GTime::utc2tt()
 - GPhases::length()
 - GMatrixBase::is_empty()
 - angle(GVector&, GVector&)
@@ -280,7 +282,9 @@ of functions (#3328).
 ----------------------
 Add GEphemerides class to deal with JPL ephemerides and specifically the
 DE200 ephemerides. A JPL DE200 ephemerides file was added to a refdata
-direction, the file is installed in $GAMMALIB/share/refdata (#3973).
+direction, the file is installed in $GAMMALIB/share/refdata. Add further-
+more GTime::leap_seconds() and GTime::utc2tt() methods to support time
+system conversions in GEphemerides and GPulsar classes (#3973).
  
 Add GPhases::length() method that returns the total length of the phase
 intervals (#3973).
@@ -266,6 +266,9 @@ void GEphemerides::load(const GFilename&amp; filename)
  *
  * Get ephemeris vector and TBD-TT value for a given time. Information is
  * only returned for pointers that are not NULL.
+ *
+ * The code was inspired from the ftools routine xtereadeph.f and the COMPASS
+ * routine bvceph.f.
  ***************************************************************************/
 void GEphemerides::ephemeris(const GTime& time,
                              GVector*     rce,
@@ -276,10 +279,16 @@ void GEphemerides::ephemeris(const GTime&amp; time,
     // Fetch ephemerides
     const_cast<GEphemerides*>(this)->fetch_data();
  
-    // Compute jd_utc as int and dt as day fraction, between -0.5 and 0.5
+    // Compute jd_utc as int
     double jd     = time.jd("UTC");
     int    jd_utc = int(jd);
-    double dt     = jd - double(jd_utc);
+
+    // Compute dt in TT as day fraction, between -0.5 and 0.5. While the
+    // time system in the xtereadeph.f routine is not specified, the
+    // bvceph.f routine explicitly adds the UTC->TT conversion coefficient
+    // to the delta
+    //double dt     = jd - double(jd_utc);
+    double dt = time.jd("TT") - double(int(time.jd("TT")));
     if (dt > 0.5) {
         jd_utc += 1;
         dt     -= 1.0;
@@ -326,8 +335,8 @@ void GEphemerides::ephemeris(const GTime&amp; time,
     }
  
     // Get Sun vector
-    if (rce != NULL) {
-        *rce = m_sun[index];
+    if (rcs != NULL) {
+        *rcs = m_sun[index];
     }
  
     // Get TBD-TT (seconds)
@@ -380,7 +389,7 @@ double GEphemerides::geo2ssb(const GTime&amp; time, const GSkyDir&amp; srcdir) const
     // Special handling if sky direction is inside the Sun
     /*
     if (cth + 1.0 < 0.5 * sunsiz * sunsiz) {
-        
+
     }
     */
  
@@ -435,7 +435,7 @@ void GPulsar::init_members(void)
     // Initialise members
     m_name.clear();
     m_ephemerides.clear();
-    
+
     // Return
     return;
 }
@@ -785,7 +785,7 @@ void GPulsar::load_psrtime(const GFilename&amp; filename, const std::string&amp; name)
                 throw GException::invalid_argument(G_LOAD_PSRTIME, msg);
             }
         }
-        
+
         // Save the pulsar name
         psrname = psrb;
         m_name  = "PSR B" + psrb;
@@ -808,7 +808,7 @@ void GPulsar::load_psrtime(const GFilename&amp; filename, const std::string&amp; name)
         double ra  = (ra_h + ra_m/60.0 + ra_s/3600.0) * 15.0;
         double dec = (dec_d < 0.0) ? -(-dec_d + dec_m/60.0 + dec_s/3600.0)
                                    : (dec_d + dec_m/60.0 + dec_s/3600.0);
-       
+
         // Set sky direction
         GSkyDir dir;
         dir.radec_deg(ra, dec);
@@ -831,14 +831,14 @@ void GPulsar::load_psrtime(const GFilename&amp; filename, const std::string&amp; name)
         toa.mjd(t0geo, "UTC");
  
         // Compute phase of first pulse according to formulae given in
-        // COM-RP-DOL-DRG-065.Note that the phase is negative, see for
+        // COM-RP-DOL-DRG-065. Note that the phase is negative, see for
         // example Eq. (3) in Yan et al. (2017), ApJ, 845, 119
         double geo2ssb     = ephemerides.geo2ssb(toa, dir);
         double utc2tt      = toa.utc2tt();
         double dt          = (t0geo - t0nom) * gammalib::sec_in_day + geo2ssb + utc2tt;
         const double c1    = 1.0/2.0;
         const double c2    = 1.0/6.0;
-        double       phase = -((f0 + (f1 * c1 + f2 * dt * c2) * dt) * dt);
+        double       phase = ((f0 + (f1 * c1 + f2 * dt * c2) * dt) * dt);
         phase             -= floor(phase);
  
         // Allocate ephemeris
@@ -851,7 +851,7 @@ void GPulsar::load_psrtime(const GFilename&amp; filename, const std::string&amp; name)
         ephemeris.tstop(tstop);
         ephemeris.timesys("UTC");  // psrtime ephemerides are in UTC
         ephemeris.t0(t0);
-        ephemeris.phase(phase);
+        ephemeris.phase(-phase);
         ephemeris.f0(f0);
         ephemeris.f1(f1);
         ephemeris.f2(f2);