diff --git a/README.md b/README.md index 8a907e0b..5f4d04aa 100644 --- a/README.md +++ b/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) @@ -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). + ### Calculating positions for a sidereal source @@ -671,6 +673,62 @@ Just make sure that you: `config.mk` or in your equivalent build setup. + +### 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 +``` + + ----------------------------------------------------------------------------- diff --git a/include/novas.h b/include/novas.h index abfb9f8e..6fcb0d19 100644 --- a/include/novas.h +++ b/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 IAU 2012 Resolution B2. /// @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 ) -/// [m3/s2] Heliocentric gravitational constant in -/// meters^3 / second^2, from DE-405. +/// [m3/s2] Heliocentric gravitational constant, from DE-405. #define NOVAS_G_SUN 1.32712440017987e+20 -/// [m3/s2] Geocentric gravitational constant in -/// meters^3 / second^2, from DE-405. +/// [m3/s2] 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 diff --git a/src/novas.c b/src/novas.c index f1e11713..ab783672 100644 --- a/src/novas.c +++ b/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, diff --git a/src/solsys-calceph.c b/src/solsys-calceph.c index 0f416bf6..899e22c1 100644 --- a/src/solsys-calceph.c +++ b/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) diff --git a/src/solsys-cspice.c b/src/solsys-cspice.c index 6320af3b..a74d432e 100644 --- a/src/solsys-cspice.c +++ b/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) diff --git a/src/super.c b/src/super.c index a73cee61..9386ad44 100644 --- a/src/super.c +++ b/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; diff --git a/test/src/test-super.c b/test/src/test-super.c index a5886f3a..8a4ca871 100644 --- a/test/src/test-super.c +++ b/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;