Tweak refrated RA/Dec calculation in equ2hor()

Smithsonian · Feb 22, 2024 · 609917a · 609917a
1 parent d01538b
commit 609917a
Show file tree

Hide file tree

Showing 3 changed files with 35 additions and 28 deletions.
diff --git a/src/novas.c b/src/novas.c
@@ -2200,11 +2200,10 @@ int equ2hor(double jd_ut1, double ut1_to_tt, enum novas_accuracy accuracy, doubl
   if(ref_option) {
     // Get refraction in zenith distance.
     const double zd0 = *zd;
-    const double refr = refract_astro(location, ref_option, *zd);
-    *zd -= refr;
+    const double refr = refract_astro(location, ref_option, zd0);
 
     // Apply refraction to celestial coordinates of object.
-    if((refr > 0.0) && (*zd > 3.0e-4)) {
+    if(refr > 0.0) {
       // Shift position vector of object in celestial system to account
       // for refraction (see USNO/AA Technical Note 1998-09).
       const double sinzd = sin(*zd * DEGREE);
@@ -2214,6 +2213,8 @@ int equ2hor(double jd_ut1, double ut1_to_tt, enum novas_accuracy accuracy, doubl
 
       int j;
 
+      *zd -= refr;
+
       // Compute refracted position vector.
       for(j = 3; --j >= 0;)
         pr[j] = ((p[j] - coszd0 * uz[j]) / sinzd0) * sinzd + uz[j] * coszd;

diff --git a/test/reference/equ2hor-refract.out b/test/reference/equ2hor-refract.out
@@ -3,38 +3,38 @@
 
 -10000.0 16-20      S2 O1 A0:   115.693174   103.225975    16.000000   -20.000000 
 
-     0.0 22+20      S2 O1 A0:    60.201360    77.812648    21.998067    20.005882 
+     0.0 22+20      S2 O1 A0:    60.201360    77.812648    22.000000    20.000000 
 
-     0.0 16-20      S2 O1 A0:    47.114909   213.442328    16.000674   -19.985612 
+     0.0 16-20      S2 O1 A0:    47.114909   213.442328    16.000000   -20.000000 
 
- 10000.0 22+20      S2 O1 A0:    66.960184   283.939460    22.002574    20.008999 
+ 10000.0 22+20      S2 O1 A0:    66.960184   283.939460    22.000000    20.000000 
 
  10000.0 16-20      S2 O1 A0:   163.034601   266.888368    16.000000   -20.000000 
 
- 10000.0 22+20      S2 O1 A0:    66.960184   283.939460    22.002574    20.008999 
+ 10000.0 22+20      S2 O1 A0:    66.960184   283.939460    22.000000    20.000000 
 
  10000.0 16-20      S2 O1 A0:   163.034601   266.888368    16.000000   -20.000000 
 
- 10000.0 22+20      S2 O1 A0:    70.239185   284.844639    22.003023    20.011289 
+ 10000.0 22+20      S2 O1 A0:    70.239185   284.844639    22.000000    20.000000 
 
  10000.0 16-20      S2 O1 A0:   166.422699   267.516844    16.000000   -20.000000 
 
 -10000.0 22+20      S2 O1 A1:   135.019352   331.536703    22.000000    20.000000 
 
 -10000.0 16-20      S2 O1 A1:   115.693174   103.225975    16.000000   -20.000000 
 
-     0.0 22+20      S2 O1 A1:    60.201360    77.812648    21.998067    20.005882 
+     0.0 22+20      S2 O1 A1:    60.201360    77.812648    22.000000    20.000000 
 
-     0.0 16-20      S2 O1 A1:    47.114909   213.442328    16.000674   -19.985612 
+     0.0 16-20      S2 O1 A1:    47.114909   213.442328    16.000000   -20.000000 
 
- 10000.0 22+20      S2 O1 A1:    66.960184   283.939460    22.002574    20.008999 
+ 10000.0 22+20      S2 O1 A1:    66.960184   283.939460    22.000000    20.000000 
 
  10000.0 16-20      S2 O1 A1:   163.034601   266.888368    16.000000   -20.000000 
 
- 10000.0 22+20      S2 O1 A1:    66.960184   283.939460    22.002574    20.008999 
+ 10000.0 22+20      S2 O1 A1:    66.960184   283.939460    22.000000    20.000000 
 
  10000.0 16-20      S2 O1 A1:   163.034601   266.888368    16.000000   -20.000000 
 
- 10000.0 22+20      S2 O1 A1:    70.239185   284.844639    22.003023    20.011289 
+ 10000.0 22+20      S2 O1 A1:    70.239185   284.844639    22.000000    20.000000 
 
  10000.0 16-20      S2 O1 A1:   166.422699   267.516844    16.000000   -20.000000 
diff --git a/test/src/test-super.c b/test/src/test-super.c
@@ -143,13 +143,13 @@ static int test_equ2hor() {
 
   if(obs.where != NOVAS_OBSERVER_ON_EARTH) return 0;
 
-  for(a = 0; a < 24.0; a += 3) {
+  for(a = 0; a < 24.0; a += 2) {
     int d;
     for(d = -90; d <= 90; d += 30) {
       double ra = a, dec = d, az, za, rar, decr;
 
-      if(!is_ok("itrs_to_hor:rar:null", equ2hor(tdb, 0.0, NOVAS_REDUCED_ACCURACY, 0.0, 0.0, &obs.on_surf, NOVAS_STANDARD_ATMOSPHERE, ra, dec, &az, &za, NULL, &decr))) return 1;
-      if(!is_ok("itrs_to_hor:decr:null", equ2hor(tdb, 0.0, NOVAS_REDUCED_ACCURACY, 0.0, 0.0, &obs.on_surf, NOVAS_STANDARD_ATMOSPHERE, ra, dec, &az, &za, &rar, NULL))) return 1;
+      if(!is_ok("equ2hor:rar:null", equ2hor(tdb, 0.0, NOVAS_REDUCED_ACCURACY, 0.0, 0.0, &obs.on_surf, NOVAS_STANDARD_ATMOSPHERE, ra, dec, &az, &za, NULL, &decr))) return 1;
+      if(!is_ok("equ2hor:decr:null", equ2hor(tdb, 0.0, NOVAS_REDUCED_ACCURACY, 0.0, 0.0, &obs.on_surf, NOVAS_STANDARD_ATMOSPHERE, ra, dec, &az, &za, &rar, NULL))) return 1;
     }
   }
 
@@ -228,13 +228,16 @@ static int test_tod_vs_cirs() {
 static int test_equ_ecl() {
   int a;
 
-  for(a = 0; a < 24; a += 3) {
-    double ra0 = a, dec0 = 0.0, elon, elat, ra, dec;
+  for(a = 0; a < 24; a += 2) {
+    int d;
+    for(d = -90; d <= 90; d += 30) {
+      double ra0 = a, dec0 = d, elon, elat, ra, dec;
 
-    if(!is_ok("equ2ecl", equ2ecl(tdb, NOVAS_GCRS_EQUATOR, NOVAS_FULL_ACCURACY, ra0, dec0, &elon, &elat))) return 1;
-    if(!is_ok("ecl2equ", ecl2equ(tdb, NOVAS_GCRS_EQUATOR, NOVAS_FULL_ACCURACY, elon, elat, &ra, &dec))) return 1;
-    if(!is_ok("equ_ecl_equ:ra", fabs(remainder(ra - ra0, 24.0)) > 1e-8)) return 1;
-    if(!is_ok("equ_ecl_equ:dec", fabs(dec - dec0) > 1e-7)) return 1;
+      if(!is_ok("equ2ecl", equ2ecl(tdb, NOVAS_GCRS_EQUATOR, NOVAS_FULL_ACCURACY, ra0, dec0, &elon, &elat))) return 1;
+      if(!is_ok("ecl2equ", ecl2equ(tdb, NOVAS_GCRS_EQUATOR, NOVAS_FULL_ACCURACY, elon, elat, &ra, &dec))) return 1;
+      if(!is_ok("equ_ecl_equ:ra", fabs(remainder((ra - ra0), 24.0) * cos(dec0 * DEG2RAD)) > 1e-8)) return 1;
+      if(!is_ok("equ_ecl_equ:dec", fabs(dec - dec0) > 1e-7)) return 1;
+    }
   }
 
   return 0;
@@ -244,13 +247,16 @@ static int test_equ_ecl() {
 static int test_equ_gal() {
   int a;
 
-  for(a = 0; a < 24; a += 3) {
-    double ra0 = a, dec0 = 0.0, glon, glat, ra, dec;
+  for(a = 0; a < 24; a += 2) {
+    int d;
+    for(d = -90; d <= 90; d += 30) {
+      double ra0 = a, dec0 = d, glon, glat, ra, dec;
 
-    if(!is_ok("equ2gal", equ2gal(ra0, dec0, &glon, &glat))) return 1;
-    if(!is_ok("gal2equ", gal2equ(glon, glat, &ra, &dec))) return 1;
-    if(!is_ok("equ_gal_equ:ra", fabs(remainder(ra - ra0, 24.0)) > 1e-8)) return 1;
-    if(!is_ok("equ_gal_equ:dec", fabs(dec - dec0) > 1e-7)) return 1;
+      if(!is_ok("equ2gal", equ2gal(ra0, dec0, &glon, &glat))) return 1;
+      if(!is_ok("gal2equ", gal2equ(glon, glat, &ra, &dec))) return 1;
+      if(!is_ok("equ_gal_equ:ra", fabs(remainder((ra - ra0), 24.0) * cos(dec0 * DEG2RAD)) > 1e-8)) return 1;
+      if(!is_ok("equ_gal_equ:dec", fabs(dec - dec0) > 1e-7)) return 1;
+    }
   }
 
   return 0;