From 04b40859e54e563ea6c96efa389f77c190ee1261 Mon Sep 17 00:00:00 2001 From: Attila Kovacs <attipaci@gmail.com> Date: Sun, 19 Jan 2025 22:29:31 +0100 Subject: [PATCH] Added functions for v1.3 --- CHANGELOG.md | 57 ++++ Makefile | 2 + README.md | 5 +- include/novas.h | 90 ++++++- include/solarsystem.h | 10 + src/frames.c | 585 +++++++++++++++++++++++++++++++++++++++++ src/novas.c | 67 ++++- src/super.c | 436 +++++++++++++++++++++++++++++- test/src/test-errors.c | 216 +++++++++++++++ test/src/test-super.c | 379 ++++++++++++++++++++++++++ 10 files changed, 1838 insertions(+), 9 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index c0cefc60..2f6e31c5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,6 +7,63 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), [Semantic Versioning](https://semver.org/spec/v2.0.0.html). +## [Unreleased] + +Upcoming feature release, expected around 1 May 2025. + +### Added + + - New `novas_hms_hours(const char *str)` and `novas_dms_degrees(const char *str)` convenience functions to make it + easier to parse HMS or DMS based time/angle values, returning the result in units of hours or degrees, + appropriately for use in SuperNOVAS. + + - New `novas_frame_lst()` convenience function to readily return the Local (apparent) Sidereal Time for a given + Earth-based observing frame. + + - New `novas_rises_above()` and `novas_sets_below()` functions to return the date/time a source rises above or sets + below a specific elevation on a given date. (Useful for Earth-based observers only). + + - New `novas_helio_dist()` function to calculate the heliocentric distance of a Solar-system body on a given date. + The `novas_solar_power()` function can be used to estimate the incident Solar power on a Solar-system body, while + `novas_solar_illum()` can be used to calculate the fraction of a spherical body that is illuminated by the Sun seen + from the observer location. + + - New `novas_hpa()` and `novas_epa()` functions to calculate the parallactic angle (a.k.a. vertical position angle) + for a given location on sky, using the local horizontal coordinates, or else the equatorial position, respectively. + The parallactic angle (PA) can be useful to convert local Cartesian offsets (e.g. from a flat image or detector + array) between the local horizontal and equatorial orientations, e.g. via the newly added `novas_h2e_offset()` or + `novas_e2j_offset()` functions. The conversion between offsets and absolute coordinates usually requires a WCS + projections, such as described in Calabretta & Greisen 2002. + + - New `novas_sep()`, `novas_equ_sep()`, and `novas_object_sep()` functions can be used to calculate the precise + apparent distance between to spherical or equatorial locations, or between two sources, respectively. + `novas_sun_angle()` and `novas_moon_angle()` can be used to calculate the apparent angular distance of sources from + the Sun and Moon, respectively. + + - New `novas_observable` and `novas_track` data structures to provide second order Taylor series expansion of the + apparent horizontal or equatorial positions, distances and redshifts for sources. They can be calculated with the + newly added `novas_hor_track()` or `novas_equ_track()` functions. Such tracking values, including rates and + accelerations can be directly useful for controlling telescope drives in horizontal or equatorial mounts to track + sources (hence the name). You can also obtain instantaneous projected (extrapolated) positions from the tracking + parameters via `novas_track_pos()` at low computational cost. + + - New `novas_xyz2uvw()` function to convert ITRS Earth locations (absolute or differential) to equatorial projections + along a line of sight in the direction of a source. Such projections are oft used in interferometry. + + - #114: New `novas_lsr2ssb_vel()` can be used to convert velocity vectors referenced to the LSR to Solar-System + Barycentric velocities. (The new function is used by `place()`, and `novas_geom_posvel()` internally.) + +### Changed + + - #114: Velocities returned by `novas_geom_posvel()` and radial velocity measures returned in `sky_pos` (e.g. by + `novas_sky_pos()` or `place()`) were not observer-based velocities for catalog sources. Instead, they were + referenced to the LSR frame, with corrections for observer motion about the SSB, but not for the motion of the SSB + in the LSR frame. As such, the radial velocities for catalog sources did not provide the expected spectroscopic + relation between observed and rest wavelengths i.e., _f_<sub>obs</sub> = (1 + _rv_ / _c_) _f_<sub>rest</sub>. This + is a NOVAS C issue, which affected prior SuperNOVAS releases also. + + + ## [1.2.0] - 2025-01-15 Feature release. New easy to use adapter modules for CALCEPH or the NAIF CSPICE Toolkit to provide precise positions diff --git a/Makefile b/Makefile index a4b2dde7..fa6d2478 100644 --- a/Makefile +++ b/Makefile @@ -88,6 +88,8 @@ solsys: $(SOLSYS_TARGETS) all: distro static test coverage analyze # Run regression tests +export + .PHONY: test test: $(MAKE) -C test run diff --git a/README.md b/README.md index 2da99801..ee34a96d 100644 --- a/README.md +++ b/README.md @@ -33,7 +33,7 @@ SuperNOVAS is entirely free to use without licensing restrictions. Its source c standard, and hence should be suitable for old and new platforms alike. It is light-weight and easy to use, with full support for the IAU 2000/2006 standards for sub-microarcsecond position calculations. -This document has been updated for the `v1.2` and later releases. +This document has been updated for the `v1.3` and later releases. ## Table of Contents @@ -150,6 +150,9 @@ provided by SuperNOVAS over the upstream NOVAS C 3.1 code: - [__v1.1__] The NOVAS C 3.1 implementation of `rad_vel()` has a number of issues that produce inaccurate results. The errors are typically at or below the tens of m/s level for objects not moving at relativistic speeds. + - [__v1.3__] The radial velocity measure returned in `sky_pos.rv` for catalog sources was not expressed relative the + observer. As such, it did not provide an observable spectroscopic measure for catalog sources. This issue affected + NOVAS C 3.1 and prior SuperNOVAS releases alike. ----------------------------------------------------------------------------- diff --git a/include/novas.h b/include/novas.h index ecd23157..1b95169e 100644 --- a/include/novas.h +++ b/include/novas.h @@ -62,10 +62,10 @@ #define SUPERNOVAS_MAJOR_VERSION 1 /// API minor version -#define SUPERNOVAS_MINOR_VERSION 2 +#define SUPERNOVAS_MINOR_VERSION 3 /// Integer sub version of the release -#define SUPERNOVAS_PATCHLEVEL 1 +#define SUPERNOVAS_PATCHLEVEL 0 /// Additional release information in version, e.g. "-1", or "-rc1", or empty string "" for releases. #define SUPERNOVAS_RELEASE_STRING "-devel" @@ -181,6 +181,12 @@ /// [s] TT - TAI time offset #define NOVAS_TAI_TO_TT 32.184 +/// [W/m<sup>2</sup>] The Solar Constant i.e., typical incident Solar power on Earth. +/// The value of 1367 Wm<sup>−2</sup> was adopted by the World Radiation Center +/// (Gueymard, 2004). +/// @since 1.3 +#define NOVAS_SOLAR_CONSTANT 1367.0 + #if !COMPAT // If we are not in the strict compatibility mode, where constants are defined @@ -1015,13 +1021,16 @@ typedef struct { * * @sa place() * @sa SKY_POS_INIT + * @sa novas_z_lsr() */ typedef struct { double r_hat[3]; ///< unit vector toward object (dimensionless) double ra; ///< [h] apparent, topocentric, or astrometric right ascension (hours) double dec; ///< [deg] apparent, topocentric, or astrometric declination (degrees) double dis; ///< [AU] true (geometric, Euclidian) distance to solar system body or 0.0 for star (AU) - double rv; ///< [km/s] radial velocity (km/s) + double rv; ///< [km/s] radial velocity (km/s). As of SuperNOVAS v1.3, this is always a proper + ///< observer-based spectroscopic velocity measure, which relates the observed wavelength + ///< to the rest wavelength as λ<sub>obs</sub> = (1 + rv / c) λ<sub>rest</sub>. } sky_pos; /** @@ -1156,6 +1165,7 @@ typedef struct { novas_matrix nutation; ///< nutation matrix (Lieske 1977 method) novas_matrix gcrs_to_cirs; ///< GCRS to CIRS conversion matrix novas_planet_bundle planets; ///< Planet positions and velocities (ICRS) + // TODO [v2] add ra_cio } novas_frame; /** @@ -1178,6 +1188,7 @@ typedef struct { novas_matrix matrix; ///< Transformation matrix elements } novas_transform; + /** * The type of elevation value for which to calculate a refraction. * @@ -1259,6 +1270,35 @@ extern int grav_bodies_full_accuracy; */ typedef double (*RefractionModel)(double jd_tt, const on_surface *loc, enum novas_refraction_type type, double el); +/** + * Spherical and spectral coordinate set. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_track + */ +typedef struct { + double lon; ///< [deg] apparent longitude coordinate in coordinate system + double lat; ///< [deg] apparent latitude coordinate in coordinate system + double dist; ///< [AU] apparent distance to source from observer + double z; ///< redshift +} novas_observable; + +/** + * The spherical and spectral tracking position of a source, and its first and second time derivatives. + * + * @since 1.3 + * @author Attila Kovacs + */ +typedef struct { + novas_timespec time; ///< The astronomical time for which the track is calculated. + novas_observable pos; ///< [deg,AU,1] Apparent source position + novas_observable rate; ///< [deg/s,AU/s,1/s] Apparent position rate of change + novas_observable accel; ///< [deg/s<sup>2</sup>,AU/s<sup>2</sup>,1/s<sup>2</sup>] Apparent position acceleration. +} novas_track; + + short app_star(double jd_tt, const cat_entry *star, enum novas_accuracy accuracy, double *ra, double *dec); @@ -1651,8 +1691,48 @@ int gcrs_to_mod(double jd_tdb, const double *in, double *out); int mod_to_gcrs(double jd_tdb, const double *in, double *out); -// <================= END of SuperNOVAS API =====================> +// ---------------------- Added in 1.3.0 ------------------------- +int novas_lsr2ssb_vel(const double *vLSR, double unit, double *vSSB); + +double novas_hms_hours(const char *hms); + +double novas_dms_degrees(const char *dms); + +double novas_hpa(double az, double el, double lat); + +double novas_epa(double ha, double dec, double lat); + +int novas_h2e_offset(double daz, double del, double pa, double *dra, double *ddec); + +int novas_e2h_offset(double dra, double ddec, double pa, double *daz, double *del); + +double novas_sep(double lon1, double lat1, double lon2, double lat2); + +double novas_equ_sep(double ra1, double dec1, double ra2, double dec2); + +int novas_xyz2uvw(const double *xyz, double ha, double dec, double *uvw); + +double novas_frame_lst(const novas_frame *frame); + +double novas_rises_above(double el, const object *source, const novas_frame *frame, RefractionModel ref_model); + +double novas_sets_below(double el, const object *source, const novas_frame *frame, RefractionModel ref_model); + +double novas_object_sep(const object *source1, const object *source2, const novas_frame *frame); + +double novas_sun_angle(const object *source, const novas_frame *frame); + +double novas_moon_angle(const object *source, const novas_frame *frame); + +int novas_equ_track(const object *source, const novas_frame *frame, double dt, novas_track *track); + +int novas_hor_track(const object *source, const novas_frame *frame, RefractionModel ref_model, novas_track *track); + +int novas_track_pos(const novas_track *track, const novas_timespec *time, double *lon, double *lat, double *dist, double *z); + + +// <================= END of SuperNOVAS API =====================> #include <solarsystem.h> @@ -1736,6 +1816,8 @@ int novas_error(int ret, int en, const char *from, const char *desc, ...); return __ret; \ } +double novas_add_beta(double beta1, double beta2); + double novas_vlen(const double *v); double novas_vdist(const double *v1, const double *v2); double novas_vdot(const double *v1, const double *v2); diff --git a/include/solarsystem.h b/include/solarsystem.h index a51a7a9e..eaf88598 100644 --- a/include/solarsystem.h +++ b/include/solarsystem.h @@ -361,6 +361,16 @@ long novas_to_naif_planet(enum novas_planet id); long novas_to_dexxx_planet(enum novas_planet id); + +// Added in v1.3 ---------------------------------> + +double novas_helio_dist(double jd_tdb, const object *source, double *rate); + +double novas_solar_power(double jd_tdb, const object *source); + +double novas_solar_illum(const object *source, const novas_frame *frame); + + /// \cond PRIVATE diff --git a/src/frames.c b/src/frames.c index bdd62424..ca64a599 100644 --- a/src/frames.c +++ b/src/frames.c @@ -34,6 +34,8 @@ #define FRAME_INITIALIZED 0xdeadbeadcafeba5e ///< frame.state for a properly initialized frame. #define GEOM_TO_APP 1 ///< Geometric to apparent conversion #define APP_TO_GEOM (-1) ///< Apparent to geometric conversion + +#define NOVAS_TRACK_DELTA 30.0 ///< [s] Time step for evaluation horizontal tracking derivatives. /// \endcond static int cmp_sys(enum novas_reference_system a, enum novas_reference_system b) { @@ -470,6 +472,9 @@ static int icrs_to_sys(const novas_frame *frame, double *pos, enum novas_referen * 2006) method is used, with the Lieske et al. 1977 nutation model, matching the behavior of the * original NOVAS C place() for that system. To obtain more precise TOD coordinates, set `sys` to * `NOVAS_CIRS` here, and follow with cirs_to_tod() after.</li> + * <li>As of SuperNOVAS v1.3, the returned velocity vector is a proper observer-based + * velocity measure. In prior releases, and in NOVAS C 3.1, this was inconsistent, with + * pseudo LSR-based measures being returned for catalog sources.</li> * </ol> * * @param source Pointer to a celestial source data structure that is observed @@ -528,6 +533,9 @@ int novas_geom_posvel(const object *source, const novas_frame *frame, enum novas // Get position of star wrt observer (corrected for parallax). bary2obs(pos1, frame->obs_pos, pos1, &t_light); + + // Change velocity reference from LSR to SSB. + novas_lsr2ssb_vel(vel1, NOVAS_AU / NOVAS_DAY, vel1); } else { int got = 0; @@ -590,6 +598,9 @@ int novas_geom_posvel(const object *source, const novas_frame *frame, enum novas * `NOVAS_CIRS` here, and follow with cirs_to_tod() / cirs_to_app_ra() on the `out->r_hat` / * `out->ra` respectively after (or you can use just convert one of the quantities, and use * radec2vector() or vector2radec() to get the other even faster).</li> + * <li>As of SuperNOVAS v1.3, the returned radial velocity component is a proper observer-based + * spectroscopic measure. In prior releases, and in NOVAS C 3.1, this was inconsistent, with + * LSR-based measures being returned for catalog sources.</li> * </ol> * * @param object Pointer to a celestial object data structure that is observed @@ -1212,3 +1223,577 @@ int novas_transform_sky_pos(const sky_pos *in, const novas_transform *transform, return 0; } + +/** + * Returns the Local (apparent) Sidereal Time for an observing frame of an Earth-bound observer. + * + * @param frame Observer frame, defining the location and time of observation + * @return [h] The LST for an Earth-bound observer [0.0--24.0), or NAN otherwise. If NAN is + * returned errno will indicate the type of error. + * + * @since 1.3 + * @author Attila Kovacs + */ +double novas_frame_lst(const novas_frame *frame) { + static const char *fn = "novas_frame_lst"; + double lst; + + if(!frame) + return novas_error(-1, EINVAL, fn, "input frame is NULL"); + + if(!is_frame_initialized(frame)) + return novas_error(-1, EINVAL, fn, "input frame is not initialized"); + + if(frame->observer.where != NOVAS_OBSERVER_ON_EARTH && frame->observer.where != NOVAS_AIRBORNE_OBSERVER) + return novas_error(-1, EINVAL, fn, "Not an Earth-bound observer: where=%d", frame->observer.where); + + lst = remainder(frame->gst + frame->observer.on_surf.longitude / 15.0, 24.0); + if(lst < 0.0) lst += 24.0; + + return lst; +} + +/** + * Returns the hourangle at which the object crosses the specified elevation angle. + * + * + * @param el [rad] Elevation angle. + * @param dec [rad] Apparent Source declination. + * @param lat [rad] Geodetic latitude of observer. + * @return [h] the hour angle at which the source crosses the specified elevation + * angle, or else NAN if the source stays above or below the specified + * elevation at all times. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_sets_below() + */ +static double calc_lha(double el, double dec, double lat) { + double c = (sin(el) - sin(lat) * sin(dec)) / (cos(lat) * cos(dec)); + return acos(c) / NOVAS_HOURANGLE; +} + +/** + * Returns the UTC date at which a source appears cross the specified elevation angle. The calculated + * time will account for the motion of the source (for Solar-system objects), and optionally for atmospheric + * refraction also. + * + * @param el [deg] Elevation angle. + * @param sign 1 for rise time, or -1 for setting time. + * @param source Observed source + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @param ref_model Refraction model, or NULL to calculate unrefracted rise time. + * @return [day] UTC-based Julian date at which the object crosses the specified elevation + * in the 24 hour period after the specified date, or else NAN if the source stays + * above or below the given elevation for the entire 24-hour period. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_sets_below() + */ +static double novas_cross_el_date(double el, int sign, const object *source, const novas_frame *frame, RefractionModel ref_model) { + static const char *fn = "novas_cross_el_time"; + + const on_surface *loc; + novas_frame frame1; + sky_pos pos = {}; + double lst, utc2tt, lastRA = NAN; + int i; + + if(!source) { + novas_error(0, EINVAL, fn, "input source is NULL"); + return NAN; + } + + if(!frame) { + novas_error(0, EINVAL, fn, "input frame is NULL"); + return NAN; + } + + if(!is_frame_initialized(frame)) { + novas_error(0, EINVAL, fn, "input frame is not initialized"); + return NAN; + } + + lst = novas_frame_lst(frame); + if(isnan(lst)) + return novas_trace_nan(fn); + + if(ref_model) { + // Apply refraction correction + double ref = ref_model(novas_get_time(&frame->time, NOVAS_TT), &frame->observer.on_surf, NOVAS_REFRACT_OBSERVED, el); + if(isnan(ref)) + return novas_trace_nan(fn); + if(ref > 0.0) el -= ref / 3600.0; + } + + el *= DEGREE; // convert to degrees. + frame1 = *frame; // Time shifted frame + loc = (on_surface *) &frame->observer.on_surf; // Earth-bound location + utc2tt = (frame->time.dut1 + frame->time.ut1_to_tt) / DAY; + + for(i = 0; i < novas_inv_max_iter; i++) { + double lha, tUTC; + + prop_error(fn, novas_sky_pos(source, &frame1, NOVAS_TOD, &pos), 0); + + // Hourangle when source crosses nominal elevation + lha = calc_lha(el, pos.dec * NOVAS_DEGREE, loc->latitude * NOVAS_DEGREE); + if(isnan(lha)) + return novas_trace_nan(fn); + + // Calculate transit UTC at observer location + tUTC = remainder(pos.ra - novas_frame_lst(frame), 24.0); + if(tUTC < 0.0) tUTC += 24.0; + + // Adjusted frame time for last crossing time estimate + frame1.time.ijd_tt = frame->time.ijd_tt; + frame1.time.fjd_tt = utc2tt + (tUTC + sign * lha) / 24.0; + + if(frame1.time.fjd_tt < frame->time.fjd_tt) frame1.time.ijd_tt++; // Make sure to check rise/set time after input frame time. + + if(source->type == NOVAS_CATALOG_OBJECT) + break; // That's it for catalog sources + if(fabs(remainder(pos.ra - lastRA, 24.0)) < 1e-8) + break; // Check if converged (ms precision) + + lastRA = pos.ra; + } + + if(i >= novas_inv_max_iter) { + novas_error(0, ECANCELED, fn, "failed to converge"); + return NAN; + } + + return novas_get_time(&frame1.time, NOVAS_UTC); +} + +/** + * Returns the UTC date at which a source appears to rise above the specified elevation angle. The + * calculated time will account for the motion of the source (for Solar-system objects), and optionally + * for atmospheric refraction also. + * + * + * @param el [deg] Elevation angle. + * @param source Observed source + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @param ref_model Refraction model, or NULL to calculate unrefracted rise time. + * @return [day] UTC-based Julian date at which the object rises above the specified elevation + * in the 24 hour period after the specified date, or else NAN if the source stays + * above or below the given elevation for the entire 24-hour period. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_sets_below() + */ +double novas_rises_above(double el, const object *source, const novas_frame *frame, RefractionModel ref_model) { + double utc = novas_cross_el_date(el, -1, source, frame, ref_model); + if(isnan(utc)) + return novas_trace_nan("novas_rises_above"); + return utc; +} + +/** + * Returns the UTC date at which a source appears to set below the specified elevation angle. The + * calculated time will account for the motion of the source (for Solar-system objects), and optionally + * for atmopsheric refraction also. + * + * + * @param el [deg] Elevation angle. + * @param source Observed source + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @param ref_model Refraction model, or NULL to calculate unrefracted setting time. + * @return [day] UTC-based Julian date at which the object sets below the specified elevation + * in the 24 hour period after the specified date, or else NAN if the source stays + * above or below the given elevation for the entire 24-hour day.. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_rises_above() + */ +double novas_sets_below(double el, const object *source, const novas_frame *frame, RefractionModel ref_model) { + double utc = novas_cross_el_date(el, 1, source, frame, ref_model); + if(isnan(utc)) + return novas_trace_nan("novas_sets_below"); + return utc; +} + +/** + * Returns the Solar illumination fraction of a source, assuming a spherical geometry for the observed body. + * + * @param source Observed source. Usually a Solar-system source. (For other source types, 1.0 + * is returned by default.) + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @return Solar illumination fraction [0.0:1.0] of a spherical body observed at the + * source location from the given observer location, or NAN if there was an error + * (errno will indicate the type of error). + * + * @since 1.3 + * @author Attila Kovacs + */ +double novas_solar_illum(const object *source, const novas_frame *frame) { + static const char *fn = "novas_solar_illum"; + + double pos[3], dSrc, dObs, dSun; + int i; + + if(!source) { + novas_error(0, EINVAL, fn, "input source is NULL"); + return NAN; + } + + if(!frame) { + novas_error(0, EINVAL, fn, "input frame is NULL"); + return NAN; + } + + if(!is_frame_initialized(frame)) { + novas_error(0, EINVAL, fn, "input frame is not initialized"); + return NAN; + } + + if(source->type == NOVAS_CATALOG_OBJECT) + return 1.0; + + if(novas_geom_posvel(source, frame, NOVAS_ICRS, pos, NULL) != 0) + return novas_trace_nan(fn); + + dSrc = novas_vlen(pos); + + for(i = 3; --i >= 0; ) pos[i] += frame->obs_pos[i]; + + dSun = novas_vdist(pos, frame->sun_pos); + dObs = novas_vdist(frame->obs_pos, frame->sun_pos); + + return 0.5 + 0.5 * (dSrc * dSrc + dSun * dSun - dObs * dObs) / (2.0 * dSun * dSrc); +} + +/** + * Returns the angular separation of two objects from the observer's point of view. + * The calculated separation includes light-time corrections, aberration and gravitational + * deflection for both sources, and thus represents a precise observed separation between + * the two sources. + * + * @param source1 An observed source + * @param source2 Another observed source + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @return [deg] Apparent angular separation between the two observed sources + * from the observer's point-of-view. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_sun_angle() + * @sa novas_moon_angle() + * @sa novas_sep() + */ +double novas_object_sep(const object *source1, const object *source2, const novas_frame *frame) { + static const char *fn = "novas_object_sep"; + + sky_pos p1 = {}, p2 = {}; + + if(!source1) { + novas_error(0, EINVAL, fn, "input source1 is NULL"); + return NAN; + } + + if(!source2) { + novas_error(0, EINVAL, fn, "input source2 is NULL"); + return NAN; + } + + if(novas_sky_pos(source1, frame, NOVAS_GCRS, &p1) != 0) + return novas_trace_nan(fn); + + if(novas_sky_pos(source2, frame, NOVAS_GCRS, &p2) != 0) + return novas_trace_nan(fn); + + if(p1.dis < 1e-11 || p2.dis < 1e-11) { + novas_error(0, EINVAL, fn, "source is at observer location"); + return NAN; + } + + return novas_equ_sep(p1.ra, p1.dec, p2.ra, p2.dec); +} + +/** + * Returns the apparent angular distance of a source from the Sun from the observer's + * point of view. + * + * @param source An observed source + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @return [deg] the apparent angular distance between the source an the Sun, from + * the observer's point of view + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_moon_angle() + */ +double novas_sun_angle(const object *source, const novas_frame *frame) { + object sun = NOVAS_SUN_INIT; + double d = novas_object_sep(source, &sun, frame); + if(isnan(d)) + return novas_trace_nan("novas_sun_angle"); + return d; +} + +/** + * Returns the apparent angular distance of a source from the Moon from the observer's + * point of view. + * + * @param source An observed source + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @return [deg] Apparent angular distance between the source an the Moon, from + * the observer's point of view + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_sun_angle() + */ +double novas_moon_angle(const object *source, const novas_frame *frame) { + object moon = NOVAS_MOON_INIT; + double d = novas_object_sep(source, &moon, frame); + if(isnan(d)) + return novas_trace_nan("novas_moon_angle"); + return d; +} + +/** + * Calculates equatorial tracking position and motion (first and second time derivatives) for the + * specified source in the given observing frame. The position and its derivatives are calculated + * via the more precise IAU2006 method, and CIRS. + * + * @param source Observed source + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @param dt [s] Time step used for calculating derivatives. + * @param[out] track Output tracking parameters to populate + * @return 0 if successful, or else -1 if any of the pointer arguments are NULL, + * or else an error code from cio_ra() or from novas_sky_pos(). + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_hor_track() + * @sa novas_track_pos() + */ +int novas_equ_track(const object *source, const novas_frame *frame, double dt, novas_track *track) { + static const char *fn = "novas_equ_track"; + + novas_timespec time1; + novas_frame frame1; + sky_pos pos0 = {}, posm = {}, posp = {}; + double ra_cio; + double idt2; + + if(dt <= 0.0) dt = NOVAS_TRACK_DELTA; + idt2 = 1.0 / (dt * dt); + + if(!source) + return novas_error(-1, EINVAL, fn, "input source is NULL"); + + if(!frame) + return novas_error(-1, EINVAL, fn, "input frame is NULL"); + + if(!is_frame_initialized(frame)) + return novas_error(-1, EINVAL, fn, "input frame is not initialized"); + + if(!track) + return novas_error(-1, EINVAL, fn, "output track is NULL"); + + track->time = frame->time; + prop_error(fn, cio_ra(frame->time.ijd_tt + frame->time.fjd_tt, frame->accuracy, &ra_cio), 0); + + prop_error(fn, novas_sky_pos(source, frame, NOVAS_CIRS, &pos0), 0); + pos0.ra += ra_cio; + pos0.rv = novas_v2z(pos0.rv); + + track->pos.lon = 15.0 * pos0.ra; + track->pos.lat = pos0.dec; + track->pos.dist = pos0.dis; + track->pos.z = pos0.rv; + + time1 = frame->time; + time1.fjd_tt -= dt / DAY; + prop_error(fn, novas_make_frame(frame->accuracy, &frame->observer, &time1, frame->dx, frame->dy, &frame1), 0); + prop_error(fn, novas_sky_pos(source, &frame1, NOVAS_CIRS, &posm), 0); + posm.ra += ra_cio; + posm.rv = novas_v2z(posm.rv); + + time1.fjd_tt += 2.0 * dt / DAY; + prop_error(fn, novas_make_frame(frame->accuracy, &frame->observer, &time1, frame->dx, frame->dy, &frame1), 0); + prop_error(fn, novas_sky_pos(source, &frame1, NOVAS_CIRS, &posp), 0); + posp.ra += ra_cio; + posp.rv = novas_v2z(posp.rv); + + // Careful with RA wraps. + if(posm.ra > 18.0 || posp.ra > 18.0 || pos0.ra > 18.0) { + if(posm.ra < 6.0) posm.ra += 24.0; + if(posp.ra < 6.0) posp.ra += 24.0; + if(pos0.ra < 6.0) pos0.ra += 24.0; + } + + track->rate.lon = 7.5 * (posp.ra - posm.ra) / dt; + track->rate.lat = 0.5 * (posp.dec - posm.dec) / dt; + track->rate.dist = 0.5 * (posp.dis - posm.dis) / dt; + track->rate.z = 0.5 * (posp.rv - posm.rv) / dt; + + track->accel.lon = 15.0 * (0.5 * (posp.ra + posm.ra) - pos0.ra) * idt2; + track->accel.lat = (0.5 * (posp.dec + posm.dec) - pos0.dec) * idt2; + track->accel.dist = (0.5 * (posp.dis + posm.dis) - pos0.dis) * idt2; + track->accel.z = (0.5 * (posp.rv + posm.rv) - pos0.rv) * idt2; + + return 0; +} + +/** + * Calculates horizontal tracking position and motion (first and second time derivatives) for the + * specified source in the given observing frame. The position and its derivatives are calculated + * via the more precise IAU2006 method, and CIRS, and then converted to local horizontal + * coordinates using the specified refraction model (if any). + * + * @param source Observed source + * @param frame Observing frame, defining the observer location and astronomical time + * of observation. + * @param ref_model Refraction model to use, or NULL for an unrefracted track. + * @param[out] track Output tracking parameters to populate + * @return 0 if successful, or else -1 if any of the pointer arguments are NULL, + * or else an error code from cio_ra() or from novas_sky_pos(), or from + * novas_app_hor(). + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_equ_track() + * @sa novas_track_pos() + */ +int novas_hor_track(const object *source, const novas_frame *frame, RefractionModel ref_model, novas_track *track) { + static const char *fn = "novas_equ_track"; + + novas_timespec time1; + novas_frame frame1; + sky_pos pos = {}; + double ra_cio; + double az0, el0, azm, elm, dm, zm, azp, elp, dp, zp; + const double idt2 = 1.0 / (NOVAS_TRACK_DELTA * NOVAS_TRACK_DELTA); + + if(!source) + return novas_error(-1, EINVAL, fn, "input source is NULL"); + + if(!frame) + return novas_error(-1, EINVAL, fn, "input frame is NULL"); + + if(!is_frame_initialized(frame)) + return novas_error(-1, EINVAL, fn, "input frame is not initialized"); + + if(frame->observer.where != NOVAS_OBSERVER_ON_EARTH && frame->observer.where != NOVAS_AIRBORNE_OBSERVER) + return novas_error(-1, EINVAL, fn, "observer is not Earth-bound: where = %d", frame->observer.where); + + if(!track) + return novas_error(-1, EINVAL, fn, "output track is NULL"); + + track->time = frame->time; + prop_error(fn, cio_ra(frame->time.ijd_tt + frame->time.fjd_tt, frame->accuracy, &ra_cio), 0); + + prop_error(fn, novas_sky_pos(source, frame, NOVAS_CIRS, &pos), 0); + prop_error(fn, novas_app_to_hor(frame, NOVAS_TOD, pos.ra + ra_cio, pos.dec, ref_model, &az0, &el0), 0); + track->pos.lon = az0; + track->pos.lat = el0; + track->pos.dist = pos.dis; + track->pos.z = novas_v2z(pos.rv); + + time1 = frame->time; + time1.fjd_tt -= NOVAS_TRACK_DELTA / DAY; + prop_error(fn, novas_make_frame(frame->accuracy, &frame->observer, &time1, frame->dx, frame->dy, &frame1), 0); + prop_error(fn, novas_sky_pos(source, &frame1, NOVAS_CIRS, &pos), 0); + prop_error(fn, novas_app_to_hor(&frame1, NOVAS_TOD, pos.ra + ra_cio, pos.dec, ref_model, &azm, &elm), 0); + dp = pos.dis; + zm = novas_v2z(pos.rv); + + time1.fjd_tt += 2.0 * NOVAS_TRACK_DELTA / DAY; + prop_error(fn, novas_make_frame(frame->accuracy, &frame->observer, &time1, frame->dx, frame->dy, &frame1), 0); + prop_error(fn, novas_sky_pos(source, &frame1, NOVAS_CIRS, &pos), 0); + prop_error(fn, novas_app_to_hor(&frame1, NOVAS_TOD, pos.ra + ra_cio, pos.dec, ref_model, &azp, &elp), 0); + dm = pos.dis; + zp = novas_v2z(pos.rv); + + // Careful with Az wraps + if(azm > 270.0 || azp > 270.0 || az0 > 270.0) { + if(azm < 90.0) azm += 360.0; + if(azp < 90.0) azp += 360.0; + if(az0 < 90.0) az0 += 360.0; + } + + track->rate.lon = 0.5 * (azp - azm) / NOVAS_TRACK_DELTA; + track->rate.lat = 0.5 * (elp - elm) / NOVAS_TRACK_DELTA; + track->rate.dist = 0.5 * (dp - dm) / NOVAS_TRACK_DELTA; + track->rate.z = 0.5 * (zp - zm) / NOVAS_TRACK_DELTA; + + track->accel.lon = (0.5 * (azp + azm) - az0) * idt2; + track->accel.lat = (0.5 * (elp + elm) - el0) * idt2; + track->accel.dist = (0.5 * (dp + dm) - track->pos.dist) * idt2; + track->accel.z = (0.5 * (zp + zm) - track->pos.z) * idt2; + + return 0; +} + +/** + * Calculates a projected position and redshift for a source, given the available tracking position and + * derivatives. Using 'tracks' to project positions can be much faster than the repeated full recalculation + * of the source position over some short period. + * + * In SuperNOVAS terminology a 'track' is a 2nd order Taylor series expansion of the observed position and + * redshift in time. For most but the fastest moving sources, horizontal (Az/El) tracks are sufficiently + * precise on minute timescales, whereas depending on the type of source equatorial tracks can be precise for + * up to days. + * + * @param track Tracking position and motion (first and second derivatives) + * @param time Astrometric time of observation + * @param[out] lon [deg] projected observed Eastward longitude in tracking coordinate system + * @param[out] lat [deg] projected observed latitude in tracking coordinate system + * @param[out] dist [AU] projected apparent distance to source from observer + * @param[out] z projected observed redshift + * @return 0 if successful, or else -1 if either input pointer is NULL + * (errno is set to EINVAL). + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_equ_track() + * @sa novas_hor_track() + * @sa novas_z2v() + */ +int novas_track_pos(const novas_track *track, const novas_timespec *time, double *lon, double *lat, double *dist, double *z) { + static const char *fn = "novas_track_pos"; + + double dt, dt2; + + if(!time) + return novas_error(-1, EINVAL, fn, "input time is NULL"); + + if(!track) + return novas_error(-1, EINVAL, fn, "input track is NULL"); + + dt = novas_diff_time(time, &track->time); + dt2 = dt * dt; + + if(lon) *lon = remainder(track->pos.lon + track->rate.lon * dt + track->accel.lon * dt2, 360.0); + if(lat) *lat = track->pos.lat + track->rate.lat * dt + track->accel.lat * dt2; + if(dist) *dist = track->pos.dist + track->rate.dist * dt + track->accel.dist * dt2; + if(z) *z = track->pos.z + track->rate.z * dt + track->accel.z * dt2; + + return 0; +} diff --git a/src/novas.c b/src/novas.c index 8d5a4574..3a82a6ce 100644 --- a/src/novas.c +++ b/src/novas.c @@ -301,7 +301,7 @@ double novas_vdot(const double *v1, const double *v2) { return (v1[0] * v2[0]) + (v1[1] * v2[1]) + (v1[2] * v2[2]); } -static double novas_add_beta(double beta1, double beta2) { +double novas_add_beta(double beta1, double beta2) { return (beta1 + beta2) / (1 + beta1 * beta2); } @@ -766,6 +766,13 @@ novas_planet_provider_hp get_planet_provider_hp() { * given its catalog mean place, proper motion, parallax, and radial velocity. See `place()` * for more information. * + * NOTES: + * <ol> + * <li>As of SuperNOVAS v1.3, the returned radial velocity component is a proper observer-based + * spectroscopic measure. In prior releases, and in NOVAS C 3.1, this was inconsistent, with + * pseudo LSR-based measures being returned for catalog sources.</li> + * </ol> + * * REFERENCES: * <ol> * <li>Kaplan, G. H. et. al. (1989). Astron. Journ. 97, 1197-1210.</li> @@ -1547,6 +1554,54 @@ int obs_posvel(double jd_tdb, double ut1_to_tt, enum novas_accuracy accuracy, co return 0; } +/** + * Converts a Local Standard of Rest (LSR) based velocity vector to a Solar-System Barycentric velocity + * vector. It is used internally to convert LSR-based velocities ro spectroscopic observable radial + * velocity measures calculated for `sky_pos` data structures, in `novas_geom_posvel()` and `place()` alike. + * + * The SSB motion w.r.t. the barycenter is assumed to be (11.1, 12.24, 7.25) km/s in ICRS (Shoenrich et al. + * 2010). + * + * REFERENCES: + * <ol> + * <li>Ralph Schoenrich, James Binney, Walter Dehnen, Monthly Notices of the Royal Astronomical Society, + * Volume 403, Issue 4, April 2010, Pages 1829–1833, https://doi.org/10.1111/j.1365-2966.2010.16253.x</li> + * </ol> + * + * @param[in] vLSR [see unit below] Velocity 3-vector relative to the Local Standard of Rest (LSR). + * @param unit Physical unit in which velocity is expressed as a multiplicative factor to use to + * get S.I. quantities. E.g. (NOVAS_AU / NOVAS_DAY) if velocities are in AU/day. + * @param[out] vSSB [see unit above] Velocity 3-vector relative to the Solar System Barycenter (SSB). + * + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa rad_vel2() + * @sa novas_geom_posvel() + * @sa place() + * @sa sky_pos + */ +int novas_lsr2ssb_vel(const double *vLSR, double unit, double *vSSB) { + static const char *fn = "novas_lsr2ssb_vel"; + + static const double beta0[3] = { 11.1 * NOVAS_KM / C, 12.24 * NOVAS_KM / C, 7.25 * NOVAS_KM / C }; + int i; + + if(!vLSR) + return novas_error(-1, EINVAL, fn, "input vLSR is NULL"); + + if(!vSSB) + return novas_error(-1, EINVAL, fn, "output vSSB is NULL"); + + if(unit <= 0.0) + return novas_error(-1, EINVAL, fn, "invalid velocity unit conversion: %g", unit); + + for(i = 3; -- i >= 0; ) vSSB[i] = novas_add_beta(vLSR[i] * unit / C, -beta0[i]) * C / unit; + + return 0; +} + /** * Computes the apparent direction of a celestial object at a specified time and in a specified * coordinate system and a specific near-Earth origin. @@ -1577,6 +1632,9 @@ int obs_posvel(double jd_tdb, double ut1_to_tt, enum novas_accuracy accuracy, co * 2006) method is used, with the Lieske et al. 1977 nutation model, matching the behavior of the * original NOVAS C place() for that system. To obtain more precise TOD coordinates, set `sys` to * `NOVAS_CIRS` here, and follow with cirs_to_tod() after.</li> + * <li>As of SuperNOVAS v1.3, the returned radial velocity component is a proper observer-based + * spectroscopic measure. In prior releases, and in NOVAS C 3.1, this was inconsistent, with + * pseudo LSR-based measures being returned for catalog sources.</li> * </ol> * * REFERENCES: @@ -1701,6 +1759,9 @@ short place(double jd_tt, const object *source, const observer *location, double output->dis = 0.0; d_sb = novas_vlen(pos); + + // Change velocity reference from LSR to SSB. + novas_lsr2ssb_vel(vel, NOVAS_AU / NOVAS_DAY, vel); } else { // Get position of body wrt observer, antedated for light-time. @@ -4446,7 +4507,7 @@ double rad_vel2(const object *source, const double *pos_emit, const double *vel_ du[1] = uk[1] - (cosdec * sin(ra)); du[2] = uk[2] - sin(dec); - beta_src += novas_vdot(vel_src, du) / C_AUDAY; + beta_src = novas_add_beta(beta_src, novas_vdot(vel_src, du) / C_AUDAY); } break; @@ -6793,7 +6854,7 @@ short make_object(enum novas_object_type type, long number, const char *name, co if(!source) return novas_error(-1, EINVAL, fn, "NULL input source"); - // FIXME will not need special case in v2.x + // FIXME [v2] will not need special case in v2.x memset(source, 0, type == NOVAS_ORBITAL_OBJECT ? sizeof(object) : offsetof(object, orbit)); // Set the object type. diff --git a/src/super.c b/src/super.c index d7531f69..a15fd9bf 100644 --- a/src/super.c +++ b/src/super.c @@ -1428,7 +1428,7 @@ enum novas_planet novas_planet_for_name(const char *name) { * @param type Coordinate reference system in which `ra` and `dec` are defined (e.g. NOVAS_GCRS). * @param ra [h] the R.A. of the pole of the oribtal reference plane. * @param dec [deg] the declination of the pole of the oribtal reference plane. - * @param[out] sys The orbital system + * @param[out] sys Orbital system * @return 0 if successful, or else -1 (errno will be set to EINVAL) if the output `sys` * pointer is NULL. * @@ -1448,3 +1448,437 @@ int novas_set_orbsys_pole(enum novas_reference_system type, double ra, double de return 0; } + +/** + * Returns the decimal hours for a HMS string specification. The hour, minute, and second + * components may be separated by spaces, tabs, colons `:`, or a combination thereof. + * Additionally, the hour and minutes may be semarated by the letter `h` or `H`, and the + * minutes and seconds may be separated by `m` or `M`, or a single quote `'`. For example, all + * of the lines below specify the same time: + * + * <pre> + * 23:59:59.999 + * 23h 59m 59.999 + * 23h 59' 59.999 + * 23H59'59.999 + * </pre> + * + * + * @param hms String specifying hours, minutes, and seconds, which correspond to + * a time between 0 and 24 h. Time in any range is permitted, but the minutes and + * seconds must be >=0 and <60. + * @return [hours] Corresponding decimal time value, or else NAN if there was an + * error parsing the string (errno will be set to EINVAL). + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_dms_degrees() + */ +double novas_hms_hours(const char *hms) { + static const char *fn = "novas_hms_hours"; + + int h = 0, m = 0; + double s = NAN; + + if(!hms) { + novas_error(-1, EINVAL, fn, "input string is NULL"); + return NAN; + } + if(!hms[0]) { + novas_error(-1, EINVAL, fn, "input string is empty"); + return NAN; + } + + if(sscanf(hms, "%d%*[:hH \t]%d%*[:mM' \t]%lf", &h, &m, &s) != 3) { + novas_error(-1, EINVAL, fn, "not in HMS format: '%s'", hms); + return NAN; + } + + if(m < 0 || m >= 60) { + novas_error(-1, EINVAL, fn, "invalid minutes: got %d, expected 0-59", m); + return NAN; + } + + if(s < 0.0 || s >= 60.0) { + novas_error(-1, EINVAL, fn, "invalid seconds: got %f, expected [0.0:60.0)", s); + return NAN; + } + + s = abs(h) + (m / 60.0) + (s / 3600.0); + return h < 0 ? -s : s; +} + +/** + * Returns the decimal degrees for a DMS string specification. The degree, (arc)minute, and + * (arc)second components may be separated by spaces, tabs, colons `:`, or a combination thereof. + * Additionally, the degree and minutes may be semarated by the letter `d` or `D`, and the + * minutes and seconds may be separated by `m` or `M`, or a single quote `'`. For example, all of + * the lines below specify the same angle: + * + * <pre> + * -179:59:59.999 + * -179 59m 59.999 + * -179d 59' 59.999 + * -179D59'59.999 + * </pre> + * + * + * @param dms String specifying degrees, minutes, and seconds, which correspond to + * an angle. Angles in any range are permitted, but the minutes and + * seconds must be >=0 and <60. + * @return [deg] Corresponding decimal angle value, or else NAN if there was + * an error parsing the string (errno will be set to EINVAL). + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_hms_hours() + */ +double novas_dms_degrees(const char *dms) { + static const char *fn = "novas_dms_degrees"; + + int d = 0, m = 0; + double s = NAN; + + if(!dms) { + novas_error(-1, EINVAL, fn, "input string is NULL"); + return NAN; + } + if(!dms[0]) { + novas_error(-1, EINVAL, fn, "input string is empty"); + return NAN; + } + + if(sscanf(dms, "%d%*[:dD \t]%d%*[:mM' \t]%lf", &d, &m, &s) != 3) { + novas_error(-1, EINVAL, fn, "not in DMS format: '%s'", dms); + return NAN; + } + + if(m < 0 || m >= 60) { + novas_error(-1, EINVAL, fn, "invalid minutes: got %d, expected 0-59", m); + return NAN; + } + + if(s < 0.0 || s >= 60.0) { + novas_error(-1, EINVAL, fn, "invalid seconds: got %f, expected [0.0:60.0)", s); + return NAN; + } + + s = abs(d) + (m / 60.0) + (s / 3600.0); + + return d < 0 ? -s : s; +} + +/** + * Returns the horizontal Parallactic Angle (PA) calculated for a gorizontal Az/El location of the sky. The PA + * is the angle between the local horizontal coordinate directions and the local true-of-date equatorial + * coordinate directions at the given location. The polar wobble is not included in the calculation. + * + * The Parallactic Angle is sometimes referrred to as the Vertical Position Angle (VPA). Both define the + * same quantity. + * + * @param az [deg] Azimuth angle + * @param el [deg] Elevation angle + * @param lat [deg] Geodetic latitude of observer + * @return [deg] Parallactic Angle (PA). I.e., the clockwise position angle of the declination direction + * w.r.t. the elevation axis in the horizontal system. Same as the the clockwise position angle + * of the elevation direction w.r.t. the declination axis in the equatorial system. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_epa() + * @sa novas_h2e_offset() + */ +double novas_hpa(double az, double el, double lat) { + double s, c; + + lat *= DEGREE; + az *= DEGREE; + el *= DEGREE; + + s = sin(lat); + c = cos(lat); + + return atan2(-c * sin(az), s * cos(el) - c * sin(el) * cos(az)) / DEGREE; +} + +/** + * Returns the equatorial Parallactic Angle (PA) calculated for an R.A./Dec location of the sky at a given + * sidereal time. The PA is the angle between the local horizontal coordinate directions and the local + * true-of-date equatorial coordinate directions, at the given location and time. The polar wobble is not + * included in the calculation. + * + * The Parallactic Angle is sometimes referrred to as the Vertical Position Angle (VPA). Both define the + * same quantity. + * + * @param ha [h] Hour angle (LST - RA) i.e., the difference between the Local (apparent) Sidereal Time + * and the apparent (true-of-date) Right Ascension of observed source. + * @param dec [deg] Apparent (true-of-date) declination of observed source + * @param lat [deg] Geodetic latitude of observer + * @return [deg] Parallactic Angle (PA). I.e., the clockwise position angle of the elevation direction + * w.r.t. the declination axis in the equatorial system. Same as the clockwise position angle + * of the declination direction w.r.t. the elevation axis, in the horizontal system. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_hpa() + * @sa novas_lst() + * @sa novas_e2h_offset() + */ +double novas_epa(double ha, double dec, double lat) { + double coslat; + + ha *= HOURANGLE; + lat *= DEGREE; + dec *= DEGREE; + + coslat = cos(lat); + return atan2(coslat * sin(ha), sin(lat) * cos(dec) - coslat * sin(dec) * cos(ha)) / DEGREE; +} + +/** + * Converts coordinate offsets, from the local horizontal system to local equatorial offsets. + * Converting between local flat projections and spherical coordinates usually requires a WCS + * projection. + * + * REFERENCES: + * <ol> + * <li>Calabretta, M.R., & Greisen, E.W., (2002), Astronomy & Astrophysics, 395, 1077-1122.</li> + * </ol> + * + * @param daz [arcsec] Projected offset position in the azimuth direction. The projected + * offset between two azimuth positions at the same reference elevation is + * δAz = (Az2 - Az1) * cos(El<sub>0</sub>). + * @param del [arcsec] projected offset position in the elevation direction + * @param pa [deg] Parallactic Angle + * @param[out] dra [arcsec] Output offset position in the local true-of-date R.A. direction. It + * can be a pointer to one of the input coordinates, or NULL if not required. + * @param[out] ddec [arcsec] Output offset position in the local true-of-date declination + * direction. It can be a pointer to one of the input coordinates, or NULL if not + * required. + * @return 0 + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_e2h_offset() + * @sa novas_hpa() + */ +int novas_h2e_offset(double daz, double del, double pa, double *dra, double *ddec) { + double dx = daz, dy = del, c, s; + + pa *= DEGREE; + c = cos(pa); + s = sin(pa); + + if(dra) *dra = s * dy - c * dx; + if(ddec) *ddec = s * dx + c * dy; + + return 0; +} + +/** + * Converts coordinate offsets, from the local equatorial system to local horizontal offsets. + * Converting between local flat projections and spherical coordinates usually requires a WCS + * projection. + * + * REFERENCES: + * <ol> + * <li>Calabretta, M.R., & Greisen, E.W., (2002), Astronomy & Astrophysics, 395, 1077-1122.</li> + * </ol> + * + * @param dra [arcsec] Projected ffset position in the apparent true-of-date R.A. direction. + * E.g. The projected offset between two RA coordinates at a same reference + * declination, is + * δRA = (RA2 - RA1) * cos(Dec<sub>0</sub>) + * @param ddec [arcsec] Projected offset position in the apparent true-of-date declination + * direction. + * @param pa [deg] Parallactic Angle + * @param[out] daz [arcsec] Output offset position in the local azimuth direction. It can be a pointer + * to one of the input coordinates, or NULL if not required. + * @param[out] del [arcsec] Output offset position in the local elevation direction. It can be a + * pointer to one of the input coordinates, or NULL if not required. + * @return 0 + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_h2e_offset() + * @sa novas_epa() + */ +int novas_e2h_offset(double dra, double ddec, double pa, double *daz, double *del) { + return novas_h2e_offset(dra, ddec, pa, daz, del); +} + +/** + * Returns a Solar-system body's distance from the Sun, and optionally also the rate of recession. + * It may be useful, e.g. to calculate the body's heating from the Sun. + * + * @param jd_tdb [day] Barycentric Dynamical Time (TDB) based Julian date. You may want to + * use a time that is antedated to when the observed light originated from the + * source. + * @param source Observed Solar-system source + * @param[out] rate [AU/day] (optional) Returned rate of recession from Sun + * @return [AU] Distance from the Sun, or NAN if not a Solar-system source. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_solar_power() + * @sa novas_solar_illum() + */ +double novas_helio_dist(double jd_tdb, const object *source, double *rate) { + static const char *fn = "novas_helio_dist"; + + const double jd2[2] = { jd_tdb, 0.0 }; + double pos[3], vel[3], d; + + if(!source) { + novas_error(0, EINVAL, fn, "input source is NULL"); + if(rate) *rate = NAN; + return NAN; + } + + if(source->type == NOVAS_CATALOG_OBJECT) { + novas_error(0, EINVAL, fn, "input source is not a Solar-system body: type %d", source->type); + if(rate) *rate = NAN; + return NAN; + } + + if(ephemeris(jd2, source, NOVAS_HELIOCENTER, NOVAS_REDUCED_ACCURACY, pos, vel) != 0) return novas_trace_nan(fn); + + d = novas_vlen(pos); + if(!d) { + // The Sun itself... + if(rate) *rate = 0.0; + return 0.0; + } + + if(rate) *rate = novas_vlen(vel); + return d; + +} + +/** + * Returns the incident Solar power on a Solar-system body at the time of observation. + * + * @param jd_tdb [day] Barycentric Dynamical Time (TDB) based Julian date. You may want to + * use a time that is antedated to when the observed light originated from the + * source. + * @param source Observed Solar-system source + * @return [W/m<sup>2</sup>] Incident Solar power on the illuminated side of the object, + * or NAN if not a Solar-system source or if the source is the Sun itself. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_solar_illum() + */ +double novas_solar_power(double jd_tdb, const object *source) { + double d = novas_helio_dist(jd_tdb, source, NULL); + return NOVAS_SOLAR_CONSTANT / (d * d); +} + +/** + * Returns the angular separation of two locations on a sphere. + * + * @param lon1 [deg] longitude of first location + * @param lat1 [deg] latitude of first location + * @param lon2 [deg] longitude of second location + * @param lat2 [deg] latitude of second location + * @return [deg] the angular separation of the two locations. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_equ_sep() + * @sa novas_sun_angle() + * @sa novas_moon_angle() + */ +double novas_sep(double lon1, double lat1, double lon2, double lat2) { + double c = sin(lat1 * DEGREE) * sin(lat2 * DEGREE) + cos(lat1 * DEGREE) * cos(lat2 * DEGREE) * cos((lon1 - lon2) * DEGREE); + return atan2(sqrt(1.0 - c * c), c) / DEGREE; +} + +/** + * Returns the angular separation of two equatorial locations on a sphere. + * + * @param ra1 [h] right ascension of first location + * @param dec1 [deg] declination of first location + * @param ra2 [h] right ascension of second location + * @param dec2 [deg] declination of second location + * @return [deg] the angular separation of the two locations. + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_sep() + * @sa novas_sun_angle() + * @sa novas_moon_angle() + */ +double novas_equ_sep(double ra1, double dec1, double ra2, double dec2) { + return novas_sep(15.0 * ra1, dec1, 15.0 * ra2, dec2); +} + +/** + * Converts rectangular telescope x,y,z (absolute or relative) coordinates (in ITRS) toequatorial + * u,v,w projected coordinates for a specified line of sight. + * + * x,y,z are Cartesian coordinates w.r.t the Greenwitch meridian. The directions are x: long=0, lat=0; + * y: long=90, lat=0; z: lat=90. + * + * u,v,w are Cartesian coordinates (u,v) along the local equatorial R.A. and declination directions as + * seen from a direction on the sky (w). + * + * @param xyz [arb.u.] Absolute or relative x,y,z coordinates (double[3]). + * @param ha [h] Hourangle (LST - RA) i.e., the difference between the Local + * (apparent) Sidereal Time and the apparent (true-of-date) Right + * Ascension of observed source. + * @param dec [deg] Apparent (true-of-date) declination of source + * @param[out] uvw [arb.u.] Converted u,v,w coordinates (double[3]) in same units as xyz. + * It may be the same vector as the input. + * + * @return 0 if successful, or else -1 if either vector argument is NULL + * (errno will be set to EINVAL) + * + * @since 1.3 + * @author Attila Kovacs + * + * @sa novas_xyz2enu() + */ +int novas_xyz2uvw(const double *xyz, double ha, double dec, double *uvw) { + static const char *fn = "novas_xyz2uvw"; + + double dX, dY, dZ, sHA, cHA, sDec, cDec; + + if(!xyz) + return novas_error(-1, EINVAL, fn, "input xyz vector is NULL"); + if(!uvw) + return novas_error(-1, EINVAL, fn, "output uvw vector is NULL"); + + dX = xyz[0]; + dY = xyz[1]; + dZ = xyz[2]; + + ha *= HOURANGLE; + dec *= DEGREE; + + sHA = sin(ha); + cHA = cos(ha); + + sDec = sin(dec); + cDec = cos(dec); + + uvw[0] = sHA * dX + cHA * dY; + uvw[1] = -cHA * sDec * dX + sHA * sDec * dY + cDec * dZ; + uvw[2] = cHA * cDec * dX - sHA * cDec * dY + sDec * dZ; + + return 0; +} + + + + diff --git a/test/src/test-errors.c b/test/src/test-errors.c index 717768f4..a5c703d0 100644 --- a/test/src/test-errors.c +++ b/test/src/test-errors.c @@ -1603,6 +1603,210 @@ static int test_set_obsys_pole() { return n; } +static int test_lsr2ssb_vel() { + int n = 0; + double v[3] = {}; + + if(check("lsr2ssb_vel:vLSR", -1, novas_lsr2ssb_vel(NULL, 1.0, v))) n++; + if(check("lsr2ssb_vel:vSSB", -1, novas_lsr2ssb_vel(v, 1.0, NULL))) n++; + if(check("lsr2ssb_vel:unit:0", -1, novas_lsr2ssb_vel(v, 0.0, v))) n++; + if(check("lsr2ssb_vel:unit:neg", -1, novas_lsr2ssb_vel(v, -1.0, v))) n++; + + return n; +} + +static int test_hms_hours() { + int n = 0; + + if(check_nan("hms_hours:null", novas_hms_hours(NULL))) n++; + if(check_nan("hms_hours:empty", novas_hms_hours(""))) n++; + if(check_nan("hms_hours:empty", novas_hms_hours(""))) n++; + if(check_nan("hms_hours:few", novas_hms_hours("12 39"))) n++; + if(check_nan("hms_hours:dms", novas_hms_hours("12d 39m 33.0"))) n++; + if(check_nan("hms_hours:sep", novas_hms_hours("12,39,33.0"))) n++; + if(check_nan("hms_hours:min:neg", novas_hms_hours("12 -1 33.0"))) n++; + if(check_nan("hms_hours:min:60", novas_hms_hours("12 60 33.0"))) n++; + if(check_nan("hms_hours:sec:neg", novas_hms_hours("12 39 -0.1"))) n++; + if(check_nan("hms_hours:min:60", novas_hms_hours("12 39 60.0"))) n++; + + return n; +} + +static int test_dms_degrees() { + int n = 0; + + if(check_nan("dms_degrees:null", novas_dms_degrees(NULL))) n++; + if(check_nan("dms_degrees:empty", novas_dms_degrees(""))) n++; + if(check_nan("dms_degrees:empty", novas_dms_degrees(""))) n++; + if(check_nan("dms_degrees:few", novas_dms_degrees("122 39"))) n++; + if(check_nan("dms_degrees:hms", novas_dms_degrees("122h 39m 33.0"))) n++; + if(check_nan("dms_degrees:sep", novas_dms_degrees("122,39,33.0"))) n++; + if(check_nan("dms_degrees:min:neg", novas_dms_degrees("122 -1 33.0"))) n++; + if(check_nan("dms_degrees:min:60", novas_dms_degrees("122 60 33.0"))) n++; + if(check_nan("dms_degrees:sec:neg", novas_dms_degrees("122 39 -0.1"))) n++; + if(check_nan("dms_degrees:min:60", novas_dms_degrees("122 39 60.0"))) n++; + + return n; +} + +static int test_helio_dist() { + int n = 0; + double rate = 0.0; + object star = {}; + + star.type = NOVAS_CATALOG_OBJECT; + + if(check_nan("helio_dist:null", novas_helio_dist(JD_J2000, NULL, NULL))) n++; + if(check_nan("helio_dist:cat_object", novas_helio_dist(JD_J2000, &star, NULL))) n++; + if(check_nan("helio_dist:null:rate", novas_helio_dist(JD_J2000, NULL, &rate))) n++; + if(check_nan("helio_dist:null:rate:nan", rate)) n++; + if(check_nan("helio_dist:cat_object:rate", novas_helio_dist(JD_J2000, &star, &rate))) n++; + if(check_nan("helio_dist:cat_object:rate:nan", rate)) n++; + + return n; +} + +static int test_xyz2uvw() { + int n = 0; + double p[3] = {}; + + if(check("xyz2uvw:xyz", -1, novas_xyz2uvw(NULL, 0.0, 0.0, p))) n++; + if(check("xyz2uvw:uvw", -1, novas_xyz2uvw(p, 0.0, 0.0, NULL))) n++; + + return n; +} + +static int test_frame_lst() { + int n = 0; + novas_timespec time = {}; + observer obs = {}; + novas_frame frame = {}; + + if(check("frame_lst:frame:null", -1, novas_frame_lst(NULL))) n++; + if(check("frame_lst:frame:init", -1, novas_frame_lst(&frame))) n++; + + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + make_observer_at_geocenter(&obs); + if(check("frame_lst:make_frame", 0, novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + + if(check("frame_lst:obs:invalid", -1, novas_frame_lst(&frame))) n++; + + return n; +} + +static int test_rise_set() { + int n = 0; + object sun = NOVAS_SUN_INIT; + novas_timespec time = {}; + observer obs = {}; + novas_frame frame = {}; + + if(check_nan("rise_set:rises_above:frame:null", novas_rises_above(0.0, &sun, NULL, NULL))) n++; + if(check_nan("rise_set:rises_above:frame:init", novas_rises_above(0.0, &sun, &frame, NULL))) n++; + if(check_nan("rise_set:sets_below:frame:null", novas_sets_below(0.0, &sun, NULL, NULL))) n++; + if(check_nan("rise_set:sets_below:frame:init", novas_sets_below(0.0, &sun, &frame, NULL))) n++; + + // noon (near transit) + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + make_observer_on_surface(0.0, 0.0, 0.0, 0.0, 0.0, &obs); + if(check("rise_set:make_frame", 0, novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + + if(check_nan("rise_set:rises_above:source:null", novas_rises_above(0.0, NULL, &frame, NULL))) n++; + if(check_nan("rise_set:sets_below:source:null", novas_sets_below(0.0, NULL, &frame, NULL))) n++; + + return n; +} + +static int test_tracks() { + int n = 0; + object sun = NOVAS_SUN_INIT; + novas_timespec time = {}; + observer obs = {}; + novas_frame frame = {}; + novas_track track = {}; + double vel[3] = {}, x; + + if(check("equ_track:frame:null", -1, novas_equ_track(&sun, NULL, 1000.0, &track))) n++; + if(check("equ_track:frame:init", -1, novas_equ_track(&sun, &frame, 1000.0, &track))) n++; + if(check("hor_track:frame:null", -1, novas_hor_track(&sun, NULL, NULL, &track))) n++; + if(check("hor_track:frame:init", -1, novas_hor_track(&sun, &frame, NULL, &track))) n++; + + // noon (near transit) + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + make_observer_on_surface(0.0, 0.0, 0.0, 0.0, 0.0, &obs); + if(check("rise_set:make_frame", 0, novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + + if(check("equ_track:source:null", -1, novas_equ_track(NULL, &frame, 1000.0, &track))) n++; + if(check("equ_track:track:null", -1, novas_equ_track(&sun, &frame, 1000.0, NULL))) n++; + if(check("hor_track:source:null", -1, novas_hor_track(NULL, &frame, NULL, &track))) n++; + if(check("hor_track:track:null", -1, novas_hor_track(&sun, &frame, NULL, NULL))) n++; + + if(check("hor_track:make_airborne_observer", 0, make_airborne_observer(&obs.on_surf, vel, &obs))) n++; + if(check("rise_set:make_frame:airborne", 0, novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(check("hor_track:track:null", 0, novas_hor_track(&sun, &frame, NULL, &track))) n++; + + make_observer_at_geocenter(&obs); + if(check("rise_set:make_frame:geocenter", 0, novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(check("hor_track:track:null", -1, novas_hor_track(&sun, &frame, NULL, &track))) n++; + + if(check("track_pos:track:null", -1, novas_track_pos(NULL, &time, &x, NULL, NULL, NULL))) n++; + if(check("track_pos:time:null", -1, novas_track_pos(&track, NULL, &x, NULL, NULL, NULL))) n++; + + return n; +} + +int test_solar_illum() { + int n = 0; + object earth = NOVAS_EARTH_INIT; + object jupiter = NOVAS_JUPITER_INIT; + novas_timespec time = {}; + observer obs; + novas_frame frame = {}; + double pos[3] = {1.0, 1.0, 0.0}, vel[3] = {}; + + make_solar_system_observer(pos, vel, &obs); + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + + if(check_nan("solar_illum:frame:null", novas_solar_illum(&earth, NULL))) n++; + if(check_nan("solar_illum:frame:init", novas_solar_illum(&earth, &frame))) n++; + + if(check("solar_illum:make_frame", 0, novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(check_nan("solar_illum:source:null", novas_solar_illum(NULL, &frame))) n++; + if(check_nan("solar_illum:source:invalid", novas_solar_illum(&jupiter, &frame))) n++; + + return n; +} + +int test_object_sep() { + int n = 0; + object sun = NOVAS_SUN_INIT; + object earth = NOVAS_EARTH_INIT; + object jupiter = NOVAS_JUPITER_INIT; + novas_timespec time = {}; + observer obs; + novas_frame frame = {}; + + make_observer_at_geocenter(&obs); + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + + if(check_nan("object_sep:frame:null", novas_object_sep(&sun, &sun, NULL))) n++; + if(check_nan("object_sep:frame:init", novas_object_sep(&sun, &sun, &frame))) n++; + + if(check("object_sep:make_frame", 0, novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + + if(check_nan("object_sep:a:null", novas_object_sep(NULL, &sun, &frame))) n++; + if(check_nan("object_sep:b:null", novas_object_sep(&sun, NULL, &frame))) n++; + + if(check_nan("object_sep:a:invalid", novas_object_sep(&jupiter, &sun, &frame))) n++; + if(check_nan("object_sep:b:invalid", novas_object_sep(&sun, &jupiter, &frame))) n++; + + if(check_nan("object_sep:a=obs", novas_object_sep(&earth, &sun, &frame))) n++; + if(check_nan("object_sep:b=obs", novas_object_sep(&sun, &earth, &frame))) n++; + if(check_nan("object_sep:a=b=obs", novas_object_sep(&earth, &earth, &frame))) n++; + + return n; +} + int main() { int n = 0; @@ -1741,6 +1945,18 @@ int main() { if(test_gcrs_to_mod()) n++; if(test_mod_to_gcrs()) n++; + if(test_lsr2ssb_vel()) n++; + if(test_hms_hours()) n++; + if(test_dms_degrees()) n++; + if(test_helio_dist()) n++; + if(test_xyz2uvw()) n++; + + if(test_frame_lst()) n++; + if(test_rise_set()) n++; + if(test_tracks()) n++; + if(test_solar_illum()) n++; + if(test_object_sep()) n++; + if(n) fprintf(stderr, " -- FAILED %d tests\n", n); else fprintf(stderr, " -- OK\n"); diff --git a/test/src/test-super.c b/test/src/test-super.c index cddb7d01..b3fbc989 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -2423,6 +2423,369 @@ static int test_orbit_posvel_callisto() { return n; } +static int test_hms_hours() { + int n = 0; + double hours = 23.0 + 59.0/60.0 + 59.999/3600.0; + + if(!is_equal("hms_hours:colons", novas_hms_hours("23:59:59.999"), hours, 1e-10)) n++; + if(!is_equal("hms_hours:spaces", novas_hms_hours("23 59 59.999"), hours, 1e-10)) n++; + if(!is_equal("hms_hours:hm", novas_hms_hours("23h59m59.999"), hours, 1e-10)) n++; + if(!is_equal("hms_hours:HM", novas_hms_hours("23H59M59.999"), hours, 1e-10)) n++; + if(!is_equal("hms_hours:hprime", novas_hms_hours("23h59'59.999"), hours, 1e-10)) n++; + if(!is_equal("hms_hours:combo", novas_hms_hours("23h 59' 59.999"), hours, 1e-10)) n++; + + return n; +} + +static int test_dms_degrees() { + int n = 0; + double degs = 179.0 + 59.0/60.0 + 59.999/3600.0; + + if(!is_equal("dms_degrees:colons", novas_dms_degrees("179:59:59.999"), degs, 1e-9)) n++; + if(!is_equal("dms_degrees:spaces", novas_dms_degrees("179 59 59.999"), degs, 1e-9)) n++; + if(!is_equal("dsm_degrees:dm", novas_dms_degrees("179d59m59.999"), degs, 1e-9)) n++; + if(!is_equal("dms_degrees:DM", novas_dms_degrees("179D59M59.999"), degs, 1e-9)) n++; + if(!is_equal("dms_degrees:dprime", novas_dms_degrees("179d59'59.999"), degs, 1e-9)) n++; + if(!is_equal("dms_degrees:combo", novas_dms_degrees("179d 59' 59.999"), degs, 1e-9)) n++; + + if(!is_equal("dms_degrees:neg:colons", novas_dms_degrees("-179:59:59.999"), -degs, 1e-9)) n++; + if(!is_equal("dms_degrees:neg:spaces", novas_dms_degrees("-179 59 59.999"), -degs, 1e-9)) n++; + if(!is_equal("dsm_degrees:neg:dm", novas_dms_degrees("-179d59m59.999"), -degs, 1e-9)) n++; + if(!is_equal("dms_degrees:neg:DM", novas_dms_degrees("-179D59M59.999"), -degs, 1e-9)) n++; + if(!is_equal("dms_degrees:neg:dprime", novas_dms_degrees("-179d59'59.999"), -degs, 1e-9)) n++; + if(!is_equal("dms_degrees:neg:combo", novas_dms_degrees("-179d 59' 59.999"), -degs, 1e-9)) n++; + + return n; +} + +static int test_hpa() { + int n = 0; + + if(!is_equal("hpa:S", novas_hpa(180, 60, 45), 0.0, 1e-9)) n++; + if(!is_equal("hpa:E", novas_hpa(90, 60, 0), -90.0, 1e-9)) n++; + if(!is_equal("hpa:W", novas_hpa(-90, 60, 0), 90.0, 1e-9)) n++; + if(!is_equal("hpa:N1", remainder(novas_hpa(0, 60, 45) - 180.0, 360.0), 0.0, 1e-9)) n++; + if(!is_equal("hpa:N2", novas_hpa(0, 30, 45), 0.0, 1e-9)) n++; + + return n; +} + +static int test_epa() { + int n = 0; + + if(!is_equal("epa:ra=0:transit:S", novas_epa(0, 30, 45), 0.0, 1e-9)) n++; + if(!is_equal("epa:ra=0:transit:N", remainder(novas_epa(0, 60, 45) - 180.0, 360.0), 0.0, 1e-9)) n++; + if(!is_equal("epa:ra=0:rise", novas_epa(-6, 30, 0), -90.0, 1e-9)) n++; + if(!is_equal("epa:ra=0:set", novas_epa(6, 30, 0), 90.0, 1e-9)) n++; + + return n; +} + + +static int test_helio_dist() { + int n = 0; + object earth = NOVAS_EARTH_INIT; + object sun = NOVAS_SUN_INIT; + double rate; + + if(!is_equal("helio_dist:earth", novas_helio_dist(NOVAS_JD_J2000, &earth, &rate), 1.0, 0.03)) n++; + if(!is_equal("helio_dist:earth:rate", rate, 0.0, 0.03)) n++; + if(!is_equal("helio_dist:earth:rate:NULL", novas_helio_dist(NOVAS_JD_J2000, &earth, NULL), 1.0, 0.03)) n++; + + if(!is_equal("helio_dist:sun", novas_helio_dist(NOVAS_JD_J2000, &sun, &rate), 0.0, 1e-9)) n++; + if(!is_equal("helio_dist:sun:rate", rate, 0.0, 1e-9)) n++; + if(!is_equal("helio_dist:sun:rate:NULL", novas_helio_dist(NOVAS_JD_J2000, &sun, NULL), 0.0, 1e-9)) n++; + + return n; +} + +static int test_solar_power() { + int n = 0; + object earth = NOVAS_EARTH_INIT; + + if(!is_equal("solar_power:earth", novas_solar_power(NOVAS_JD_J2000, &earth), 1360.8, 130.0)) n++; + + return n; +} + +static int test_solar_illum() { + int n = 0; + object cat, earth = NOVAS_EARTH_INIT; + novas_timespec time = {}; + observer obs; + novas_frame frame = {}; + double pos[3] = {}, vel[3] = {}; + int i; + + make_redshifted_object("test", 0.0, 0.0, 0.0, &cat); + make_solar_system_observer(pos, vel, &obs); + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + + if(!is_ok("solar_illum:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(!is_equal("solar_illum:source:sidereal", 1.0, novas_solar_illum(&cat, &frame), 1e-12)) n++; + if(!is_equal("solar_illum:source:earth:ssb", 1.0, novas_solar_illum(&earth, &frame), 1e-3)) n++; + + for(i = 3; --i >= 0; ) obs.near_earth.sc_pos[i] = 1.1 * frame.earth_pos[i]; + if(!is_ok("solar_illum:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(!is_equal("solar_illum:source:earth:beyond", 0.0, novas_solar_illum(&earth, &frame), 1e-3)) n++; + + for(i = 3; --i >= 0; ) obs.near_earth.sc_pos[i] = frame.earth_pos[i] + frame.earth_vel[i]; + if(!is_ok("solar_illum:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(!is_equal("solar_illum:source:earth:beyond", 0.5, novas_solar_illum(&earth, &frame), 1e-3)) n++; + + return n; +} + +static int test_equ_sep() { + int n = 0; + + if(!is_equal("equ_sep:dec=0:ra+1", novas_equ_sep(5.5, 0.0, 6.5, 0.0), 15.0, 1e-9)) n++; + if(!is_equal("equ_sep:dec=0:ra-1", novas_equ_sep(5.5, 0.0, 6.5, 0.0), 15.0, 1e-9)) n++; + if(!is_equal("equ_sep:dec=60:ra+0.01", novas_equ_sep(5.5, 60.0, 5.51, 60.0), 0.075, 1e-5)) n++; + if(!is_equal("equ_sep:dec+1", novas_equ_sep(5.5, 15.3, 5.5, 16.3), 1.0, 1e-9)) n++; + if(!is_equal("equ_sep:poles", novas_equ_sep(1.0, -90.0, 3.0, 90.0), 180.0, 1e-9)) n++; + if(!is_equal("equ_sep:pole:equ", novas_equ_sep(1.0, -90.0, 3.0, 0.0), 90.0, 1e-9)) n++; + return n; +} + +static int test_h2e_offset() { + int pa, n = 0; + + for(pa = -180; pa <= 180; pa += 15) { + double s = sin(pa * DEGREE); + double c = cos(pa * DEGREE); + int daz; + + for(daz = -100; daz < 100; daz += 10) { + int del; + for(del = -100; del <= 100; del += 10) { + char label[80]; + double dra, ddec, dAZ, dEL; + + sprintf(label, "h2e_offset:PA=%d:az=%d:el=%d:dra", pa, daz, del); + novas_h2e_offset(daz, del, pa, &dra, NULL); + if(!is_equal(label, dra, -c * daz + s * del, 1e-9)) n++; + + sprintf(label, "h2e_offset:PA=%d:az=%d:el=%d:ddec", pa, daz, del); + novas_h2e_offset(daz, del, pa, NULL, &ddec); + if(!is_equal(label, ddec, s * daz + c * del, 1e-9)) n++; + + sprintf(label, "h2e_offset:PA=%d:az=%d:el=%d:daz", pa, daz, del); + novas_e2h_offset(dra, ddec, pa, &dAZ, NULL); + if(!is_equal(label, dAZ, daz, 1e-9)) n++; + + sprintf(label, "h2e_offset:PA=%d:az=%d:el=%d:del", pa, daz, del); + novas_e2h_offset(dra, ddec, pa, NULL, &dEL); + if(!is_equal(label, dEL, del, 1e-9)) n++; + } + } + } + + return n; +} + +static int test_object_sep() { + int n = 0; + object a = {}, b = {}; + novas_timespec time = {}; + observer obs = {}; + novas_frame frame = {}; + + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + make_observer_at_geocenter(&obs); + if(!is_ok("object_sep:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + + make_redshifted_object("a", 0.0, 60.0, 0.0, &a); + make_redshifted_object("b", 0.01, 60.0, 0.0, &b); + + if(!is_equal("object_sep:same", novas_object_sep(&a, &a, &frame), 0.0, 1e-12)) n++; + + if(!is_equal("object_sep:dra=0.01:dec=60", novas_object_sep(&a, &b, &frame), 0.075, 1e-5)) n++; + + b.star.ra = 0.0; + b.star.dec = 61.0; + if(!is_equal("object_sep:ddec=1:ra=1", novas_object_sep(&a, &b, &frame), 1.0, 1e-4)) n++; + + b.star.ra = 0.02/15.0; + b.star.dec = 60.01; + if(!is_equal("object_sep:ddra=ddec=0.01", novas_object_sep(&a, &b, &frame), 0.01 * sqrt(2.0), 1e-4)) n++; + + + return n; +} + +static int test_frame_lst() { + int n = 0; + observer obs = {}; + novas_timespec time = {}; + novas_frame frame = {}; + on_surface loc = {}; + double vel[3] = {}; + + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + make_observer_on_surface(33.0, 15.0, 0.0, 0.0, 0.0, &obs); + if(!is_ok("frame_lst:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(!is_equal("frame_lst:lst", novas_frame_lst(&frame), frame.gst + 1.0, 1e-9)) n++; + + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000 - 0.5, 32.0, 0.0, &time); + make_observer_on_surface(33.0, 15.0, 0.0, 0.0, 0.0, &obs); + if(!is_ok("frame_lst:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(!is_equal("frame_lst:lst", novas_frame_lst(&frame), frame.gst + 13.0, 1e-9)) n++; + + loc = obs.on_surf; + make_airborne_observer(&loc, vel, &obs); + if(!is_ok("frame_lst:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(!is_equal("frame_lst:lst", novas_frame_lst(&frame), frame.gst + 13.0, 1e-9)) n++; + + return n; +} + +static int test_rise_set() { + int n = 0; + observer obs = {}; + novas_timespec time = {}; + novas_frame frame = {}; + object sun = NOVAS_SUN_INIT; + double refr; + + // noon (near transit) + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + make_observer_on_surface(0.0, 0.0, 0.0, 0.0, 0.0, &obs); + if(!is_ok("rise_set:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + + // 6am next day... + if(!is_equal("rise_set:rise", novas_rises_above(0.0, &sun, &frame, NULL), NOVAS_JD_J2000 + 0.75, 0.01)) n++; + + // 6pm same day... + if(!is_equal("rise_set:set", novas_sets_below(0.0, &sun, &frame, NULL), NOVAS_JD_J2000 + 0.25, 0.01)) n++; + + refr = refract_astro(&obs.on_surf, NOVAS_STANDARD_ATMOSPHERE, 90.0); + + // 6am next day... + if(!is_equal("rise_set:rise", novas_rises_above(refr, &sun, &frame, novas_standard_refraction), NOVAS_JD_J2000 + 0.75, 0.01)) n++; + + // 6pm same day... + if(!is_equal("rise_set:set", novas_sets_below(refr, &sun, &frame, novas_standard_refraction), NOVAS_JD_J2000 + 0.25, 0.01)) n++; + + return n; +} + +static int test_equ_track() { + int n = 0; + observer obs = {}; + novas_timespec time = {}; + novas_frame frame = {}; + object sun = NOVAS_SUN_INIT; + novas_track track = {}; + sky_pos pos = {}; + double x = 0.0; + + // noon (near transit) + novas_set_time(NOVAS_TDB, NOVAS_JD_J2000, 32.0, 0.0, &time); + make_observer_on_surface(0.0, 0.0, 0.0, 0.0, 0.0, &obs); + if(!is_ok("equ_track:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + + if(!is_ok("equ_track:sky_pos", novas_sky_pos(&sun, &frame, NOVAS_TOD, &pos))) n++; + if(!is_ok("equ_track", novas_equ_track(&sun, &frame, 3600.0, &track))) n++; + + if(!is_equal("equ_track:ra", track.pos.lon / 15.0, pos.ra, 1e-9)) n++; + if(!is_equal("equ_track:dec", track.pos.lat, pos.dec, 1e-9)) n++; + if(!is_equal("equ_track:dis", track.pos.dist, pos.dis, 1e-9)) n++; + if(!is_equal("equ_track:z", track.pos.z, novas_v2z(pos.rv), 1e-9)) n++; + + if(!is_equal("equ_track:rate", hypot(track.rate.lon, track.rate.lat), (360.0 / 365.25) / DAY, 0.2 / DAY)) n++; + if(!is_equal("equ_track:rate:z", track.rate.z, 0.0, 1e-9)) n++; + if(!is_equal("equ_track:rate:dist", track.rate.dist, 0.0, 1e-9)) n++; + + if(!is_equal("equ_track:accel", hypot(track.accel.lon, track.accel.lat), 0.0, 0.03 / (DAY * DAY))) n++; + if(!is_equal("equ_track:accel:z", track.accel.z, 0.0, 1e-16)) n++; + if(!is_equal("equ_track:accel:dist", track.accel.dist, 0.0, 1e-12)) n++; + + time.fjd_tt += 0.01; + if(!is_ok("equ_track:make_frame:shifted", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(!is_ok("equ_track:sky_pos", novas_sky_pos(&sun, &frame, NOVAS_TOD, &pos))) n++; + + if(!is_ok("equ_track:track_pos:lon", novas_track_pos(&track, &time, &x, NULL, NULL, NULL))) n++; + if(!is_equal("equ_track:track_pos:lon:check", x, remainder(15.0 * pos.ra, 360.0), 1e-5)) n++; + + if(!is_ok("equ_track:track_pos:lat", novas_track_pos(&track, &time, NULL, &x, NULL, NULL))) n++; + if(!is_equal("equ_track:track_pos:lat:check", x, pos.dec, 1e-5)) n++; + + if(!is_ok("equ_track:track_pos:dist", novas_track_pos(&track, &time, NULL, NULL, &x, NULL))) n++; + if(!is_equal("equ_track:track_pos:dist:check", x, pos.dis, 1e-9)) n++; + + if(!is_ok("equ_track:track_pos:z", novas_track_pos(&track, &time, NULL, NULL, NULL, &x))) n++; + if(!is_equal("equ_track:track_pos:dist:z", x, novas_v2z(pos.rv), 1e-9)) n++; + + return n; +} + + +static int test_hor_track() { + int n = 0; + observer obs = {}; + novas_timespec time = {}; + novas_frame frame = {}; + object source = {}; + novas_track track = {}; + sky_pos pos = {}; + double az0 = 0.0, el0 = 0.0, x = 0.0; + + // noon (near transit) + novas_set_time(NOVAS_TT, NOVAS_JD_J2000, 32.0, 0.0, &time); + make_observer_on_surface(0.0, 0.0, 0.0, 0.0, 0.0, &obs); + if(!is_ok("hor_track:make_frame", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + + // A source that transits at 60 deg elevation + make_redshifted_object("Test", frame.gst, -60.0, 0.0, &source); + + if(!is_ok("hor_track:sky_pos", novas_sky_pos(&source, &frame, NOVAS_TOD, &pos))) n++; + if(!is_ok("hor_track:app_to_hor", novas_app_to_hor(&frame, NOVAS_TOD, pos.ra, pos.dec, NULL, &az0, &el0))) n++; + if(!is_ok("hor_track", novas_hor_track(&source, &frame, NULL, &track))) n++; + + if(!is_equal("hor_track:ra", track.pos.lon, az0, 1e-9)) n++; + if(!is_equal("hor_track:dec", track.pos.lat, el0, 1e-9)) n++; + if(!is_equal("hor_track:dis", track.pos.dist, pos.dis, 1e-12 * pos.dis)) n++; + if(!is_equal("hor_track:z", track.pos.z, novas_v2z(pos.rv), 1e-9)) n++; + + if(!is_equal("hor_track:rate:lat", track.rate.lat, 0.0, 1e-5)) n++; + if(!is_equal("hor_track:rate:z", track.rate.z, 0.0, 1e-9)) n++; + if(!is_equal("hor_track:rate:dist", track.rate.dist, 0.0, 1e-3)) n++; + + if(!is_equal("hor_track:accel:lon", track.accel.lon, 0.0, 1e-9)) n++; + if(!is_equal("hor_track:rate:lat", track.rate.lat, 0.0, 1e-3)) n++; + if(!is_equal("hor_track:accel:z", track.accel.z, 0.0, 1e-16)) n++; + if(!is_equal("hor_track:accel:dist", track.accel.dist, 0.0, 1.0)) n++; + + if(!is_ok("hor_track:app_to_hor:ref", novas_app_to_hor(&frame, NOVAS_TOD, pos.ra, pos.dec, novas_standard_refraction, &az0, &el0))) n++; + if(!is_ok("hor_track:ref", novas_hor_track(&source, &frame, novas_standard_refraction, &track))) n++; + + if(!is_equal("hor_track:ra:ref", track.pos.lon, az0, 1e-9)) n++; + if(!is_equal("hor_track:dec:ref", track.pos.lat, el0, 1e-9)) n++; + if(!is_equal("hor_track:dis:ref", track.pos.dist, pos.dis, 1e-12 * pos.dis)) n++; + if(!is_equal("hor_track:z:ref", track.pos.z, novas_v2z(pos.rv), 1e-9)) n++; + + time.fjd_tt += 10.0 / DAY; + if(!is_ok("hor_track:make_frame:shifted", novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &time, 0.0, 0.0, &frame))) n++; + if(!is_ok("hor_track:sky_pos", novas_sky_pos(&source, &frame, NOVAS_TOD, &pos))) n++; + if(!is_ok("hor_track:app_to_hor", novas_app_to_hor(&frame, NOVAS_TOD, pos.ra, pos.dec, novas_standard_refraction, &az0, &el0))) n++; + + if(!is_ok("hor_track:track_pos:lon", novas_track_pos(&track, &time, &x, NULL, NULL, NULL))) n++; + if(!is_equal("hor_track:track_pos:lon:check", x, remainder(az0, 360.0), 1e-3)) n++; + + if(!is_ok("hor_track:track_pos:lat", novas_track_pos(&track, &time, NULL, &x, NULL, NULL))) n++; + if(!is_equal("hor_track:track_pos:lat:check", x, el0, 1e-3)) n++; + + if(!is_ok("hor_track:track_pos:dist", novas_track_pos(&track, &time, NULL, NULL, &x, NULL))) n++; + if(!is_equal("hor_track:track_pos:dist:check", x, pos.dis, 1e-12 * pos.dis)) n++; + + if(!is_ok("hor_track:track_pos:z", novas_track_pos(&track, &time, NULL, NULL, NULL, &x))) n++; + if(!is_equal("hor_track:track_pos:dist:z", x, novas_v2z(pos.rv), 1e-9)) n++; + + return n; +} + + int main(int argc, char *argv[]) { int n = 0; @@ -2492,6 +2855,22 @@ int main(int argc, char *argv[]) { if(test_orbit_place()) n++; if(test_orbit_posvel_callisto()) n++; + if(test_hms_hours()) n++; + if(test_dms_degrees()) n++; + if(test_hpa()) n++; + if(test_epa()) n++; + if(test_helio_dist()) n++; + if(test_solar_power()) n++; + if(test_solar_illum()) n++; + if(test_equ_sep()) n++; + if(test_object_sep()) n++; + if(test_h2e_offset()) n++; + + if(test_frame_lst()) n++; + if(test_rise_set()) n++; + if(test_equ_track()) n++; + if(test_hor_track()) n++; + n += test_dates(); return n;