More physical units in novas.h

README.md

@@ -349,6 +349,7 @@ switch between different planet and ephemeris calculator functions at will, duri
  - [Calculating positions for a Solar-system source](#solsys-example)
  - [Reduced accuracy shortcuts](#accuracy-notes)
  - [Performance considerations](#performance-note)
+ - [Physical units](#physical-units)
 
 
 <a name="methodologies"></a>
@@ -383,6 +384,7 @@ astrometric positions of celestial objects. The guide below is geared towards th
 NOVAS C approach remains viable also (albeit often less efficient). You may find an equivalent example usage 
 showcasing the original NOVAS method in [LEGACY.md](LEGACY.html).
 
+
 <a name="sidereal-example"></a>
 ### Calculating positions for a sidereal source
 
@@ -671,6 +673,62 @@ Just make sure that you:
    `config.mk` or in your equivalent build setup.
 
 
+<a name="physical-units"></a>
+### Physical units
+
+The NOVAS API has been using conventional units (e.g. AU, km, day, deg, h) typically for its parameters and return 
+values alike. Hence, SuperNOVAS follows the same conventions for its added functions and data structures also. 
+However, when interfacing SuperNOVAS with other programs, libraries, or data files, it is often necessary to use
+quantities that are expressed in different units, such as SI or CGS. To facilitate such conversions, `novas.h` 
+provides a set of unit constants, which can be used for converting to/from SI units (and radians). For example, 
+`novas.h` contains the following definitions:
+
+```c
+  /// [s] The length of a synodic day, that is 24 hours exactly. @since 1.2
+  #define NOVAS_DAY                 86400.0
+
+  /// [rad] A degree expressed in radians. @since 1.2
+  #define NOVAS_DEGREE              (M_PI / 180.0)
+
+  /// [rad] An hour of angle expressed in radians. @since 1.2
+  #define NOVAS_HOURANGLE           (M_PI / 12.0)
+```
+
+You can use these, for example, to convert quantities expressed in conventional units for NOVAS to standard (SI) 
+values, by multiplying NOVAS quantities with the corresponding unit definition. E.g.:
+
+```c
+  // A difference in Julian Dates [day] in seconds.
+  double delta_t = (tjd - tjd0) * NOVAS_DAY;
+
+  // R.A. [h] / declination [deg] converted radians (e.g. for trigonometric functions).
+  double ra_rad = ra_h * NOVAS_HOURANGLE;
+  double dec_rad = dec_d * NOVAS_DEGREE; 
+```
+
+And vice-versa: to convert values expressed in standard (SI) units, you can divide by the appropriate constant to
+'cast' an SI value into the particular physical unit, e.g.:
+
+```c
+  // Increment a Julian Date [day] with some time differential [s].
+  double tjd = tjd0 + delta_t / NOVAS_DAY;
+
+  // convert R.A. / declination in radians to hours and degrees
+  double ra_h = ra_rad / NOVAS_HOURANGLE;
+  double dec_d = dec_rad / NOVAS_DEGREE;
+```
+
+Finally, you can combine them to convert between two different conventional units, e.g.:
+
+```c
+  // Convert angle from [h] -> [rad] -> [deg]
+  double lst_d = lst_h * HOURANGLE / DEGREE; 
+
+  // Convert [AU/day] -> [m/s] (SI) -> [km/s]
+  double v_kms = v_auday * (NOVAS_AU / NOVAS_DAY) / NOVAS_KM
+```
+
+
 -----------------------------------------------------------------------------
 
 <a name="precision"></a>

include/novas.h

@@ -114,13 +114,30 @@
 /// [day] Julian date at B1900
 #define NOVAS_JD_B1900            15019.81352
 
-/// [day] Julian date for J1991.25, which the Hipparcos catalog is
-/// referred to
+/// [day] Julian date for J1991.25, which the Hipparcos catalog is referred to.
 #define NOVAS_JD_HIP              2448349.0625
 
 /// [m/s] Speed of light in meters/second is a defining physical constant.
 #define NOVAS_C                   299792458.0
 
+/// [s] The length of a synodic day, that is 24 hours exactly. @since 1.2
+#define NOVAS_DAY                 86400.0
+
+/// [rad] A degree expressed in radians. @since 1.2
+#define NOVAS_DEGREE              (M_PI / 180.0)
+
+/// [rad] An arc minute expressed in radians. @since 1.2
+#define NOVAS_ARCMIN              (NOVAS_DEGREE / 60.0)
+
+/// [rad] An arc second expressed in radians. @since 1.2
+#define NOVAS_ARCSEC              (NOVAS_ARCMIN / 60.0)
+
+/// [rad] An hour of angle expressed in radians. @since 1.2
+#define NOVAS_HOURANGLE           (M_PI / 12.0)
+
+/// [m] A kilometer (km) in meters. @since 1.2
+#define NOVAS_KM                  1000.0
+
 /// [m] Astronomical unit (AU). IAU definition.
 /// See <a href="https://www.iau.org/static/resolutions/IAU2012_English.pdf">IAU 2012 Resolution B2</a>.
 /// @sa DE405_AU
@@ -133,18 +150,16 @@
 /// [s] Light-time for one astronomical unit (AU) in seconds.
 #define NOVAS_AU_SEC              ( NOVAS_AU / NOVAS_C )
 
-/// [AU/day] Speed of light in AU/day.  Value is 86400 / AU_SEC.
-#define NOVAS_C_AU_PER_DAY        ( 86400.0 / AU_SEC )
+/// [AU/day] Speed of light in units of AU/day.
+#define NOVAS_C_AU_PER_DAY        ( NOVAS_DAY / AU_SEC )
 
-/// [km] Astronomical Unit in kilometers.
+/// [km] Astronomical Unit (AU) in kilometers.
 #define NOVAS_AU_KM               ( 1e-3 * NOVAS_AU )
 
-/// [m<sup>3</sup>/s<sup>2</sup>] Heliocentric gravitational constant in
-/// meters^3 / second^2, from DE-405.
+/// [m<sup>3</sup>/s<sup>2</sup>] Heliocentric gravitational constant, from DE-405.
 #define NOVAS_G_SUN               1.32712440017987e+20
 
-/// [m<sup>3</sup>/s<sup>2</sup>] Geocentric gravitational constant in
-/// meters^3 / second^2, from DE-405.
+/// [m<sup>3</sup>/s<sup>2</sup>] Geocentric gravitational constant, from DE-405.
 #define NOVAS_G_EARTH             3.98600433e+14
 
 /// [m] Solar radius (photosphere)
@@ -154,12 +169,10 @@
 /// [m] Radius of Earth in meters from IERS Conventions (2003).
 #define NOVAS_EARTH_RADIUS        6378136.6
 
-/// Earth ellipsoid flattening from IERS Conventions (2003). Value is
-/// 1 / 298.25642.
+/// Earth ellipsoid flattening from IERS Conventions (2003). Value is 1 / 298.25642.
 #define NOVAS_EARTH_FLATTENING    (1.0 / 298.25642)
 
-/// [rad/s] Rotational angular velocity of Earth in radians/sec from IERS
-/// Conventions (2003).
+/// [rad/s] Rotational angular velocity of Earth from IERS Conventions (2003).
 #define NOVAS_EARTH_ANGVEL        7.2921150e-5
 
 /// [s] TAI - GPS time offset
@@ -1671,14 +1684,14 @@ int mod_to_gcrs(double jd_tdb, const double *in, double *out);
 #define ANGVEL              NOVAS_EARTH_ANGVEL
 
 // Various locally used physical units
-#define DAY                 86400.0               ///< [s] seconds in a day
+#define DAY                 NOVAS_DAY
 #define DAY_HOURS           24.0
 #define DEG360              360.0
 #define JULIAN_YEAR_DAYS    365.25
 #define JULIAN_CENTURY_DAYS 36525.0
-#define ARCSEC              ASEC2RAD
-#define DEGREE              DEG2RAD
-#define HOURANGLE           (M_PI / 12.0)
+#define ARCSEC              NOVAS_ARCSEC
+#define DEGREE              NOVAS_DEGREE
+#define HOURANGLE           NOVAS_HOURANGLE
 #define MAS                 (1e-3 * ASEC2RAD)
 
 // On some older platform NAN may not be defined, so define it here if need be

src/novas.c

@@ -344,7 +344,7 @@ double novas_z2v(double z) {
   }
   z += 1.0;
   z *= z;
-  return 1e-3 * (z - 1.0) / (z + 1.0) * C;
+  return (z - 1.0) / (z + 1.0) * C / NOVAS_KM;
 }
 
 /**
@@ -2894,9 +2894,9 @@ int terra(const on_surface *location, double lst, double *pos, double *vel) {
   cosphi = cos(phi);
   c = 1.0 / sqrt(cosphi * cosphi + df2 * sinphi * sinphi);
   s = df2 * c;
-  ht_km = location->height / 1000.0;
-  ach = 1e-3 * ERAD * c + ht_km;
-  ash = 1e-3 * ERAD * s + ht_km;
+  ht_km = location->height / NOVAS_KM;
+  ach = ERAD * c / NOVAS_KM + ht_km;
+  ash = ERAD / NOVAS_KM * s + ht_km;
 
   // Compute local sidereal time factors at the observer's longitude.
   stlocl = lst * HOURANGLE + location->longitude * DEGREE;
@@ -4405,7 +4405,7 @@ double rad_vel2(const object *source, const double *pos_emit, const double *vel_
 
       // Compute radial velocity measure of sidereal source rel. barycenter
       // Including proper motion
-      beta_src = 1e3 * star->radialvelocity / C;
+      beta_src = NOVAS_KM * star->radialvelocity / C;
 
       if(star->parallax > 0.0) {
         double du[3];
@@ -5041,13 +5041,13 @@ int starvectors(const cat_entry *star, double *pos, double *vel) {
   if(vel) {
     // Compute Doppler factor, which accounts for change in
     // light travel time to star.
-    const double k = 1.0 / (1.0 - 1000.0 * star->radialvelocity / C);
+    const double k = 1.0 / (1.0 - NOVAS_KM * star->radialvelocity / C);
 
     // Convert proper motion and radial velocity to orthogonal components of
     // motion with units of AU/day.
     const double pmr = k * star->promora / (paralx * JULIAN_YEAR_DAYS);
     const double pmd = k * star->promodec / (paralx * JULIAN_YEAR_DAYS);
-    const double rvl = k * 1000.0 * star->radialvelocity / (AU / DAY);
+    const double rvl = k * NOVAS_KM * star->radialvelocity / (AU / DAY);
 
     // Transform motion vector to equatorial system.
     vel[0] = -pmr * sra - pmd * sdc * cra + rvl * cdc * cra;
@@ -6218,7 +6218,7 @@ short transform_cat(enum novas_transform_type option, double jd_tt_in, const cat
   pos[2] = dist * sdc;
 
   // Compute Doppler factor, which accounts for change in light travel time to star.
-  k = 1.0 / (1.0 - in->radialvelocity / C * 1000.0);
+  k = 1.0 / (1.0 - in->radialvelocity / C * NOVAS_KM);
 
   // Convert proper motion and radial velocity to orthogonal components
   // of motion, in spherical polar system at star's original position,

src/solsys-calceph.c

@@ -40,7 +40,7 @@
 #define CALCEPH_UNITS               (CALCEPH_UNIT_KM | CALCEPH_UNIT_DAY)
 
 /// Multiplicative normalization for the positions returned by CALCEPH to AU
-#define NORM_POS                    (1e3 / NOVAS_AU)
+#define NORM_POS                    (NOVAS_KM / NOVAS_AU)
 
 /// Multiplicative normalization for the velocities returned by CALCEPH to AU/day
 #define NORM_VEL                    (NORM_POS)

src/solsys-cspice.c

@@ -39,7 +39,7 @@
 #include "cspice/SpiceZpr.h"        // for reset_c
 
 /// Multiplicative normalization for the positions returned by km to AU
-#define NORM_POS                    (1e3 / NOVAS_AU)
+#define NORM_POS                    (NOVAS_KM / NOVAS_AU)
 
 /// Multiplicative normalization for the velocities returned by km/s to AU/day
 #define NORM_VEL                    (NORM_POS * 86400.0)

src/super.c

@@ -1182,9 +1182,9 @@ int make_solar_system_observer(const double *sc_pos, const double *sc_vel, obser
  * @since 1.2
  */
 double novas_v2z(double vel) {
-  vel *= 1e3 / C;   // [km/s] -> beta
+  vel *= NOVAS_KM / C;   // [km/s] -> beta
   if(fabs(vel) > 1.0) {
-    novas_error(-1, EINVAL, "novas_v2z", "velocity exceeds speed of light v=%g km/s", 1e-3 * vel * C);
+    novas_error(-1, EINVAL, "novas_v2z", "velocity exceeds speed of light v=%g km/s", vel * C / NOVAS_KM);
     return NAN;
   }
   return sqrt((1.0 + vel) / (1.0 - vel)) - 1.0;

test/src/test-super.c

@@ -26,8 +26,8 @@
 #define __NOVAS_INTERNAL_API__      ///< Use definitions meant for internal use by SuperNOVAS only
 #include "novas.h"
 
-#define J2000   2451545.0
-#define DAY     86400.0
+#define J2000   NOVAS_JD_J2000
+
 
 static char *workPath;