Tests for refract.c

Author: attipaci

include/novas.h

@@ -1144,6 +1144,8 @@ double novas_optical_refraction(double jd_tt, const on_surface *loc, enum novas_
 
 double novas_radio_refraction(double jd_tt, const on_surface *loc, enum novas_refraction_type type, double el);
 
+double novas_inv_refract(RefractionModel model, double jd_tt, const on_surface *loc, enum novas_refraction_type type, double el0);
+
 // <================= END of SuperNOVAS API =====================>
 
 #include "solarsystem.h"

src/novas.c

@@ -6521,8 +6521,8 @@ double refract(const on_surface *location, enum novas_refraction_model option, d
 
   zd_obs = fabs(zd_obs);
 
-  // Compute refraction up to zenith distance 90.1 degrees.
-  if(zd_obs > 90.1)
+  // Compute refraction up to zenith distance 91 degrees.
+  if(zd_obs > 91.0)
     return 0.0;
 
   // If observed weather data are available, use them.  Otherwise, use

src/refract.c

@@ -21,13 +21,18 @@
 
 
 static double novas_refraction(enum novas_refraction_model model, const on_surface *loc, enum novas_refraction_type type, double el) {
-  if(!loc)
-    return novas_error(-1, EINVAL, "novas_refraction", "NULL on surface observer location");
+  if(!loc) {
+    novas_error(-1, EINVAL, "novas_refraction", "NULL on surface observer location");
+    return NAN;
+  }
 
   if(type == NOVAS_REFRACT_OBSERVED)
-    return refract(loc, NOVAS_WEATHER_AT_LOCATION, 90.0 - el);
+    return refract(loc, model, 90.0 - el);
+
+  if(type == NOVAS_REFRACT_ASTROMETRIC)
+    return refract_astro(loc, model, 90.0 - el);
 
-  return refract_astro(loc, NOVAS_WEATHER_AT_LOCATION, 90.0 - el);
+  return NAN;
 }
 
 /**
@@ -120,6 +125,8 @@ double novas_optical_refraction(double jd_tt, const on_surface *loc, enum novas_
  * data is fully defined, and that the humidity was explicitly set after calling
  * `make_on_surface()`.
  *
+ * Adapted from FORTAN code provided by Berman & Rockwell 1976.
+ *
  * REFERENCES:
  * <ol>
  * <li>Berman, Allan L., and Rockwell, Stephen T. (1976), NASA JPL Technical Report 32-1601</li>
@@ -140,7 +147,6 @@ double novas_optical_refraction(double jd_tt, const on_surface *loc, enum novas_
 double novas_radio_refraction(double jd_tt, const on_surface *loc, enum novas_refraction_type type, double el) {
   // Various coefficients...
   static const double E[] = { 0.0, 46.625, 45.375, 4.1572, 1.4468, 0.25391, 2.2716, -1.3465, -4.3877, 3.1484, 4.520, -1.8982, 0.8900 };
-  //static const double EULER = 2.71828182845905; // Euler's number
 
   double E0, TK;
   double y = 1.0, z;
@@ -149,23 +155,27 @@ double novas_radio_refraction(double jd_tt, const on_surface *loc, enum novas_re
   double refraction;
   int j;
 
-  // The Bermann
-  if(type == NOVAS_REFRACT_OBSERVED) {
-    return novas_inv_refract(novas_radio_refraction, jd_tt, loc, NOVAS_REFRACT_ASTROMETRIC, el);
+  if(!loc) {
+    novas_error(-1, EINVAL, "novas_radio_refraction", "NULL on surface observer location");
+    return NAN;
   }
 
-  // TODO check which type of argument it takes.
+  if(type != NOVAS_REFRACT_OBSERVED && type != NOVAS_REFRACT_ASTROMETRIC) {
+    novas_error(-1, EINVAL, "novas_radio_refraction", "invalid refraction type: %d", type);
+    return NAN;
+  }
+
+  // The Bermann
+  if(type == NOVAS_REFRACT_OBSERVED)
+    return novas_inv_refract(novas_radio_refraction, jd_tt, loc, NOVAS_REFRACT_ASTROMETRIC, el);
 
   // Zenith angle in degrees
   z = (90.0 - el) / DEGREE;
 
   // Temperature in Kelvin
   TK = loc->temperature + 273.16;
-
   fptem = (loc->pressure / 1000.) * (273.16 / TK);
-
   E0 = (z - E[1]) / E[2];
-
   poly = E[11];
 
   for(j = 1; j <= 8; j++)
@@ -174,9 +184,7 @@ double novas_radio_refraction(double jd_tt, const on_surface *loc, enum novas_re
   if(poly <= -80.) poly = 0.;
 
   poly = exp(poly) - E[12];
-
   refraction = poly * fptem;
-
   y = exp(((TK * 17.149) - 4684.1) / (TK - 38.45));
 
   return refraction * 1.0 + (y * loc->humidity * 71.) / (TK * loc->pressure * 0.760);

src/timescale.c

@@ -74,6 +74,10 @@ int novas_set_time(enum novas_timescale timescale, double jd, int leap, double d
  * precision is reached if the fractional part is on the order of ≤= 1 day. In that case, the
  * time may be specified to picosecond accuracy, if needed.
  *
+ * The accuracy of Barycentric Time measures (TDB and TCB) relative to other time measures is
+ * limited by the precision of `tbd2tt()` implementation, to around 10 μs.
+ *
+ *
  * REFERENCES:
  * <ol>
  * <li>IAU 1991, RECOMMENDATION III. XXIst General Assembly of the
@@ -223,6 +227,9 @@ double novas_get_time(const novas_timespec *time, enum novas_timescale timescale
  * at the picosecond level, as opposed to `novas_set_time()`, whose precision is limited to a
  * ew microseconds. Ultimately, the accutacy of the time will depend on how it was set.
  *
+ * The accuracy of Barycentric Time measures (TDB and TCB) relative to other time measures is
+ * limited by the precision of the `tbd2tt()` implemenation, to around 10 &mu;s.
+ *
  * REFERENCES:
  * <ol>
  * <li>IAU 1991, RECOMMENDATION III. XXIst General Assembly of the

test/src/test-errors.c

@@ -162,7 +162,7 @@ static int test_refract() {
   if(check("refract:model", 1, r == 0.0 && errno == EINVAL)) n++;
 
   errno = 0;
-  r = refract(&o, NOVAS_STANDARD_ATMOSPHERE, 91.0);
+  r = refract(&o, NOVAS_STANDARD_ATMOSPHERE, 91.01);
   if(check("refract:zd", 1, r == 0.0)) n++;
 
   return n;
@@ -905,6 +905,28 @@ static int test_time() {
   return n;
 }
 
+static int test_refraction() {
+  int n = 0;
+  on_surface obs = {};
+
+  if(check_nan("stardard_refraction:loc", novas_standard_refraction(NOVAS_JD_J2000, NULL, NOVAS_REFRACT_OBSERVED, 10.0))) n++;
+  if(check_nan("stardard_refraction:type:-2", novas_standard_refraction(NOVAS_JD_J2000, &obs, -2, 10.0))) n++;
+  if(check_nan("stardard_refraction:type:1", novas_standard_refraction(NOVAS_JD_J2000, &obs, 1, 10.0))) n++;
+  if(check_nan("stardard_refraction:el:neg", novas_standard_refraction(NOVAS_JD_J2000, &obs, 1, -10.0))) n++;
+
+  if(check_nan("optical_refraction:loc", novas_optical_refraction(NOVAS_JD_J2000, NULL, NOVAS_REFRACT_OBSERVED, 10.0))) n++;
+  if(check_nan("optical_refraction:type:-2", novas_optical_refraction(NOVAS_JD_J2000, &obs, -2, 10.0))) n++;
+  if(check_nan("optical_refraction:type:1", novas_optical_refraction(NOVAS_JD_J2000, &obs, 1, 10.0))) n++;
+  if(check_nan("optical_refraction:el:neg", novas_optical_refraction(NOVAS_JD_J2000, &obs, 1, -10.0))) n++;
+
+  if(check_nan("radio_refraction:loc", novas_radio_refraction(NOVAS_JD_J2000, NULL, NOVAS_REFRACT_OBSERVED, 10.0))) n++;
+  if(check_nan("radio_refraction:type:-2", novas_radio_refraction(NOVAS_JD_J2000, &obs, -2, 10.0))) n++;
+  if(check_nan("radio_refraction:type:1", novas_radio_refraction(NOVAS_JD_J2000, &obs, 1, 10.0))) n++;
+  if(check_nan("radio_refraction:el:neg", novas_radio_refraction(NOVAS_JD_J2000, &obs, 1, -10.0))) n++;
+
+  return n;
+}
+
 int main() {
   int n = 0;
 
@@ -998,6 +1020,7 @@ int main() {
 
   if(test_obs_posvel()) n++;
   if(test_time()) n++;
+  if(test_refraction()) n++;
 
   if(n) fprintf(stderr, " -- FAILED %d tests\n", n);
   else fprintf(stderr, " -- OK\n");