Fix a few tt2tdb conversions, expanded README

Smithsonian · Feb 6, 2024 · 7703532 · 7703532
1 parent abe2d58
commit 7703532
Show file tree

Hide file tree

Showing 3 changed files with 224 additions and 28 deletions.
diff --git a/README.md b/README.md
@@ -26,6 +26,7 @@ and easier to use. It is entirely free to use without any licensing restrictions
  - [Fixed NOVAS C 3.1 issues](#fixed-issues)
  - [Building and installation](#installation)
  - [Building your application with SuperNOVAS](#building-your-application)
+ - [Example usage](#examples)
  - [SuperNOVAS specific features](#supernovas-features)
  - [Notes on precision](#precision)
  - [External Solar-system ephemeris data or services](#solarsystem)
@@ -223,6 +224,195 @@ being unset in `config.mk`).
 
 -----------------------------------------------------------------------------
 
+
+<a name="examples"></a>
+## Example usage
+
+### Calculating positions for a sidereal source
+
+Sidereal sources may be anything beyond the solar-system with 'fixed' catalog coordinates. It may be a star, 
+or a galactic molecular cloud, or a distant quasar. First, you must provide the coordinates (which may include
+proper motion and parallax). Let's assume we pick a star for which we have B1950 (i.e. FK4) coordinates:
+
+```c
+  cat_entry source; // Structure to contain information on sidereal source 
+
+  // Let's assume we have B1950 (FK4) coordinates...
+  // 16h26m20.1918s, -26d19m23.138s (B1950), proper motion -12.11, -23.30 mas/year, parallax 5.89 mas, 
+  // radial velocity -3.4 km/s.
+  make_cat_entry("Antares", "FK4", 1, 16.43894213, -26.323094, -12.11, -23.30, 5.89, -3.4, &source);
+```
+
+We must convert these coordinates to the now standard ICRS system for calculations in NOVAS, first by calculating 
+equivalent J2000 coordinates, by applying the proper motion and the appropriate precession. Then, we apply a small 
+adjustment to convert from J2000 to ICRS coordinates.
+
+```c
+  // First change the catalog coordinates (in place) to the J2000... 
+  transform_cat(CHANGE_EPOCH, NOVAS_B1950, &source, NOVAS_J2000, &source);
+  
+  // Then convert J2000 coordinates to ICRS (also in place). Here the dates don't matter...
+  transform_cat(CHANGE_J2000_TO_ICRS, 0.0, &source, 0.0, &source);
+```
+
+Next, we define the location where we observe from. Here we can (but don't have to) specify local weather parameters
+(temperature and pressure) also for refraction correction later (in this example, we'll skip the weather):
+
+```c
+  observer obs;	 // Structure to contain observer location 
+
+  // Specify the location we are observing from
+  // 50.7374 deg N, 7.0982 deg E, 60m elevation
+  make_observer_on_surface(50.7374, 7.0982, 60.0, 0.0, 0.0, &obs);
+```
+
+We also need to set the time of observation. Clocks usually measure UTC, but for NOVAS we usually need time measured
+based on Terrestrial Time (TT) or Barycentric Time (TDB) or UT1. So 
+
+```c
+  // The current value for the leap seconds (UTC - TAI)
+  int leap_seconds = 37;
+
+  // Set the DUT1 = UT1 - UTC time difference in seconds (e.g. from IERS Bulletins)
+  int dut1 = ...
+
+  // Calculate the Terrestrial Time (TT) based Julian date of observation (in days)
+  // Let's say on 2024-02-06 at 14:53:06 UTC.
+  double jd_tt = julian_date(2024, 2, 6, tt_hour(14.885, leap_seconds)); 
+  
+  // We'll also need the TT - UT1 difference, which we can obtain from the above...
+  double ut1_to_tt = get_ut1_to_tt(jd_tt, dut1);
+```
+
+Now we can calculate the precise apparent position of the source, such as it's right ascension (R.A.) and declination, 
+and the equatorial _x,y,z_ unit vector pointing in the direction of the source (in the requested coordinate system 
+and for the specified observing location). We also get radial velocity (for spectroscopy), and distance (e.g. for 
+apparent-to-physical size conversion):
+
+```c
+  sky_pos pos;	// We'll return the topocantric apparent position in this structure
+
+  // We calculate the topocentric position (in True of Date [TOD] coordinates), for the above 
+  // configuration
+  int status = place(jd_tt, &source, &obs, ut1_to_tt, NOVAS_TOD, NOVAS_FULL_ACCURACY, &pos);
+
+  // You should always check that the calculation was successful...
+  if(status) {
+    // Ooops, something went wrong...
+    return -1;
+  }
+```
+
+Finally, we may want to calculate the azimuth and elevation at which the source is visible from the ground at the
+specified observing location (with or without refraction correction):
+
+```c
+  // Current polar offsets provided by the IERS Bulletins (in arcsec)
+  double dx = ... 
+  double dy = ...
+  
+  // Zenith distance and azimuth (in degrees) to populate...
+  double zd, az;
+  
+  // Convert to horizontal using the pole offsets, and a standard atmpsphere for
+  // the refraction. (We don't care to return refracted RA and dec coordinates).
+  equ2hor(jd_tt - ut1_to_tt, ut1_to_tt, NOVAS_FULL_ACCURACY, dx, dy, &obs.on_surf, 
+  	pos.ra, pos.dec, NOVAS_STANDARD_ATMOSPHERE, &zd, &az, NULL, NULL);
+``` 
+
+In the example above we have calculated the apparent coordinates in the Pseudo Earth Fixed (PEF) system, which does 
+not include polar wobble corrections, and applied the small (arcsec-level) measured variation of the pole (dx, dy) 
+only when calculating the precise actual Az/El coordinates. Alternatively, you may set the pole offsets before calling 
+`place()` using `cel_pole()`. In this case `place()` will calculate precise coordinates in the International 
+Terrestrial Reference System (ITRS), which we can convert to horizontal directly (without applying polar wobble 
+corrections twice):
+
+```c
+  cel_pole(jd_tt, POLE_OFFSETS_X_Y, dx, dy);
+  ...
+  place(...);
+  ...
+  equ2hor(... NOVAS_FULL_ACCURACY, 0.0, 0.0, &obs.on_surf ...);
+```
+
+#### Reduced accuracy
+
+When one needs positions to the arcsecond level only (and not to the microarcsecond level), some shortcuts can
+be made to the recipe above:
+
+ - You can use `NOVAS_REDUCED_ACCURACY` instead of `NOVAS_FULL_ACCURACY` for the calculations. This typically has an 
+   effect at the milliarcsecond level only, but may less computationally intense.
+ - You can skip the J2000 to ICRS conversion and use J2000 coordinates as a fair approximation (at the &lt;~ 22 mas
+   level).
+ - You can skip the pole offsets dx, dy. These may be a couple arcsec. 
+
+
+### Calculating positions for a Solar-system source
+
+Solar-system sources work similarly to the above with a few differences.
+
+First, You will have to provide one or more functions to obtain the barycentric ICRS positions for your Solar-system 
+source of interest for the specific Barycentric Dynamical Time of observation. See Section on integrating 
+[External Solar-system ephemeris data or services](#solarsystem) with SuperNOVAS.
+
+You can use `tt2tdb()` to help convert Terrstrial Time (TT) to Barycentric Dynamic Time (TDB) for your ephemeris 
+provider function (they only differ when you really need extreme precision -- for most applications you can used TT 
+and TDB interchangeably in the present era). E.g. 
+
+```c
+  double jd_tdb = jd_tt + tt2tdb(jd_tt) / 86400.0;
+```
+
+Instead of `make_cat_entry` you define your source as an `object` with an ID number that is used by the ephemeris 
+service you provided. For major planets you might want to use type `NOVAS_MAJOR_PLANET` if they use a 
+`novas_planet_calculator` function to access ephemeris data with their NOVAS IDs, or else `NOVAS_MINOR_PLANET` for 
+more generic ephemeris handling via a user-provided `novas_ephem_reader_func`. 
+
+```c
+  object mars, ceres; // Hold data on solar-system bodies.
+
+  // Mars will be handled by hte planet calculator function
+  make_object(NOVAS_MAJOR_PLANET_MARS, NOVAS_MARS, "Mars", NULL, &mars);
+
+  // Ceres will be handled by the generic ephemeris reader, which say uses the 
+  // NAIF ID of 2000001
+  make_object(NOVAS_MINOR_PLANET, 2000001, "Ceres", NULL, &ceres);
+```
+
+You will want to specify the functions that will handle the respective ephemeris data before making the NOVAS 
+calls that need them, e.g.:
+
+```c
+  // Set the function to use for regular precision planet position calculations
+  set_planet_calc(my_planet_calcular_function);
+  
+  // Set the function for high-precision planet position calculations
+  set_planet_calc(my_very_precise_planet_calculator_function);
+  
+  // Set the function to use for calculating all sorts of solar-system bodies
+  set_ephem_reader(my_ephemeris_provider_function);
+```
+
+Other than that, it's the same spiel as before. E.g.:
+
+```c
+  sky_pos pos;	// We'll return the topocentric apparent position in this structure
+
+  // We calculate the topocentric position (in True of Date [TOD] coordinates), for the above 
+  // configuration
+  int status = place(jd_tt, &mars, &obs, ut1_to_tt, NOVAS_TOD, NOVAS_FULL_ACCURACY, &pos);
+
+  // You should always check that the calculation was successful...
+  if(status) {
+    // Ooops, something went wrong...
+    return -1;
+  }
+```
+
+
+-----------------------------------------------------------------------------
+
+
 <a name="supernovas-features"></a>
 ## SuperNOVAS specific features
 

diff --git a/include/novas.h b/include/novas.h
@@ -310,15 +310,15 @@ enum novas_origin {
 enum novas_transform_type {
   /// Updates the star's data to account for the star's space motion between
   /// the first and second dates, within a fixed reference frame.
-  CHANGE_EPOCH = 1,
+  PROPER_MOTION = 1,
 
   /// applies a rotation of the reference frame
   /// corresponding to precession between the first and second dates,
   /// but leaves the star fixed in space.
-  CHANGE_EQUATOR_EQUINOX,
+  CHANGE_SYSTEM,
 
   /// The combined equivalent of CHANGE_EPOCH and CHANGE_EQUATOR_EQUINOX together.
-  CHANGE_EQUATOR_EQUINOX_EPOCH,
+  CHANGE_EPOCH,
 
   /// A fixed rotation about very small angles (<0.1 arcsecond) to take data from the
   /// dynamical system of J2000.0 to the ICRS.
@@ -732,6 +732,8 @@ int light_time2(double jd_tdb, const object *ss_object, const double *pos_obs, d
 
 double tt2tdb(double jd_tt);
 
+double tt_hour(double utchour, int leap_seconds);
+
 int cio_set_locator_file(const char *filename);
 
 int nutation_set_lp(novas_nutate_func f);

diff --git a/src/novas.c b/src/novas.c
@@ -203,7 +203,7 @@ static int j2000_to_tod(double jd_tt, enum novas_accuracy accuracy, const double
 
   if(!time_equals(jd_tt, JD_J2000)) {
     double v[3];
-    const double jd_tdb = jd_tt + tt2tdb(jd_tt);
+    const double jd_tdb = jd_tt + tt2tdb(jd_tt) / DAY;
     precession(JD_J2000, in, jd_tdb, v);
     nutation(jd_tdb, NUTATE_MEAN_TO_TRUE, accuracy, v, out);
 
@@ -240,7 +240,7 @@ static int tod_to_j2000(double jd_tt, enum novas_accuracy accuracy, const double
 
   if(!time_equals(jd_tt, JD_J2000)) {
     double v[3];
-    const double jd_tdb = jd_tt + tt2tdb(jd_tt);
+    const double jd_tdb = jd_tt + tt2tdb(jd_tt) / DAY;
     nutation(jd_tdb, NUTATE_TRUE_TO_MEAN, accuracy, in, v);
     precession(jd_tdb, v, JD_J2000, out);
   }
@@ -327,7 +327,7 @@ static int gcrs_to_cirs(double jd_tt, enum novas_accuracy accuracy, const double
   short rs;
   int error;
 
-  const double jd_tdb = jd_tt + tt2tdb(jd_tt);
+  const double jd_tdb = jd_tt + tt2tdb(jd_tt) / DAY;
 
   if(!in || !out) {
     errno = EINVAL;
@@ -376,7 +376,7 @@ static int cirs_to_gcrs(double jd_tt, enum novas_accuracy accuracy, const double
   short rs;
   int i;
 
-  const double jd_tdb = jd_tt + tt2tdb(jd_tt);
+  const double jd_tdb = jd_tt + tt2tdb(jd_tt) / DAY;
 
   if(!in || !out) {
     errno = EINVAL;
@@ -1820,7 +1820,7 @@ short ecl2equ_vec(double jd_tt, enum novas_equator_type coord_sys, enum novas_ac
  * @param xp          [arcsec] Conventionally-defined x coordinate of celestial intermediate pole with respect
  *                    to ITRS reference pole, in arcseconds.
  * @param yp          [arcsec] Conventionally-defined y coordinate of celestial intermediate pole with respect
- *                      to ITRS reference pole, in arcseconds.
+ *                    to ITRS reference pole, in arcseconds.
  * @param location    The observer location
  * @param ra          [h] Topocentric right ascension of object of interest, in hours, referred to true equator
  *                    and equinox of date.
@@ -5011,7 +5011,7 @@ int transform_hip(const cat_entry *hipparcos, cat_entry *hip_2000) {
  * are ignored.
  *
  * If 'option' is CHANGE_EPOCH, input data can be in any fixed reference
- * system. If 'option' is CHANGE_EQUATOR_EQUINOX or CHANGE_EQUATOR_EQUINOX_EPOCH, this function assumes
+ * system. If 'option' is CHANGE_SYSTEM or CHANGE_EPOCH, this function assumes
  * the input data is in the dynamical system and produces output in
  * the dynamical system. If 'option' is CHANGE_J2000_TO_ICRS, the input data must be
  * on the dynamical equator and equinox of J2000.0. And if 'option' is CHANGE_ICRS_TO_J2000,
@@ -5031,24 +5031,25 @@ int transform_hip(const cat_entry *hipparcos, cat_entry *hip_2000) {
  * @param date_newcat     [day|yr] Terrestrial Time (TT) based Julian date, or year, of output catalog data.
  * @param newcat_id       Catalog identifier (0 terminated)
  * @param[out] newcat     The transformed catalog entry, with units as given in the struct definition
- * @return                0 if successful, -1 if any of the pointer arguments is NULL or if the output
- *                        catalog is the same as the input, or else
- *                        1 if the input date is invalid for for option CHANGE_EQUATOR_EQUINOX or
- *                        CHANGE_EQUATOR_EQUINOX_EPOCH, or 2 if 'newcat_id' out of bounds.
+ * @return                0 if successful, -1 if any of the pointer arguments is NULL, or else
+ *                        1 if the input date is invalid for for option CHANGE_SYSTEM or
+ *                        CHANGE_EPOCH, or 2 if 'newcat_id' out of bounds.
  */
 short transform_cat(enum novas_transform_type option, double date_incat, const cat_entry *incat, double date_newcat, const char *newcat_id,
         cat_entry *newcat) {
 
   double jd_incat, jd_newcat, paralx, dist, r, d, cra, sra, cdc, sdc, k;
   double pos1[3], term1, pmr, pmd, rvl, vel1[3], pos2[3], vel2[3], xyproj;
 
-  if(!incat || !newcat || !newcat_id || newcat == incat) {
+  if(!incat || !newcat || !newcat_id) {
     errno = EINVAL;
     return -1;
   }
 
   if(strlen(newcat_id) >= sizeof(newcat->starname)) return 2;
 
+  // TODO move frame tie up front...
+
   // If necessary, compute Julian dates.
 
   // This function uses TDB Julian dates internally, but no distinction between TDB and TT is necessary.
@@ -5097,7 +5098,7 @@ short transform_cat(enum novas_transform_type option, double date_incat, const c
   memcpy(vel2, vel1, sizeof(vel2));
 
   // Update star's position vector for space motion (only if 'option' = 1 or 'option' = 3).
-  if((option == CHANGE_EPOCH) || (option == CHANGE_EQUATOR_EQUINOX_EPOCH)) {
+  if((option == PROPER_MOTION) || (option == CHANGE_EPOCH)) {
     int j;
     for(j = 0; j < 3; j++)
       pos2[j] = pos1[j] + vel1[j] * (jd_newcat - jd_incat);
@@ -5106,8 +5107,8 @@ short transform_cat(enum novas_transform_type option, double date_incat, const c
 
   switch(option) {
 
-    case CHANGE_EQUATOR_EQUINOX:
-    case CHANGE_EQUATOR_EQUINOX_EPOCH: {
+    case CHANGE_SYSTEM:
+    case CHANGE_EPOCH: {
       // Precess position and velocity vectors (only if 'option' = 2 or 'option' = 3).
       int error;
 
@@ -5153,7 +5154,6 @@ short transform_cat(enum novas_transform_type option, double date_incat, const c
   dist = vlen(pos2);
 
   paralx = asin(1.0 / dist) / MAS;
-  newcat->parallax = paralx;
 
   // Transform motion vector back to spherical polar system at star's new position.
   cra = cos(r);
@@ -5169,18 +5169,22 @@ short transform_cat(enum novas_transform_type option, double date_incat, const c
   newcat->promodec = pmd * paralx * JULIAN_YEAR_DAYS / k;
   newcat->radialvelocity = rvl * (AU_KM / DAY) / k;
 
-  // Take care of zero-parallax case.
-  if(newcat->parallax <= 1.01e-6) {
-    newcat->parallax = 0.0;
-    newcat->radialvelocity = incat->radialvelocity;
-  }
 
-  // Set the catalog identification code for the transformed catalog entry.
-  strncpy(newcat->catalog, newcat_id, sizeof(newcat->catalog) - 1);
+  if(newcat != incat) {
+    // Take care of zero-parallax case.
+    if(incat->parallax <= 0.0) {
+      newcat->parallax = 0.0;
+      newcat->radialvelocity = incat->radialvelocity;
+    }
+    else newcat->parallax = incat->parallax;
+
+    // Set the catalog identification code for the transformed catalog entry.
+    strncpy(newcat->catalog, newcat_id, sizeof(newcat->catalog) - 1);
 
-  // Copy unchanged quantities from the input catalog entry to the transformed catalog entry.
-  strncpy(newcat->starname, incat->starname, sizeof(newcat->starname) - 1);
-  newcat->starnumber = incat->starnumber;
+    // Copy unchanged quantities from the input catalog entry to the transformed catalog entry.
+    strncpy(newcat->starname, incat->starname, sizeof(newcat->starname) - 1);
+    newcat->starnumber = incat->starnumber;
+  }
 
   return 0;
 }