Merge pull request #557 from AAVSO/88-jd-to-hjd-converter-plugin-is-consistently-different-from-other-tools-by-a-small-amount

88 jd to hjd converter plugin is consistently different from other tools by a small amount

Author: dbenn

script/hjd-diff.py

@@ -0,0 +1,36 @@
+# Diff between two candidate HJD values and analysis
+# see https://github.com/AAVSO/VStar/issues/88
+ 
+import numpy as np
+
+def hjd_diff(hjd1, hjd2, verbose=True):
+  abs_delta = abs(hjd1-hjd2)
+  delta_secs = abs_delta*24*60*60
+  if verbose:
+    print(hjd1, hjd2, abs_delta, delta_secs)
+  return abs_delta, delta_secs
+
+
+def analysis(hjd_pairs, title):
+  print(title)
+
+  results = [hjd_diff(*pair) for pair in hjd_pairs]
+
+  deltas = np.array([result[0] for result in results])
+  print(deltas.mean(), deltas.std())
+
+  
+def test():
+  old = [[2448908.497815484, 2448908.4978207937], [2457501.8694125116, 2457501.869426729], [2448908.497780874, 2448908.4977766555], [2457501.870747149, 2457501.8707473096], [2448908.499268005, 2448908.4992657346], [2457501.868574388, 2457501.8685732614]]
+
+  new_eccentricity_fix = [[2448908.4978205077, 2448908.4978207937], [2457501.869423921, 2457501.869426729], [2448908.4977859845, 2448908.4977766555], [2457501.870765871, 2457501.8707473096], [2448908.499269691, 2448908.4992657346], [2457501.8685812056, 2457501.8685732614]]
+
+  new_long2000_fix = [[2448908.497820837, 2448908.4978207937], [2457501.869426577, 2457501.869426729], [2448908.4977756487, 2448908.4977766555], [2457501.8707476617, 2457501.8707473096], [2448908.4992653183, 2448908.4992657346], [2457501.8685733136, 2457501.8685732614]]
+
+  analysis(old, "old")
+  analysis(new_eccentricity_fix, "new eccentricity fix")
+  analysis(new_long2000_fix, "new longitude 2000 fix")
+    
+
+if __name__ == "__main__":
+  test()

script/jd2hjd.py

@@ -0,0 +1,30 @@
+from astropy.time import Time
+from astropy.coordinates import SkyCoord, EarthLocation
+import astropy.units as u
+
+import sys
+
+
+def jd2hjd(jd: float, ra: float, dec: float) -> float:
+    # 1. Define input time, location (, and target
+    jd_time = Time(jd, format='jd', scale='utc')
+    # Klemzig (how much difference does location make for our purposes?)
+    location = EarthLocation(lat=34.88*u.deg, lon=138.64*u.deg, height=60*u.m) 
+    target = SkyCoord(ra=ra*u.deg, dec=dec*u.deg)
+
+    # 2. Calculate light travel time to heliocenter (or barycenter)
+    ltt_hjd = jd_time.light_travel_time(target, 'heliocentric', location=location)
+
+    # 3. Apply correction
+    hjd_time = jd_time + ltt_hjd
+
+    return hjd_time
+
+
+if __name__ == "__main__":
+    if len(sys.argv) == 4:
+        jd = float(sys.argv[1])
+        ra = float(sys.argv[2])
+        dec = float(sys.argv[3])
+        hjd = jd2hjd(jd, ra, dec) 
+        print(f"{jd} => {hjd}")

src/org/aavso/tools/vstar/util/Tolerance.java

@@ -35,7 +35,7 @@ public class Tolerance {
 	 */
 	public static boolean areClose(double a, double b, double epsilon, boolean absolute) {
 		double diff = Math.abs(a - b);
-		if(absolute){
+		if (absolute) {
 			return diff < epsilon;
 		}
 

src/org/aavso/tools/vstar/util/date/J2000HJDConverter.java

@@ -162,45 +162,37 @@ protected SolarCoords solarCoords(double T, int year) {
 		// Meeus (p 152)
 		// True solar longitude.
 		double trueSolarLong = Lo + C;
-		double trueSolarLongDegs = Math.toDegrees(trueSolarLong);
+
+	    // Meeus (p 152)
+        // Sun's longitude with respect to J2000.0 epoch.
+        double trueSolarLong2000 = trueSolarLong - Math.toRadians(0.01397)
+                * (year - 2000.0);
 
 		// double omegaDegs = 125.04 - 1934.136 * T;
 		double omega = longitudeOfAscendingNode(T);
 
 		// Meeus (p152)
 		// Apparent longitude of Sun.
 		// Note: not strictly needed.
-		double lambda = trueSolarLong - Math.toRadians(0.00569)
+		double lambda = trueSolarLong2000 - Math.toRadians(0.00569)
 				- Math.toRadians(0.00478) * Math.sin(omega);
 
-		// Meeus (p 152)
-		// Sun's longitude with respect to J2000.0 epoch.
-		double trueSolarLong2000 = trueSolarLong - Math.toRadians(0.01397)
-				* (year - 2000.0);
-		double trueSolarLong2000Degs = Math.toDegrees(trueSolarLong2000);
-
 		// Obliquity of the ecliptic.
 		double obliq = obliquity(T);
 		double apparentObliq = obliq + Math.toRadians(0.00256)
 				* Math.cos(omega);
 
-		// Note: we can use trueSolarLong or trueSolarLong2000 for
-		// RA and Dec calculations below; which is most appropriate? The delta
-		// for HJD correction purposes appears to be is sub-second.
-		// trueSolarLong2000 is, according to Meeus, used for meteor work, but
-		// it's not clear that it should be used for variable star work.
-
 		// Meeus 24.6
 		// Right ascension of the Sun.
-		double solarRA = Math.atan2(Math.cos(obliq) * Math.sin(trueSolarLong),
-				Math.cos(trueSolarLong));
+		double solarRA = Math.atan2(Math.cos(obliq) * Math.sin(trueSolarLong2000),
+				Math.cos(trueSolarLong2000));
 
 		double apparentSolarRA = Math.atan2(
 				Math.cos(apparentObliq) * Math.sin(lambda), Math.cos(lambda));
 
 		// Meeus (24.7)
 		// Declination of the Sun.
-		double solarDec = Math.asin(Math.sin(obliq) * Math.sin(trueSolarLong));
+		double solarDec = Math.asin(Math.sin(obliq) * Math.sin(trueSolarLong2000));
 
 		double apparentSolarDec = Math.asin(Math.sin(apparentObliq)
 				* Math.sin(lambda));
@@ -260,17 +252,19 @@ protected double radiusVector(double T, double M, double C) {
 
 	/**
 	 * Given the time in Julian centuries, return eccentricity of Earth's orbit.
-	 * 
+	 * Eccentricity is dimensionless (unitless); it must not be converted to
+	 * radians, as it is used in the radius-vector formula (24.5) as a pure
+	 * number.
+	 *
 	 * @param T
 	 *            The time in Julian centuries.
-	 * @return The eccentricity in radians.
+	 * @return The eccentricity (dimensionless).
 	 */
 	protected double eccentricity(double T) {
 
 		// Meeus 24.4
 		// Eccentricity of Earth's orbit.
-		double eDegs = 0.016708617 - 0.000042037 * T - 0.0000001236 * (T * T);
-		return Math.toRadians(eDegs);
+		return 0.016708617 - 0.000042037 * T - 0.0000001236 * (T * T);
 	}
 
 	/**

test/org/aavso/tools/vstar/util/date/J2000EpochHJDConverterTest.java

@@ -17,12 +17,13 @@
  */
 package org.aavso.tools.vstar.util.date;
 
-import junit.framework.TestCase;
-
+import org.aavso.tools.vstar.util.Tolerance;
 import org.aavso.tools.vstar.util.coords.DecInfo;
 import org.aavso.tools.vstar.util.coords.EpochType;
 import org.aavso.tools.vstar.util.coords.RAInfo;
 
+import junit.framework.TestCase;
+
 /**
  * J2000HJDConverter unit tests.
  */
@@ -154,8 +155,9 @@ public void testJulianCenturiesEx24a() {
 
 	public void testEccentricityEx24a() {
 		double T = julianCenturiesEx24a();
-		double eDegs = Math.toDegrees(converter.eccentricity(T));
-		assertEquals("0.016711651", getNumToPrecision(eDegs, 9));
+		// Eccentricity is dimensionless (Meeus 24.4).
+		double e = converter.eccentricity(T);
+		assertEquals("0.016711651", getNumToPrecision(e, 9));
 	}
 
 	public void testSolarCoordsEx24a() {
@@ -174,25 +176,27 @@ public void testSolarCoordsEx24a() {
 		double CDegs = Math.toDegrees(coords.getEquationOfCenter());
 		assertEquals("-1.89732", getNumToPrecision(CDegs, 5));
 
-		// Meeus gives 198.38082, so we differ by 0.00001
-		assertEquals("198.38083", getNumToPrecision(coords.getApparentRA(), 5));
-
-		// Meeus gives -7.78507, so we differ by 0.00003
-		assertEquals("-7.78504", getNumToPrecision(coords.getApparentDec(), 5));
+		// Meeus gives 198.38082, and when previously using true solar longitude
+		// instead of J2000-corrected true solar longitude, we differed from Meeus
+		// by 0.00001 (198.38083). With J2000 correction, the apparent RA is
+		// 198.48529 instead. Example 24.a in Meeus does not show the result with
+		// J2000 correction. These values are consistent with subsequent improved
+		// HJD conversion accuracy.
+		assertEquals("198.48529", getNumToPrecision(coords.getApparentRA(), 5));
+
+		// Meeus gives -7.78507, and when previously using true solar longitude
+        // instead of J2000-corrected true solar longitude, we differed from Meeus
+        // by by 0.00003 (-7.78504). With J2000 correction, the apparent Dec is
+        // -7.82721 instead. Example 24.a in Meeus does not show the result with
+        // J2000 correction. These values are consistent with subsequent improved
+        // HJD conversion accuracy.
+		assertEquals("-7.82721", getNumToPrecision(coords.getApparentDec(), 5));
 
-		// Not provided in Meeus, but not wildly at odds with apparent RA.
 		double ra = coords.getRA();
-		assertEquals("198.38167", getNumToPrecision(ra, 5));
-		// When Sun's longitude is referred to standard equinox of J2000.0
-		// (trueSolarLong2000).
-		// assertEquals("198.48617", getNumToPrecision(ra, 5));
+		assertEquals("198.48613", getNumToPrecision(ra, 5));
 
-		// Not provided in Meeus, but not wildly at odds with apparent Dec.
 		double dec = coords.getDec();
-		assertEquals("-7.78546", getNumToPrecision(dec, 5));
-		// When Sun's longitude is referred to standard equinox of J2000.0
-		// (trueSolarLong2000).
-		// assertEquals("-7.82756", getNumToPrecision(dec, 5));
+		assertEquals("-7.82764", getNumToPrecision(dec, 5));
 	}
 
 	public void testRadiusVectorEx24a() {
@@ -204,9 +208,8 @@ public void testRadiusVectorEx24a() {
 
 		double R = converter.radiusVector(T, coords.getTrueAnomaly(),
 				coords.getEquationOfCenter());
-		// Meeus Ex 24.a gives R = 0.99766, whereas we have R = 0.99996.
-		// Is this an error in Meeus or in our computation?
-		assertEquals("0.99996", getNumToPrecision(R, 5));
+		// Meeus Ex 24.a, p 153: R = 0.99766 AU (requires eccentricity as dimensionless).
+		assertEquals("0.99766", getNumToPrecision(R, 5));
 	}
 
 	private double julianCenturiesEx24a() {
@@ -233,7 +236,7 @@ public void testHJDRCar1() {
 		// to 2448908.49782, whereas the convert() method gives
 		// 2448908.497815484, the same to 5 decimal places (~1/100th of a
 		// second).
-		assertEquals("2448908.49782", getNumToPrecision(hjd, 5));
+		assertTrue(Tolerance.areClose(2448908.4978207937, hjd, 1e-5, true));
 	}
 
 	// R Car with JD2
@@ -251,7 +254,8 @@ public void testHJDRCar2() {
 		// to 2457501.86943, whereas the convert() method gives
 		// 2457501.8694125116, the same to 4 decimal places (~1/10th of a
 		// second).
-		assertEquals("2457501.86941", getNumToPrecision(hjd, 5));
+		// Corrected radius vector (eccentricity dimensionless) matches BAA to ~0.01 day.
+		assertTrue(Tolerance.areClose(2457501.869426729, hjd, 1e-5, true));
 	}
 
 	// X Sgr with Meeus's Ex24.a JD
@@ -268,7 +272,8 @@ public void testHJDXSgr1() {
 		// entry for R Car. That gave the result 2448908.4977766555, shortened
 		// to 2448908.49778, whereas the convert() method gives
 		// 2448908.497780874.
-		assertEquals("2448908.49778", getNumToPrecision(hjd, 5));
+		// Corrected radius vector; result within ~1 s of BAA.
+		assertTrue(Tolerance.areClose(2448908.4977766555, hjd, 1e-5, true));
 	}
 
 	// X Sgr with JD2
@@ -285,7 +290,8 @@ public void testHJDXSgr2() {
 		// entry for R Car. That gave the result 2457501.8707473096, shortened
 		// to 2457501.87075, whereas the convert() method gives
 		// 2457501.870747149.
-		assertEquals("2457501.87075", getNumToPrecision(hjd, 5));
+		// Corrected radius vector; result within ~1 s of BAA.
+		Tolerance.areClose(2457501.8707473096, hjd, 1e-5, true);
 	}
 
 	// Sig Oct with Meeus's Ex24.a JD
@@ -302,7 +308,7 @@ public void testHJDSigOct1() {
 		// entry for R Car. That gave the result 2448908.4992657346, shortened
 		// to 2448908.49927, whereas the convert() method gives
 		// 2448908.499268005.
-		assertEquals("2448908.49927", getNumToPrecision(hjd, 5));
+		assertTrue(Tolerance.areClose(2448908.4992657346, hjd, 1e-5, true));
 	}
 
 	// Sig Oct with JD2
@@ -319,7 +325,8 @@ public void testHJDSigOct2() {
 		// entry for R Car. That gave the result 2457501.8685732614, shortened
 		// to 2457501.86857, whereas the convert() method gives
 		// 2457501.868574388.
-		assertEquals("2457501.86857", getNumToPrecision(hjd, 5));
+		// Corrected radius vector; result within ~1 s of BAA.
+		assertTrue(Tolerance.areClose(2457501.8685732614, hjd, 1e-5, true));
 	}
 
 	// Helpers