diff --git a/.gitignore b/.gitignore index 64383049..f865ba37 100644 --- a/.gitignore +++ b/.gitignore @@ -19,3 +19,4 @@ README-headless.md README-orig.md TODO VERSION +RELEASE-HOWTO.md diff --git a/CHANGELOG.md b/CHANGELOG.md index 65eab599..7ab3e1d6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -7,47 +7,36 @@ The format is based on [Keep a Changelog](https://keepachangelog.com/en/1.1.0/), -## [Unreleased] +## [1.1.0] -- 2024-08-11 -Changes coming to the next quarterly release, expected around 1 September 2024. Some or all of these may be readily -available on the `main` branch. +Feature release. Introducing a more efficient and elegant approach to position and velocity calculations using +observer frames, versatile handling of astronomical timescales, support for further observer locations, coordinate +reference systems and atmospheric refraction models. The release also fixes a number of bugs, of varying severity, +which affected prior SuperNOVAS releases. ### Fixed - - #41: `grav_def()` gravitating body position antedated somewhat incorrectly (in v1.0) when observed source is a - Solar-system object between the observer and the gravitating body. The resulting positional error is typically - small at below 10 uas. - - - #39: `tod_to_itrs()` used wrong Earth rotation measure (`NOVAS_ERA` instead of `NOVAS_GST`). + - #29: Fix portability to non-Intel platforms. Previously, SuperNOVAS used `char` for storing integer coefficients, + assuming `char` was signed. However, on some platforms like ARM and PowerPC `char` is unsigned, which broke + calculations badly on such platforms. As of now, we use the explicit platform independent signed `int8_t` storage + type for these coefficients. - #38: `gcrs_to_j2000` transformed in the wrong direction. - - #37: `gcrs_to_cirs()` did not handle well if input and output vectors were the same. + - #39: `tod_to_itrs()` used wrong Earth rotation measure (`NOVAS_ERA` instead of `NOVAS_GST`). - #36: `tt2tdb()` Had a wrong scaling in sinusoidal period, resulting in an error of up to +/- 1.7 ms. - - - #34: Radial velocity calculation to precede aberration and gravitational bending in `place()`, since the radial - velocity that is observed is in the geometric direction towards the source (unaffected by aberration). A precise - accounting of the gravitational effects would require figuring out the direction in which the observed light was - emitted from the source before it was bent by gravitating bodies along the way. In practice, this may be difficult - to generalize. For a single gravitationg body the geometric direction of the source is between the direction in - which the light is emitted, and the observed deflected direction. Therefore, for the time being, the radial - velocity calculated via the geometric direction is closer to the actual value. - - - #29: Fix portability to non-Intel x86 platforms (see Issue #29). Previously, SuperNOVAS used `char` for storing - integer coefficients, assuming `char` was signed. However, on some platforms like ARM and PowerPC `char` is - unsigned, which broke many calculations badly for such platforms. As of now, we use the explicit platform - independent `int8_t` storage type for these coefficients. + + - #37: `gcrs_to_cirs()` did not handle well if input and output vectors were the same. - #28: Division by zero bug in `d_light()` (since NOVAS C 3.1) if the first position argument is the ephemeris reference position (e.g. the Sun for `solsys3.c`). The bug affects for example `grav_def()`, where it effectively results in the gravitational deflection due to the Sun being skipped. See Issue #28. - - - `PSI_COR` and `EPS_COR` made globally visible again, thus improving NOVAS C 3.1 compatibility. - - - Adjusted regression testing to treat `nan` and `-nan` effectively the same. They both represent an equally invalid - result regardless of the sign. + + - #41: `grav_def()` gravitating body position antedated somewhat incorrectly (in v1.0) when observed source is a + Solar-system object between the observer and the gravitating body. The resulting positional error is typically + small at below 10 uas. - #24: Bungled definition of `SUPERNOVAS_VERSION_STRING` in `novas.h`. @@ -58,10 +47,10 @@ available on the `main` branch. - #33: New observing-frame based approach for calculations (`frames.c`). A `novas_frame` object uniquely defines both the place and time of observation, with a set of pre-calculated transformations and constants. Once the frame is - defined it can be used very efficiently to calculate positions for multiple celestial objects with minimal - additional computational cost. The frames API is also more elegant and simpler than the low-level NOVAS C approach - for performing the same kind of calculations. And, frames are inherently thread-safe since post-creation their - internal state is never modified during the calculations. The following new functions were added: + defined it can be used very efficiently to calculate positions for multiple celestial objects with minimum + additional computational cost. The frames API is also more elegant and more versatile than the low-level NOVAS C + approach for performing the same kind of calculations. And, frames are inherently thread-safe since post-creation + their internal state is never modified during the calculations. The following new functions were added: `novas_make_frame()`, `novas_change_observer()`, `novas_geom_posvel()`, `novas_geom_to_app()`, `novas_sky_pos()`, `novas_app_to_hor()`, `novas_app_to_geom()`, `novas_hor_to_app()`, `novas_make_transform()`, `novas_invert_transform()`, `novas_transform_vector()`, and `novas_transform_sky_pos()`. @@ -71,14 +60,18 @@ available on the `main` branch. (UTC, UT1, GPS, TAI, TT, TCG, TDB, or TCB), or to a UNIX time with `novas_set_unix_time()`. Once set, you can obtain an expression of that time in any timescale of choice via `novas_get_time()`, `novas_get_split_time()` or `novas_get_unix_time()`. And, you can create a new time specification by incrementing an existing one, using - `novas_increment_time()`, or measure time differences via `novas_diff_time()`. - + `novas_increment_time()`, or measure time differences via `novas_diff_time()`, `novas_diff_tcg()`, or + `novas_diff_tcb()`. + + - Added `novas_planet_bundle` structure to handle planet positions and velocities more elegantly (e.g. for + gravitational deflection calculations). + - #32: Added `grav_undef()` to undo gravitational bending of the observed light to obtain geometric positions from observed ones. - Added `obs_posvel()` to calculate the observer position and velocity relative to the Solar System Barycenter (SSB). - - Added `obs_planets()` to calculate planet positions (relative to observer) and velocities (w.r.t. SSB). + - Added `obs_planets()` to calculate apparent planet positions (relative to observer) and velocities (w.r.t. SSB). - Added new observer locations `NOVAS_AIRBORNE_OBSERVER` for an observer moving relative to the surface of Earth e.g. in an aircraft or balloon based telescope platform, and `NOVAS_SOLAR_SYSTEM_OBSERVER` for spacecraft orbiting the @@ -98,17 +91,19 @@ available on the `main` branch. - Added humidity field to `on_surface` structure, e.g. for refraction calculations at radio wavelengths. The `make_on_surface()` function will set humidity to 0.0, but the user can set the field appropriately afterwards. - - New set of built-in refraction models to use with the frame-based `novas_app_to_hor()` function. The models - `novas_standard_refraction()` and `novas_optical_refraction()` implement the same refraction model as `refract()` - in NOVAS C 3.1, with `NOVAS_STANDARD_ATMOSPHERE` and `NOVAS_WEATHER_AT_LOCATION` respectively, including the - reversed direction provided by `refract_astro()`. The user may supply their own refraction models to - `novas_app_to_hor()` also, and may make used of the generic reversal function `novas_inv_refract` to calculate - refraction in the reverse direction (observer vs astrometric elevations) as needed. + - New set of built-in refraction models to use with the frame-based `novas_app_to_hor()` / `novas_hor_to_app()` + functions. The models `novas_standard_refraction()` and `novas_optical_refraction()` implement the same refraction + model as `refract()` in NOVAS C 3.1, with `NOVAS_STANDARD_ATMOSPHERE` and `NOVAS_WEATHER_AT_LOCATION` + respectively, including the reversed direction provided by `refract_astro()`. The user may supply their own custom + refraction also, and may make use of the generic reversal function `novas_inv_refract()` to calculate refraction in + the reverse direction (observer vs astrometric elevations) as needed. - Added radio refraction model `novas_radio_refraction()` based on the formulae by Berman & Rockwell 1976. - - #42: `cio_array()` can now parse the original ASCII CIO locator data file (`data/CIO_RA.TXT`) efficiently also, - thus no longer requiring a platform-specific binary translation via the `cio_file` tool. + - Added `cirs_to_tod()` and `tod_to_cirs()` functions for efficient tranformation between True of Date (TOD) and + Celestial Intermediate Reference System (CIRS), and vice versa. + + - Added `make_cat_object()` function to create a NOVAS celestial `object` structure from existing `cat_entry` data. - `make help` to provide a brief list and explanation of the available build targets. (Thanks to `@teuben` for suggesting this.) @@ -119,6 +114,17 @@ available on the `main` branch. ### Changed + - #34: Radial velocity calculation to precede aberration and gravitational bending in `place()`, since the radial + velocity that is observed is in the geometric direction towards the source (unaffected by aberration). A precise + accounting of the gravitational effects would require figuring out the direction in which the observed light was + emitted from the source before it was bent by gravitating bodies along the way. In practice, this may be difficult + to generalize. For a single gravitationg body the geometric direction of the source is between the direction in + which the light is emitted, and the observed deflected direction. Therefore, for the time being, the radial + velocity calculated via the geometric direction is closer to the actual value. + + - #42: `cio_array()` can now parse the original ASCII CIO locator data file (`data/CIO_RA.TXT`) efficiently also, + thus no longer requiring a platform-specific binary translation via the `cio_file` tool. + - `cio_file` tool parses interval from header rather than the less precise differencing of the first two record timestamps. This leads to `cio_array()` being more accurately centered on matching date entries, e.g. J2000. @@ -136,12 +142,21 @@ available on the `main` branch. about the order in which terms are accumulated and combined, resulting in a small improvement on the few uas (micro-arcsecond) level. - - The `ra` or `dec` arguments passed to `vector2radec()` may now be NULL if not required. + - `vector2radec()`: `ra` or `dec` arguments may now be NULL if not required. - - `tt2tdb()` Now uses the same more precise series as the original NOVAS C `tdb2tt()`. + - `tt2tdb()` Now uses the same, slightly more precise series as the original NOVAS C `tdb2tt()`. + + - `PSI_COR` and `EPS_COR` made globally visible again, thus improving NOVAS C 3.1 compatibility. + + - Convergent inverse calculations now use the `novas_inv_max_iter` variable declared in `novas.c` to specify the + maximum number of iterations before inverse functions return with an error (with errno set to `ECANCELED`). Users + may adjust this limit, if they prefer some other maximum value. + + - Adjusted regression testing to treat `nan` and `-nan` effectively the same. They both represent an equally invalid + result regardless of the sign. - The default make target is now `distro`. It's similar to the deprecated `api` target from before except that it - skips building `static` libraries. + skips building `static` libraries and `cio_ra.bin`. - `make` now generates `.so` shared libraries with `SONAME` set to `lib.so.1` where the `.1` indicates that it is major version 1 of the `ABI`. All 1.x.x releases are expected to be ABI compatible with earlier releases. @@ -170,7 +185,7 @@ available on the `main` branch. - Doxygen tag file renamed to `supernovas.tag` for consistency. - - Initialize test variable for reproducibility + - Initialize test variables for reproducibility @@ -231,7 +246,7 @@ from which SuperNOVAS is forked from. requested mean equinox of date coordinates. - Some remainder calculations in NOVAS C 3.1 used the result from `fmod()` unchecked, which resulted in angles outside - of the expected [0:π] range and was also the reason why `cal_date()` did not work for negative JD values. + of the expected [0:2π] range and was also the reason why `cal_date()` did not work for negative JD values. - Fixes NOVAS C 3.1 `aberration()` returning NaN vectors if the `ve` argument is 0. It now returns the unmodified input vector appropriately instead. diff --git a/Doxyfile b/Doxyfile index 789f74eb..19133b2d 100644 --- a/Doxyfile +++ b/Doxyfile @@ -48,7 +48,7 @@ PROJECT_NAME = SuperNOVAS # could be handy for archiving the generated documentation or if some version # control system is used. -PROJECT_NUMBER = v1.0 +PROJECT_NUMBER = v1.1 # Using the PROJECT_BRIEF tag one can provide an optional one line description # for a project that appears at the top of each page and should give viewer a diff --git a/LEGACY.md b/LEGACY.md index ec87d824..6e796747 100644 --- a/LEGACY.md +++ b/LEGACY.md @@ -1,11 +1,22 @@ +# SuperNOVAS: Astrometric Positions the Old Way -### Calculating positions for a sidereal source +As of version 1.1, the SuperNOVAS library offers a new, more versatile, more intuitive, and more efficient way to +calculate the astrometric positions (and velocities) of celestial sources, via observing frames (see `README.md`). +However the old approach of the NOVAS C library remain viable also. This document demonstrates calculating the +astrometric places of sources the old way, without using the observing frames approach that is now preferred in +SuperNOVAS. + - [Calculating positions for a sidereal source](#old-sidereal-example) + - [Calculating positions for a Solar-system source](#old-solsys-example) + + + +## Calculating positions for a sidereal source A sidereal source may be anything beyond the solar-system with 'fixed' catalog coordinates. It may be a star, or a galactic molecular cloud, or a distant quasar. -#### Specify the object of interest +### Specify the object of interest First, you must provide the coordinates (which may include proper motion and parallax). Let's assume we pick a star for which we have B1950 (i.e. FK4) coordinates: @@ -34,7 +45,7 @@ adjustment to convert from J2000 to ICRS coordinates. (Naturally, you can skip the transformation steps above if you have defined your source in ICRS coordinates from the start.) -#### Spefify the observer location +### Spefify the observer location Next, we define the location where we observe from. Here we can (but don't have to) specify local weather parameters (temperature and pressure) also for refraction correction later (in this example, we'll skip the weather): @@ -47,7 +58,7 @@ Next, we define the location where we observe from. Here we can (but don't have make_observer_on_surface(50.7374, 7.0982, 60.0, 0.0, 0.0, &obs); ``` -#### Specify the time of observation +### Specify the time of observation We also need to set the time of observation. Our clocks usually measure UTC, but for astrometry we usually need time measured based on Terrestrial Time (TT) or Barycentric Time (TDB) or UT1. For a ground-based observer, you will often @@ -70,7 +81,7 @@ UT1 - UTC time difference (a.k.a. DUT1): double ut1_to_tt = get_ut1_to_tt(leap_seconds, dut1); ``` -#### Specify Earth orientation parameters +### Specify Earth orientation parameters Next, you may want to set the small diurnal (sub-arcsec level) corrections to Earth orientation, which are published in the [IERS Bulletins](https://www.iers.org/IERS/EN/Publications/Bulletins/bulletins.html). The obvious utility of @@ -87,7 +98,7 @@ early on: cel_pole(jd_tt, POLE_OFFSETS_X_Y, dx, dy); ``` -#### Calculate apparent positions on sky +### Calculate apparent positions on sky Now we can calculate the precise apparent position (CIRS or TOD) of the source, such as it's right ascension (R.A.) and declination, and the equatorial _x,y,z_ unit vector pointing in the direction of the source (in the requested @@ -116,7 +127,7 @@ the given date/time of observation. We may use it to get true apparent R.A. and azimuth and elevation at the observing location. We'll consider these two cases separately below. -##### A. True apparent R.A. and declination +#### A. True apparent R.A. and declination If you want to know the apparent R.A. and declination coordinates from the `sky_pos` structure you obtained, then you can follow with: @@ -140,7 +151,8 @@ not the true equinox of date. Thus, we must correct for the difference of the or ra = cirs_to_app_ra(jd_tt, NOVAS_FULL_ACCURACY, ra); ``` -##### B. Azimuth and elevation angles at the observing location + +#### B. Azimuth and elevation angles at the observing location If your goal is to calculate the astrometric azimuth and zenith distance (= 90° - elevation) angles of the source at the specified observing location (without refraction correction), you can proceed from the `sky_pos` data you @@ -173,8 +185,8 @@ if you want, e.g.: zd -= refract_astro(&obs.on_surf, NOVAS_STANDARD_ATMOSPHERE, zd); ``` - -### Calculating positions for a Solar-system source + +## Calculating positions for a Solar-system source Solar-system sources work similarly to the above with a few important differences. @@ -230,45 +242,3 @@ E.g.: } ``` -### Reduced accuracy shortcuts - -When one does not need positions at the microarcsecond level, some shortcuts can be made to the recipe above: - - - You can use TT and TDB timescales interchangeably in the present era unless you require the utmost precision. - - You can use `NOVAS_REDUCED_ACCURACY` instead of `NOVAS_FULL_ACCURACY` for the calculations. This typically has an - effect at the milliarcsecond level only, but may be much faster to calculate. - - You can skip the J2000 to ICRS conversion and use J2000 coordinates directly as a fair approximation (at the - <= 22 mas level). - - You might skip the pole offsets dx, dy. These are tenths of arcsec, typically. - - -### Performance considerations - -Some of the calculations involved can be expensive from a computational perspective. For the most typical use case -however, NOVAS (and SuperNOVAS) has a trick up its sleeve: it caches the last result of intensive calculations so they -may be re-used if the call is made with the same environmental parameters again (such as JD time and accuracy). -Therefore, when calculating positions for a large number of sources at different times: - - - It is best to iterate over the sources in the inner loop while keeping the time fixed. - - You probably want to stick to one accuracy mode (`NOVAS_FULL_ACCURACY` or `NOVAS_REDUCED_ACCURACY`) to prevent - re-calculating the same quantities repeatedly to alternating precision. - - If super-high accuracy is not required `NOVAS_REDUCED_ACCURACY` mode offers much faster calculations, in general. - - -### Multi-threaded calculations - -A direct consequence of the caching of results in NOVAS is that calculations are generally not thread-safe as -implemented by the original NOVAS C 3.1 library. One thread may be in the process of returning cached values for one -set of input parameters while, at the same time, another thread is saving cached values for a different set of -parameters. Thus, when running calculations in more than one thread, the results returned may at times be incorrect, -or more precisely they may not correspond to the requested input parameters. - -While you should never call NOVAS C from multiple threads simultaneously, SuperNOVAS caches the results in thread -local variables (provided your compiler supports it), and is therefore safe to use in multi-threaded applications. -Just make sure that you: - - - use a compiler which supports the C11 language standard; - - or, compile with GCC >= 3.3; - - or else, set the appropriate non-standard keyword to use for declaring thread-local variables for your compiler in - `config.mk` or in your equivalent build setup. - diff --git a/Makefile b/Makefile index a974de53..d1ea9f51 100644 --- a/Makefile +++ b/Makefile @@ -116,7 +116,7 @@ lib/libnovas.so: lib/libsupernovas.so SO_LINK := $(LDFLAGS) -lm # Share librarry recipe -lib/%.so.$(SO_VERSION) : | lib +lib/%.so.$(SO_VERSION) : | lib Makefile $(CC) -o $@ $(CPPFLAGS) $(CFLAGS) $^ -shared -fPIC -Wl,-soname,$(subst lib/,,$@) $(SO_LINK) lib/libsolsys%.so.$(SO_VERSION): SO_LINK += -Llib -lsupernovas @@ -142,7 +142,7 @@ lib/libsolsys-ephem.so.$(SO_VERSION): $(SRC)/solsys-ephem.c | lib/libsupernovas. # Static library: novas.a -lib/libnovas.a: $(OBJECTS) | lib +lib/libnovas.a: $(OBJECTS) | lib Makefile ar -rc $@ $^ ranlib $@ @@ -165,10 +165,10 @@ Doxyfile.local: # Local documentation without specialized headers. The resulting HTML documents do not have # Google Search or Analytics tracking info. +.PHONY: local-dox local-dox: README-orig.md Doxyfile.local doxygen Doxyfile.local - .PHONY: help help: @echo diff --git a/README.md b/README.md index 03021c66..35d81b24 100644 --- a/README.md +++ b/README.md @@ -27,6 +27,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.1` release. ## Table of Contents @@ -118,7 +119,7 @@ provided by SuperNOVAS over the upstream NOVAS C 3.1 code: mean equinox of date coordinates. - Some remainder calculations in NOVAS C 3.1 used the result from `fmod()` unchecked, which resulted in angles outside - of the expected [0:π] range and was also the reason why `cal_date()` did not work for negative JD values. + of the expected [0:2π] range and was also the reason why `cal_date()` did not work for negative JD values. - Fixes `aberration()` returning NaN vectors if the `ve` argument is 0. It now returns the unmodified input vector appropriately instead. @@ -135,9 +136,6 @@ provided by SuperNOVAS over the upstream NOVAS C 3.1 code: position (e.g. the Sun for `solsys3.c`). The bug affects for example `grav_def()`, where it effectively results in the gravitational deflection due to the Sun being skipped. - - [__v1.1__] Radial velocity calculation to precede aberration and gravitational bending in `place()`, since the - radial velocity that is observed is in the geometric direction towards the source (unaffected by aberration), and - `rad_vel()` requires geometric directions also to account for the gravitational effects of the Sun and Earth. ----------------------------------------------------------------------------- @@ -188,10 +186,10 @@ the necessary variables in the shell prior to invoking `make`. For example: one tries to use the functions from `solsys1.c`). Note, that a `readeph()` implementation is not always necessary and you can provide a superior ephemeris reader implementation at runtime via the `set_ephem_provider()` call. - - If you want to use the CIO locator binary file for `cio_location()`, you can specify the path to the binary file - (e.g. `/usr/local/share/novas/cio_ra.bin`) on your system e.g. by setting the `CIO_LOCATOR_FILE` shell variable - prior to calling `make`. (The CIO locator file is not at all necessary for the functioning of the library, unless - you specifically require CIO positions relative to GCRS.) + - If you want to use the CIO locator binary file for `cio_location()`, you can specify the path to the CIO locator + file (e.g. `/usr/local/share/supernovas/CIO_RA.TXT`) on your system e.g. by setting the `CIO_LOCATOR_FILE` shell + variable prior to calling `make`. (The CIO locator file is not at all necessary for the functioning of the library, + unless you specifically require CIO positions relative to GCRS.) - If your compiler does not support the C11 standard and it is not GCC >=3.3, but provides some non-standard support for declaring thread-local variables, you may want to pass the keyword to use to declare variables as @@ -206,8 +204,8 @@ Now you are ready to build the library: will compile the shared (e.g. `lib/libsupernovas.so`) libraries, produce a CIO locator data file (e.g. `tools/data/cio_ra.bin`), and compile the API documentation (into `apidoc/`) using `doxygen` (if available). -Alternatively, you can build select components of the above with the `make` targets `shared`, `cio_file`, and -`local-dox` respectively. And, if unsure, you can always call `make help` to see what build targets are available. +Alternatively, you can build select components of the above with the `make` targets `shared`, and `local-dox` +respectively. And, if unsure, you can always call `make help` to see what build targets are available. After building the library you can install the above components to the desired locations on your system. For a system-wide install you may place the static or shared library into `/usr/local/lib/`, copy the CIO locator file to @@ -377,7 +375,7 @@ UT1 - UTC time difference (a.k.a. DUT1), and the current leap seconds. int dut1 = ...; ``` -Now we can set a standard UNIX time, for example, using the current time: +Now we can set the time of observation, for example, using the current UNIX time: ```c novas_timescale t_obs; // Structure that will define astrometric time @@ -600,10 +598,10 @@ before that level of accuracy is reached. 4. __Refraction__: Ground based observations are also subject to atmospheric refraction. SuperNOVAS offers the option to include approximate _optical_ refraction corrections either for a standard atmosphere or more precisely - using the weather parameters defined in the `on_surface` data structure that specifies the observer locations. - Note, that refraction at radio wavelengths is notably different from the included optical model. In any case you - may want to skip the refraction corrections offered in this library, and instead implement your own as appropriate - (or not at all). + using the weather parameters defined in the `on_surface` data structure that specifies the observer locations. + Note, that refraction at radio wavelengths is notably different from the included optical model, and a standard + radio refraction model is included as of version 1.1. In any case you may want to skip the refraction corrections + offered in this library, and instead implement your own as appropriate (or not at all). @@ -700,14 +698,15 @@ before that level of accuracy is reached. - New `make_planet()` and `make_ephem_object()` to make it simpler to configure Solar-system objects. + #### Added in v1.1 - New observing-frame based approach for calculations (`frames.c`). A `novas_frame` object uniquely defines both the - place and time of observation, with a set of pre-calculated transformations and constants. Once the frame is defined - it can be used very efficiently to calculate positions for multiple celestial objects with minimal - additional computational cost. The frames API is also more elegant and simpler than the low-level NOVAS C approach - for performing the same kind of calculations. And, frames are inherently thread-safe since post-creation their - internal state is never modified during the calculations. The following new functions were added: + place and time of observation, with a set of pre-calculated transformations and constants. Once the frame is + defined it can be used very efficiently to calculate positions for multiple celestial objects with minimum + additional computational cost. The frames API is also more elegant and more versatile than the low-level NOVAS C + approach for performing the same kind of calculations. And, frames are inherently thread-safe since post-creation + their internal state is never modified during the calculations. The following new functions were added: `novas_make_frame()`, `novas_change_observer()`, `novas_geom_posvel()`, `novas_geom_to_app()`, `novas_sky_pos()`, `novas_app_to_hor()`, `novas_app_to_geom()`, `novas_hor_to_app()`, `novas_make_transform()`, `novas_invert_transform()`, `novas_transform_vector()`, and `novas_transform_sky_pos()`. @@ -717,7 +716,11 @@ before that level of accuracy is reached. UT1, GPS, TAI, TT, TCG, TDB, or TCB), or to a UNIX time with `novas_set_unix_time()`. Once set, you can obtain an expression of that time in any timescale of choice via `novas_get_time()`, `novas_get_split_time()` or `novas_get_unix_time()`. And, you can create a new time specification by incrementing an existing one, using - `novas_increment_time()`, or measure time differences via `novas_diff_time()`. + `novas_increment_time()`, or measure time differences via `novas_diff_time()`, `novas_diff_tcg()`, or + `novas_diff_tcb()`. + + - Added `novas_planet_bundle` structure to handle planet positions and velocities more elegantly (e.g. for + gravitational deflection calculations). - `obs_posvel()` to calculate the observer position and velocity relative to the Solar System Barycenter (SSB). @@ -744,15 +747,19 @@ before that level of accuracy is reached. - Added humidity field to `on_surface` structure, e.g. for refraction calculations at radio wavelengths. The `make_on_surface()` function will set humidity to 0.0, but the user can set the field appropriately afterwards. - - New set of built-in refraction models to use with the frame-based `novas_app_to_hor()` function. The models - `novas_standard_refraction()` and `novas_optical_refraction()` implement the same refraction model as `refract()` - in NOVAS C 3.1, with `NOVAS_STANDARD_ATMOSPHERE` and `NOVAS_WEATHER_AT_LOCATION` respectively, including the - reversed direction provided by `refract_astro()`. The user may supply their own refraction models to - `novas_app_to_hor()` also, and may make used of the generic reversal function `novas_inv_refract` to calculate - refraction in the reverse direction (observer vs astrometric elevations) as needed. + - New set of built-in refraction models to use with the frame-based `novas_app_to_hor()` / `novas_hor_to_app()` + functions. The models `novas_standard_refraction()` and `novas_optical_refraction()` implement the same refraction + model as `refract()` in NOVAS C 3.1, with `NOVAS_STANDARD_ATMOSPHERE` and `NOVAS_WEATHER_AT_LOCATION` + respectively, including the reversed direction provided by `refract_astro()`. The user may supply their own custom + refraction also, and may make use of the generic reversal function `novas_inv_refract()` to calculate refraction in + the reverse direction (observer vs astrometric elevations) as needed. - Added radio refraction model `novas_radio_refraction()` based on the formulae by Berman & Rockwell 1976. + + - Added `cirs_to_tod()` and `tod_to_cirs()` functions for efficient tranformation between True of Date (TOD) and + Celestial Intermediate Reference System (CIRS), and vice versa. + - Added `make_cat_object()` function to create a NOVAS celestial `object` structure from existing `cat_entry` data. @@ -826,6 +833,9 @@ before that level of accuracy is reached. - [__v1.1__] `grav_def()` is simplified. It no longer uses the location type argument. Instead it will skip deflections due to a body, if the observer is within ~1500 km of its center. + - [__v1.1__] Radial velocity calculation to precede aberration and gravitational bending in `place()`, since the + radial velocity that is observed is in the geometric direction towards the source (unaffected by aberration), and + `rad_vel()` requires geometric directions also to account for the gravitational effects of the Sun and Earth. ----------------------------------------------------------------------------- diff --git a/build.mk b/build.mk index fa3089b8..6d4ff585 100644 --- a/build.mk +++ b/build.mk @@ -41,7 +41,7 @@ check: # Doxygen documentation (HTML and man pages) under apidocs/ .PHONY: dox -dox: README.md | apidoc +dox: README.md Doxyfile | apidoc @echo " [doxygen]" @$(DOXYGEN) diff --git a/include/eph_manager.h b/include/eph_manager.h index 88d2fbfd..388afdb9 100644 --- a/include/eph_manager.h +++ b/include/eph_manager.h @@ -4,10 +4,15 @@ * * SuperNOVAS header for managing 1997 version of JPL ephemerides specifically for solsys2.c. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications + * + * @sa solsys1.c */ #ifndef _EPHMAN_ diff --git a/include/novas.h b/include/novas.h index f6a4f597..0bca7d15 100644 --- a/include/novas.h +++ b/include/novas.h @@ -2,16 +2,18 @@ * @file * * @author G. Kaplan and A. Kovacs - * @version 1.0.1 + * @version 1.1.0 * - * SuperNOVAS astrometry softwate based on the Naval Observatory Vector Astrometry Software (NOVAS). + * SuperNOVAS astrometry software based on the Naval Observatory Vector Astrometry Software (NOVAS). * It has been modified to fix outstanding issues and to make it easier to use. * + * Based on the NOVAS C Edition, Version 3.1: * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications */ #ifndef _NOVAS_ @@ -20,6 +22,7 @@ #include // for M_PI #include // NULL #include +#include // The upstream NOVAS library had a set of include statements that really were not necessary // First, including standard libraries here meant that those libraries were included in the @@ -61,7 +64,7 @@ #define SUPERNOVAS_PATCHLEVEL 0 /// Additional release information in version, e.g. "-1", or "-rc1". -#define SUPERNOVAS_RELEASE_STRING "-devel" +#define SUPERNOVAS_RELEASE_STRING "" @@ -201,7 +204,7 @@ enum novas_debug_mode { NOVAS_DEBUG_OFF = 0, ///< Do not print errors and traces to the standard error (default). NOVAS_DEBUG_ON, ///< Print errors and traces to the standard error. - NOVAS_DEBUG_EXTRA ///< Print all errors and traces, even if they may be 'normal' behavior, to the standard error. + NOVAS_DEBUG_EXTRA ///< Print all errors and traces to the standard error, even if they may be acceptable behavior. }; /** @@ -719,6 +722,19 @@ typedef struct { } novas_matrix; +/** + * Position and velocity data for a set of major planets (which may include the Sun and the Moon also). + * + * @since 1.1 + * + * @sa enum novas_planet + */ +typedef struct { + int mask; ///< Bitwise mask (1 << planet-number) specifying wich planets have pos/vel data + double pos[NOVAS_PLANETS][3]; ///< [AU] Apparent positions of planets w.r.t. observer antedated for light-time + double vel[NOVAS_PLANETS][3]; ///< [AU/day] Apparent velocity of planets w.r.t. barycenter antedated for light-time +} novas_planet_bundle; + /** * A set of parameters that uniquely define the place and time of observation. The user may * initialize the frame with novas_make_frame(). Once the observer frame is set up, it can be @@ -763,18 +779,16 @@ typedef struct { novas_matrix precession; ///< precession matrix novas_matrix nutation; ///< nutation matrix (Lieske 1977 method) novas_matrix gcrs_to_cirs; ///< GCRS to CIRS conversion matrix - int pl_mask; ///< Bitwise mask (1 << planet-number) specifying wich planets have pos/vel data - double pl_pos[NOVAS_PLANETS][3]; ///< [AU] Apparent positions of planets w.r.t. observer antedated for light-time - double pl_vel[NOVAS_PLANETS][3]; ///< [AU/day] Apparent velocity of planets w.r.t. barycenter antedated for light-time + novas_planet_bundle planets; ///< Planet positions and velocities } novas_frame; /** * A transformation between two astronomical coordinate systems for the same observer * location and time. This allows for more elegant, generic, and efficient coordinate - * transformations than using the low-level NOVAS functions. + * transformations than the low-level NOVAS functions. * - * The transformation can be (should be) initialized via novas_make_trasform(), or else - * modified via novas_invert_transform() or novas_change_observer(). + * The transformation can be (should be) initialized via novas_make_trasform(), or via + * novas_invert_transform(). * * @since 1.1 * @@ -792,6 +806,8 @@ typedef struct { * The type of elevation value for which to calculate a refraction. * * @sa RefractionModel + * + * @since 1.1 */ enum novas_refraction_type { NOVAS_REFRACT_OBSERVED = -1, ///< Refract observed elevation value @@ -799,7 +815,7 @@ enum novas_refraction_type { }; /** - * Default set of gravitating bodies to use for deflection calculations in reduced accuracy mode + * Default set of gravitating bodies to use for deflection calculations in reduced accuracy mode. * * @sa grav_bodies_reduced_accuracy * @@ -809,7 +825,7 @@ enum novas_refraction_type { #define DEFAULT_GRAV_BODIES_REDUCED_ACCURACY ( (1 << NOVAS_SUN) | (1 << NOVAS_EARTH) ) /** - * Default set of gravitating bodies to use for deflection calculations in full accuracy mode + * Default set of gravitating bodies to use for deflection calculations in full accuracy mode. * * @sa grav_bodies_full_accuracy * @@ -819,24 +835,34 @@ enum novas_refraction_type { #define DEFAULT_GRAV_BODIES_FULL_ACCURACY ( DEFAULT_GRAV_BODIES_REDUCED_ACCURACY | (1 << NOVAS_JUPITER) | (1 << NOVAS_SATURN) ) /** - * Current set of gravitating bodies to use for deflection calculations in reduced accuracy mode + * Current set of gravitating bodies to use for deflection calculations in reduced accuracy mode. Each + * bit signifies whether a given body is to be accounted for as a gravitating body that bends light, + * such as the bit `(1 << NOVAS_JUPITER)` indicates whether or not Jupiter is considered as a deflecting + * body. You should also be sure that you provide ephemeris data for bodies that are designated for the + * deflection calculation. * * @sa grav_def() * @sa grav_planets() + * @sa DEFAULT_GRAV_BODIES_REDUCED_ACCURACY + * @sa set_ephem_provider() * * @since 1.1 - * @author Attila Kovacs */ extern int grav_bodies_reduced_accuracy; /** - * Current set of gravitating bodies to use for deflection calculations in full accuracy mode + * Current set of gravitating bodies to use for deflection calculations in full accuracy mode. Each + * bit signifies whether a given body is to be accounted for as a gravitating body that bends light, + * such as the bit `(1 << NOVAS_JUPITER)` indicates whether or not Jupiter is considered as a deflecting + * body. You should also be sure that you provide ephemeris data for bodies that are designated for the + * deflection calculation. * * @sa grav_def() * @sa grav_planets() + * @sa DEFAULT_GRAV_BODIES_FULL_ACCURACY + * @sa set_ephem_provider_hp() * * @since 1.1 - * @author Attila Kovacs */ extern int grav_bodies_full_accuracy; @@ -845,15 +871,17 @@ extern int grav_bodies_full_accuracy; * A function that returns a refraction correction for a given date/time of observation at the * given site on earth, and for a given astrometric source elevation * - * @param j_tt [day] Terrestrial Time (TT) based Julian data of observation + * @param jd_tt [day] Terrestrial Time (TT) based Julian data of observation * @param loc Pointer to structure defining the observer's location on earth, and local weather * @param type Whether the input elevation is observed or astrometric: REFRACT_OBSERVED (-1) or * REFRACT_ASTROMETRIC (0). * @param el [deg] Astrometric (unrefracted) source elevation * @return [arcsec] Estimated refraction, or NAN if there was an error (it should * also set errno to indicate the type of error). + * + * @since 1.1 */ -typedef double (*RefractionModel)(double j_tt, const on_surface *loc, enum novas_refraction_type type, double el); +typedef double (*RefractionModel)(double jd_tt, const on_surface *loc, enum novas_refraction_type type, double el); short app_star(double jd_tt, const cat_entry *star, enum novas_accuracy accuracy, double *ra, double *dec); @@ -1125,14 +1153,13 @@ double app_to_cirs_ra(double jd_tt, enum novas_accuracy accuracy, double ra); int obs_posvel(double jd_tdb, double ut1_to_tt, enum novas_accuracy accuracy, const observer *obs, const double *geo_pos, const double *geo_vel, double *pos, double *vel); -int obs_planets(double jd_tdb, enum novas_accuracy accuracy, const double *pos_obs, int *pl_mask, double pl_pos[][3], double pl_vel[][3]); +int obs_planets(double jd_tdb, enum novas_accuracy accuracy, const double *pos_obs, int pl_mask, novas_planet_bundle *planets); int grav_undef(double jd_tdb, enum novas_accuracy accuracy, const double *pos_app, const double *pos_obs, double *out); -int grav_planets(const double *pos_src, const double *pos_obs, int pl_mask, const double pl_pos[][3], const double pl_vel[][3], double *out); +int grav_planets(const double *pos_src, const double *pos_obs, const novas_planet_bundle *planets, double *out); -int grav_undo_planets(const double *pos_app, const double *pos_obs, enum novas_accuracy accuracy, int pl_mask, const double pl_pos[][3], - const double pl_vel[][3], double *out); +int grav_undo_planets(const double *pos_app, const double *pos_obs, const novas_planet_bundle *planets, double *out); int make_airborne_observer(const on_surface *location, const double *vel, observer *obs); @@ -1144,6 +1171,10 @@ int place_mod(double jd_tt, const object *source, enum novas_accuracy accuracy, int place_j2000(double jd_tt, const object *source, enum novas_accuracy accuracy, sky_pos *pos); +int cirs_to_tod(double jd_tt, enum novas_accuracy accuracy, const double *in, double *out); + +int tod_to_cirs(double jd_tt, enum novas_accuracy accuracy, const double *in, double *out); + // in timescale.c int novas_set_time(enum novas_timescale timescale, double jd, int leap, double dut1, novas_timespec *time); @@ -1160,9 +1191,9 @@ time_t novas_get_unix_time(const novas_timespec *time, long *nanos); double novas_diff_time(const novas_timespec *t1, const novas_timespec *t2); -double novas_tcb_diff(const novas_timespec *t1, const novas_timespec *t2); +double novas_diff_tcb(const novas_timespec *t1, const novas_timespec *t2); -double novas_tcg_diff(const novas_timespec *t1, const novas_timespec *t2); +double novas_diff_tcg(const novas_timespec *t1, const novas_timespec *t2); int novas_offset_time(const novas_timespec *time, double seconds, novas_timespec *out); @@ -1196,6 +1227,7 @@ int novas_transform_vector(const double *in, const novas_transform *transform, d int novas_transform_sky_pos(const sky_pos *in, const novas_transform *transform, sky_pos *out); + // in refract.c double novas_standard_refraction(double jd_tt, const on_surface *loc, enum novas_refraction_type type, double el); @@ -1240,7 +1272,7 @@ double novas_inv_refract(RefractionModel model, double jd_tt, const on_surface * #define ANGVEL NOVAS_EARTH_ANGVEL // Various locally used physical units -#define DAY 86400.0 ///< [s] seconds in a day +#define DAY 86400.0 ///< [s] seconds in a day #define DAY_HOURS 24.0 #define DEG360 360.0 #define JULIAN_YEAR_DAYS 365.25 @@ -1250,7 +1282,10 @@ double novas_inv_refract(RefractionModel model, double jd_tt, const on_surface * #define HOURANGLE (M_PI / 12.0) #define MAS (1e-3 * ASEC2RAD) -#define INV_MAX_ITER 100 +// On some older platform NAN may not be defined, so define it here if need be +#ifndef NAN +# define NAN (0.0/0.0) +#endif # ifndef THREAD_LOCAL # if __STDC_VERSION__ >= 201112L @@ -1262,12 +1297,13 @@ double novas_inv_refract(RefractionModel model, double jd_tt, const on_surface * # endif # endif + int novas_trace(const char *loc, int n, int offset); void novas_set_errno(int en, const char *from, const char *desc, ...); int novas_error(int ret, int en, const char *from, const char *desc, ...); /** - * Propagate an error (if any) with an offset. If the error is non-zero, it returns with the offset + * Propagates an error (if any) with an offset. If the error is non-zero, it returns with the offset * error value. Otherwise it keeps going as if it weren't even there... * * @param n {int} error code or the call that produces the error code @@ -1287,6 +1323,7 @@ double novas_vdot(const double *v1, const double *v2); int polar_dxdy_to_dpsideps(double jd_tt, double dx, double dy, double *dpsi, double *deps); +extern int novas_inv_max_iter; #endif /* __NOVAS_INTERNAL_API__ */ /// \endcond diff --git a/include/novascon.h b/include/novascon.h index e1cfe6ad..1d703b30 100644 --- a/include/novascon.h +++ b/include/novascon.h @@ -11,10 +11,13 @@ * conflicts with the super-simplistic naming scheme here. You are better off * without this. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications */ #ifndef _CONSTS_ diff --git a/include/nutation.h b/include/nutation.h index f3e967cf..a25a29ad 100644 --- a/include/nutation.h +++ b/include/nutation.h @@ -7,10 +7,13 @@ * between computational cost and precision. It provides support for both the IAU 2000A and * IAU 2000B series as well as a NOVAS-specific truncated low-precision version we call NU2000K. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications */ #ifndef _NUTATION_ diff --git a/include/solarsystem.h b/include/solarsystem.h index b823a54e..f6e427da 100644 --- a/include/solarsystem.h +++ b/include/solarsystem.h @@ -17,10 +17,18 @@ * Additionally, users may set their custom choice of major planet ephemeris handler at * runtime via the set_ephem_size(). * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications + * + * @sa solsys1.c + * @sa solsys2.c + * @sa solsys3.c + * @sa solsys-ephem.c */ #ifndef _SOLSYS_ diff --git a/src/cio_file.c b/src/cio_file.c index e03e5fe8..eb32ab86 100644 --- a/src/cio_file.c +++ b/src/cio_file.c @@ -1,13 +1,20 @@ -/* - Naval Observatory Vector Astrometry Software (NOVAS) - C Edition, Version 3.1 - - cio_file.c: Produces binary file of RA values for CIO - - U. S. Naval Observatory - Astronomical Applications Dept. - Washington, DC - http://www.usno.navy.mil/USNO/astronomical-applications +/** + * @author G. Kaplan and A. Kovacs + * + * SuperNOVAS tool to Produces binary data file of RA values for CIO. The resulting binary file + * is platform-dependent. As of SuperNOVAS version 1.1, one may use the ASCII file directly with + * SuperNOVAS. As such, there is no longer the need to produce the platform-dependent binary, and + * the use of the ASCII CIO locator file is now preferred for reasons of portability. + * + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications + * + * @sa set_cio_locator_file() */ #include diff --git a/src/eph_manager.c b/src/eph_manager.c index b5efc14e..eae736c8 100644 --- a/src/eph_manager.c +++ b/src/eph_manager.c @@ -8,12 +8,15 @@ * * This module exposes a lot of its own internal state variables globally. You probably should * not access them from outside this module, but they are kept ad globals to ensure compatibility - * with existing NOVAS C applications that do tinker with those values. + * with existing NOVAS C applications that might access those values. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications * * @sa solsys2.c */ diff --git a/src/frames.c b/src/frames.c index 36e05c9f..b4a154c9 100644 --- a/src/frames.c +++ b/src/frames.c @@ -5,7 +5,12 @@ * @author Attila Kovacs * @since 1.1 * - * Routines for higher-level and efficient repeat coordinate transformations. + * SuperNOVAS routines for higher-level and efficient repeat coordinate transformations using + * observer frames. Observer frames represent an observer location at a specific astronomical + * time (instant), which can be re-used again and again to calculate or transform positions of + * celestial sources in a a range of astronomical coordinate systems. + * + * @sa timescale.c */ /// \cond PRIVATE @@ -218,7 +223,9 @@ static int set_obs_posvel(novas_frame *frame) { static int frame_aberration(const novas_frame *frame, int dir, double *pos) { - double d, p, q, r; + const double pos0[3] = { pos[0], pos[1], pos[2] }; + double d; + int i; if(frame->v_obs == 0.0) return 0; @@ -227,24 +234,34 @@ static int frame_aberration(const novas_frame *frame, int dir, double *pos) { if(d == 0.0) return 0; - p = frame->beta * novas_vdot(pos, frame->obs_vel) / (d * frame->v_obs); - q = (1.0 + p / (1.0 + frame->gamma)) * d / C_AUDAY; - r = 1.0 + p; + // Iterate as necessary (for inverse only) + for(i = 0; i < novas_inv_max_iter; i++) { + const double p = frame->beta * novas_vdot(pos, frame->obs_vel) / (d * frame->v_obs); + const double q = (1.0 + p / (1.0 + frame->gamma)) * d / C_AUDAY; + const double r = 1.0 + p; - if(dir < 0) { - // Apparent to geometric - pos[0] = (r * pos[0] - q * frame->obs_vel[0]) / frame->gamma; - pos[1] = (r * pos[1] - q * frame->obs_vel[1]) / frame->gamma; - pos[2] = (r * pos[2] - q * frame->obs_vel[2]) / frame->gamma; - } - else { - // Geometric to apparent - pos[0] = (frame->gamma * pos[0] + q * frame->obs_vel[0]) / r; - pos[1] = (frame->gamma * pos[1] + q * frame->obs_vel[1]) / r; - pos[2] = (frame->gamma * pos[2] + q * frame->obs_vel[2]) / r; + if(dir < 0) { + const double pos1[3] = { pos[0], pos[1], pos[2] }; + + // Apparent to geometric + pos[0] = (r * pos0[0] - q * frame->obs_vel[0]) / frame->gamma; + pos[1] = (r * pos0[1] - q * frame->obs_vel[1]) / frame->gamma; + pos[2] = (r * pos0[2] - q * frame->obs_vel[2]) / frame->gamma; + + // Iterate, since p, q, r are defined by unaberrated position. + if(novas_vdist(pos, pos1) < 1e-13 * d) + return 0; + } + else { + // Geometric to apparent + pos[0] = (frame->gamma * pos0[0] + q * frame->obs_vel[0]) / r; + pos[1] = (frame->gamma * pos0[1] + q * frame->obs_vel[1]) / r; + pos[2] = (frame->gamma * pos0[2] + q * frame->obs_vel[2]) / r; + return 0; // No need to iterate... + } } - return 0; + return novas_error(-1, ECANCELED, "frame_aberration", "failed to converge"); } @@ -255,7 +272,10 @@ static int is_frame_initialized(const novas_frame *frame) { /** * Sets up a observing frame for a specific observer location, time of observation, and accuracy - * requirement. + * requirement. The frame is initialized using the currently configured planet ephemeris provider + * function (see set_planet_provider() and set_planet_provider_hp()), and in case of reduced + * accuracy mode, the currently configured IAU nutation model provider (see + * set_nutation_lp_provider()). * * * @param accuracy Accuracy requirement, NOVAS_FULL_ACCURACY (0) for the utmost precision or @@ -264,10 +284,11 @@ static int is_frame_initialized(const novas_frame *frame) { * @param time Time of observation * @param dx [mas] Earth orientation parameter, polar offset in x. * @param dy [mas] Earth orientation parameter, polar offset in y. - * @param[out] frame Pointer to the observing frame to set up. - * @return 0 if successful, 10--40: error is 10 + the error ephemeris(), + * @param[out] frame Pointer to the observing frame to configure. + * @return 0 if successful, + * 10--40: error is 10 + the error from ephemeris(), * 40--50: error is 40 + the error from geo_posvel(), - * 50--80: error is 30 + the error from sidereal_time(), + * 50--80: error is 50 + the error from sidereal_time(), * 80--90 error is 80 + error from cio_location(), * 90--100 error is 90 + error from cio_basis(). * or else -1 if there was an error (errno will indicate the @@ -277,6 +298,8 @@ static int is_frame_initialized(const novas_frame *frame) { * @sa novas_sky_pos() * @sa novas_geom_posvel() * @sa novas_make_transform() + * @sa set_planet_provider() + * @sa set_nutation_lp_provider() * * @since 1.1 * @author Attila Kovacs @@ -364,6 +387,7 @@ int novas_make_frame(enum novas_accuracy accuracy, const observer *obs, const no int novas_change_observer(const novas_frame *orig, const observer *obs, novas_frame *out) { static const char *fn = "novas_change_observer"; double jd_tdb; + int pl_mask; if(!orig || !obs || !out) return novas_error(-1, EINVAL, fn, "NULL parameter: orig=%p, obs=%p, out=%p", orig, obs, out); @@ -376,12 +400,13 @@ int novas_change_observer(const novas_frame *orig, const observer *obs, novas_fr out->state = FRAME_DEFAULT; out->observer = *obs; - out->pl_mask = (out->accuracy == NOVAS_FULL_ACCURACY) ? grav_bodies_full_accuracy : grav_bodies_reduced_accuracy; + + pl_mask = (out->accuracy == NOVAS_FULL_ACCURACY) ? grav_bodies_full_accuracy : grav_bodies_reduced_accuracy; prop_error(fn, set_obs_posvel(out), 0); jd_tdb = novas_get_time(&out->time, NOVAS_TDB); - prop_error(fn, obs_planets(jd_tdb, out->accuracy, out->obs_pos, &out->pl_mask, out->pl_pos, out->pl_vel), 0); + prop_error(fn, obs_planets(jd_tdb, out->accuracy, out->obs_pos, pl_mask, &out->planets), 0); out->state = FRAME_INITIALIZED; return 0; @@ -488,8 +513,20 @@ int novas_geom_posvel(const object *source, const novas_frame *frame, enum novas bary2obs(pos1, frame->obs_pos, pos1, &t_light); } else { - // Get position of body wrt observer, antedated for light-time. - prop_error(fn, light_time2(jd_tdb, source, frame->obs_pos, 0.0, frame->accuracy, pos1, vel1, &t_light), 50); + int got = 0; + + // If we readily have the requested planet data in the frame, use it. + if(source->type == NOVAS_PLANET) + if(frame->planets.mask & (1 << source->number)) { + memcpy(pos1, &frame->planets.pos[source->number][0], sizeof(pos1)); + memcpy(vel1, &frame->planets.vel[source->number][0], sizeof(vel1)); + got = 1; + } + + // Otherwise, get the position of body wrt observer, antedated for light-time. + if(!got) { + prop_error(fn, light_time2(jd_tdb, source, frame->obs_pos, 0.0, frame->accuracy, pos1, vel1, &t_light), 50); + } } if(pos) { @@ -625,7 +662,7 @@ int novas_geom_to_app(const novas_frame *frame, const double *pos, enum novas_re return novas_error(-1, EINVAL, fn, "invalid accuracy: %d", frame->accuracy); // Compute gravitational deflection and aberration. - prop_error(fn, grav_planets(pos, frame->obs_pos, frame->pl_mask, frame->pl_pos, frame->pl_vel, pos1), 0); + prop_error(fn, grav_planets(pos, frame->obs_pos, &frame->planets, pos1), 0); // Aberration correction frame_aberration(frame, GEOM_TO_APP, pos1); @@ -886,7 +923,7 @@ int novas_app_to_geom(const novas_frame *frame, enum novas_reference_system sys, frame_aberration(frame, APP_TO_GEOM, app_pos); // Undo gravitational deflection and aberration. - prop_error(fn, grav_undo_planets(app_pos, frame->obs_pos, frame->accuracy, frame->pl_mask, frame->pl_pos, frame->pl_vel, geom_icrs), 0); + prop_error(fn, grav_undo_planets(app_pos, frame->obs_pos, &frame->planets, geom_icrs), 0); return 0; } diff --git a/src/novas.c b/src/novas.c index d4703867..76804184 100644 --- a/src/novas.c +++ b/src/novas.c @@ -2,46 +2,45 @@ * @file * * @author G. Kaplan and A. Kovacs - * @version 1.0.1 + * @version 1.1.0 * - * SuperNOVAS astrometry softwate based on the Naval Observatory Vector Astrometry Software (NOVAS). + * SuperNOVAS astrometry software based on the Naval Observatory Vector Astrometry Software (NOVAS). * It has been modified to fix outstanding issues and to make it easier to use. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications */ +#include +#include +#include +#include +#include // stdarg.h before stdio.h (for older gcc...) + #if !COMPAT # include # include # include #endif -#include -#include -#include -#include -#include - /// \cond PRIVATE -#define __NOVAS_INTERNAL_API__ ///< Use definitions meant for internal use by SuperNOVAS only +#define __NOVAS_INTERNAL_API__ ///< Use definitions meant for internal use by SuperNOVAS only #include "novas.h" #define CIO_INTERP_POINTS 6 ///< Number of points to load from CIO interpolation table at once. -// On some older platform NAN may not be defined, so define it here if need be -#ifndef NAN -# define NAN (0.0/0.0) -#endif - #ifndef DEFAULT_SOLSYS /// Will use solarsystem() and solarsystem_hp() that is linked with application # define DEFAULT_SOLSYS 0 #endif +// <---------- GLOBAL VARIABLES --------------> + #if !DEFAULT_SOLSYS novas_planet_provider planet_call = (novas_planet_provider) solarsystem; novas_planet_provider_hp planet_call_hp = (novas_planet_provider_hp) solarsystem_hp; @@ -50,7 +49,8 @@ novas_planet_provider_hp planet_call_hp = (novas_planet_provider_hp) solarsystem /// \endcond /** - * Celestial pole offset ψ for high-precision applications. + * Celestial pole offset ψ for high-precision applications. It was visible to users in NOVAS C 3.1, + * hence we continue to expose it also for back compatibility. * * @sa EPS_COR * @sa cel_pole() @@ -58,19 +58,36 @@ novas_planet_provider_hp planet_call_hp = (novas_planet_provider_hp) solarsystem double PSI_COR = 0.0; /** - * Celestial pole offset ε for high-precision applications. + * Celestial pole offset ε for high-precision applications. It was visible to users in NOVAS C 3.1, + * hence we continue to expose it also for back compatibility. * * @sa PSI_COR * @sa cel_pole() */ double EPS_COR = 0.0; +/** + * Maximum number of iterations for convergent inverse calculations. Most iterative inverse functions should + * normally converge in a handful of iterations. In some pathological cases more iterations may be required. + * This variable sets an absolute maximum for the number of iterations in order to avoid runaway (zombie) + * behaviour. If inverse functions faile to converge, they will return a value indicating an error, and + * errno should be set to ECANCELED. + * + * @since 1.1 + */ +int novas_inv_max_iter = 100; + +// Defined in novas.h int grav_bodies_reduced_accuracy = DEFAULT_GRAV_BODIES_REDUCED_ACCURACY; +// Defined in novas.h int grav_bodies_full_accuracy = DEFAULT_GRAV_BODIES_FULL_ACCURACY; + +// <---------- LOCAL VARIABLES --------------> + /// Current debugging state for reporting errors and traces to stderr. -static enum novas_debug_mode novas_debug_state = 0; +static enum novas_debug_mode novas_debug_state = NOVAS_DEBUG_OFF; ///< Opened CIO locator data file, or NULL. static FILE *cio_file; @@ -349,6 +366,9 @@ int j2000_to_tod(double jd_tdb, enum novas_accuracy accuracy, const double *in, * * @sa j2000_to_tod() * @sa j2000_to_gcrs() + * @sa tod_to_gcrs() + * @sa tod_to_cirs() + * @sa tod_to_itrs() * * @since 1.0 * @author Attila Kovacs @@ -375,6 +395,7 @@ int tod_to_j2000(double jd_tdb, enum novas_accuracy accuracy, const double *in, * @param[out] out GCRS output 3-vector * @return 0 if successful, or else an error from frame_tie() * + * @sa j2000_to_tod() * @sa gcrs_to_j2000() * * @since 1.0 @@ -393,6 +414,7 @@ int j2000_to_gcrs(const double *in, double *out) { * @return 0 if successful, or else an error from frame_tie() * * @sa j2000_to_gcrs() + * @sa tod_to_j2000() * * @since 1.0 * @author Attila Kovacs @@ -444,7 +466,9 @@ static int gcrs_to_tod(double jd_tdb, enum novas_accuracy accuracy, const double * @return 0 if successful, or -1 if either of the vector arguments is NULL. * * @sa j2000_to_tod() - * @sa tod_to_gcrs() + * @sa tod_to_cirs() + * @sa tod_to_j2000() + * @sa tod_to_itrs() * * @since 1.0 * @author Attila Kovacs @@ -520,6 +544,8 @@ int gcrs_to_cirs(double jd_tdb, enum novas_accuracy accuracy, const double *in, * * @sa tod_to_gcrs() * @sa gcrs_to_cirs() + * @sa cirs_to_itrs() + * @sa cirs_to_tod() * * * @since 1.0 @@ -551,6 +577,78 @@ int cirs_to_gcrs(double jd_tdb, enum novas_accuracy accuracy, const double *in, return 0; } +/** + * Transforms a rectangular equatorial (x, y, z) vector from the Celestial Intermediate + * Reference System (CIRS) at the given epoch to the True of Date (TOD) reference + * system. + * + * @param jd_tt [day] Terrestrial Time (TT) based Julian date that defines + * the output epoch. Typically it does not require much precision, and + * Julian dates in other time measures will be unlikely to affect the + * result + * @param accuracy NOVAS_FULL_ACCURACY (0) or NOVAS_REDUCED_ACCURACY (1) + * @param in CIRS Input (x, y, z) position or velocity vector + * @param[out] out Output position or velocity 3-vector in the True of Date (TOD) frame. + * It can be the same vector as the input. + * @return 0 if successful, or -1 if either of the vector arguments is NULL + * or the accuracy is invalid, or 10 + the error from cio_location(), or + * else 20 + the error from cio_basis(). + * + * @sa tod_to_cirs() + * @sa cirs_to_app_ra() + * @sa cirs_to_gcrs() + * @sa cirs_to_itrs() + * + * + * @since 1.1 + * @author Attila Kovacs + */ +int cirs_to_tod(double jd_tt, enum novas_accuracy accuracy, const double *in, double *out) { + double ra_cio; // [h] R.A. of the CIO (from the true equinox) we'll calculate + + // Obtain the R.A. [h] of the CIO at the given date + prop_error("cirs_to_tod", cio_ra(jd_tt, NOVAS_FULL_ACCURACY, &ra_cio), 0); + + return spin(-15.0 * ra_cio, in, out); +} + + +/** + * Transforms a rectangular equatorial (x, y, z) vector from the True of Date (TOD) reference + * system to the Celestial Intermediate Reference System (CIRS) at the given epoch to the . + * + * @param jd_tt [day] Terrestrial Time (TT) based Julian date that defines + * the output epoch. Typically it does not require much precision, and + * Julian dates in other time measures will be unlikely to affect the + * result + * @param accuracy NOVAS_FULL_ACCURACY (0) or NOVAS_REDUCED_ACCURACY (1) + * @param in CIRS Input (x, y, z) position or velocity vector + * @param[out] out Output position or velocity 3-vector in the True of Date (TOD) frame. + * It can be the same vector as the input. + * @return 0 if successful, or -1 if either of the vector arguments is NULL + * or the accuracy is invalid, or 10 + the error from cio_location(), or + * else 20 + the error from cio_basis(). + * + * @sa cirs_to_tod() + * @sa app_to_cirs_ra() + * @sa tod_to_gcrs() + * @sa tod_to_j2000() + * @sa tod_to_itrs() + * + * + * @since 1.1 + * @author Attila Kovacs + */ +int tod_to_cirs(double jd_tt, enum novas_accuracy accuracy, const double *in, double *out) { + double ra_cio; // [h] R.A. of the CIO (from the true equinox) we'll calculate + + // Obtain the R.A. [h] of the CIO at the given date + prop_error("tod_to_cirs", cio_ra(jd_tt, NOVAS_FULL_ACCURACY, &ra_cio), 0); + + return spin(15.0 * ra_cio, in, out); +} + + /** * Set a custom function to use for regular precision (see NOVAS_REDUCED_ACCURACY) * ephemeris calculations instead of the default solarsystem() routine. @@ -1389,7 +1487,7 @@ short mean_star(double jd_tt, double tra, double tdec, enum novas_accuracy accur // Iteratively find ICRS coordinates that produce input apparent place // of star at date 'jd_tt'. - for(iter = INV_MAX_ITER; --iter >= 0;) { + for(iter = novas_inv_max_iter; --iter >= 0;) { double ra1, dec1; prop_error(fn, app_star(jd_tt, &star, accuracy, &ra1, &dec1), 20); @@ -3085,6 +3183,8 @@ short cel2ter(double jd_ut1_high, double jd_ut1_low, double ut1_to_tt, enum nova * @sa tod_to_itrs() * @sa itrs_to_cirs() * @sa gcrs_to_cirs() + * @sa cirs_to_gcrs() + * @sa cirs_to_tod() * * @since 1.0 * @author Attila Kovacs @@ -3133,6 +3233,9 @@ int cirs_to_itrs(double jd_tt_high, double jd_tt_low, double ut1_to_tt, enum nov * @sa cirs_to_itrs() * @sa itrs_to_tod() * @sa j2000_to_tod() + * @sa tod_to_gcrs() + * @sa tod_to_j2000() + * @sa tod_to_cirs() * * @since 1.0 * @author Attila Kovacs @@ -3949,7 +4052,7 @@ short geo_posvel(double jd_tt, double ut1_to_tt, enum novas_accuracy accuracy, c break; } - case NOVAS_AIRBORNE_OBSERVER: { + case NOVAS_AIRBORNE_OBSERVER: { // Airborne observer const double ivu = DAY / AU_KM; observer surf = *obs; int i; @@ -3966,7 +4069,7 @@ short geo_posvel(double jd_tt, double ut1_to_tt, enum novas_accuracy accuracy, c break; } - case NOVAS_SOLAR_SYSTEM_OBSERVER: { + case NOVAS_SOLAR_SYSTEM_OBSERVER: { // Observer in Solar orbit const object earth = { NOVAS_PLANET, NOVAS_EARTH, "Earth" }; const double tdb[2] = { jd_tdb, 0.0 }; int i; @@ -4190,12 +4293,8 @@ double d_light(const double *pos_src, const double *pos_body) { * @param pos_obs [AU] Position 3-vector of observer (or the geocenter), with respect to * origin at solar system barycenter, referred to ICRS axes, * components in AU. - * @param pl_mask Bitwise (1 << planet-number) mask indicating which planets to use, - * and have apparent positions and velocities defined. - * @param pl_pos [AU] Apparent geometric position of planets rel. to observer, antedated - * for light travel time to observer. - * @param pl_vel [AU/day] Apparent velocities of planets rel. to observer, antedated - * for light travel time to observer. + * @param planets Apparent planet data containing positions and velocities for the major + * gravitating bodies in the solar-system. * @param[out] out [AU] Position vector of observed object, with respect to origin at * observer (or the geocenter), referred to ICRS axes, corrected * for gravitational deflection, components in AU. It can be the same @@ -4213,7 +4312,7 @@ double d_light(const double *pos_src, const double *pos_body) { * @since 1.1 * @author Attila Kovacs */ -int grav_planets(const double *pos_src, const double *pos_obs, int pl_mask, const double pl_pos[][3], const double pl_vel[][3], double *out) { +int grav_planets(const double *pos_src, const double *pos_obs, const novas_planet_bundle *planets, double *out) { static const char *fn = "grav_planets"; static const double rmass[] = NOVAS_RMASS_INIT; @@ -4226,11 +4325,8 @@ int grav_planets(const double *pos_src, const double *pos_obs, int pl_mask, cons if(!out) return novas_error(-1, EINVAL, fn, "NULL output 3-vector"); - if(!pl_pos || !pl_vel) - return novas_error(-1, EINVAL, fn, "NULL input planet posvel: pl_pos=%p, pl_vel=%p", pl_pos, pl_vel); - - if(pl_pos == pl_vel) - return novas_error(-1, EINVAL, fn, "Indinstinct arrays: planet pos = vel"); + if(!planets) + return novas_error(-1, EINVAL, fn, "NULL input planet data"); // Initialize output vector of observed object to equal input vector. if(out != pos_src) @@ -4238,22 +4334,22 @@ int grav_planets(const double *pos_src, const double *pos_obs, int pl_mask, cons tsrc = novas_vlen(pos_src) / C_AUDAY; - for(i = 0; i < NOVAS_PLANETS; i++) { + for(i = 1; i < NOVAS_PLANETS; i++) { double lt, dlt, dpl, p1[3]; int k; - if((pl_mask & (1 << i)) == 0) + if((planets->mask & (1 << i)) == 0) continue; // Distance from observer to gravitating body - dpl = novas_vlen(&pl_pos[i][0]); + dpl = novas_vlen(&planets->pos[i][0]); // If observing from within ~1500 km of the gravitating body, then skip deflections by it... if(dpl < 1e-5) continue; // Calculate light time to the point where incoming geometric light ray is closest to gravitating body. - lt = d_light(pos_src, &pl_pos[i][0]); + lt = d_light(pos_src, &planets->pos[i][0]); // If gravitating body is in opposite direction from the source then use the gravitating // body position at the time the light is observed. @@ -4270,7 +4366,7 @@ int grav_planets(const double *pos_src, const double *pos_obs, int pl_mask, cons // Calculate planet position at the time it is gravitationally acting on light. for(k = 3; --k >= 0;) - p1[k] = pos_obs[k] + pl_pos[i][k] - dlt * pl_vel[i][k]; + p1[k] = pos_obs[k] + planets->pos[i][k] - dlt * planets->vel[i][k]; // Compute deflection due to gravitating body. grav_vec(out, pos_obs, p1, rmass[i], out); @@ -4294,14 +4390,8 @@ int grav_planets(const double *pos_src, const double *pos_obs, int pl_mask, cons * @param pos_obs [AU] Position 3-vector of observer (or the geocenter), with respect to * origin at solar system barycenter, referred to ICRS axes, * components in AU. - * @param accuracy NOVAS_FULL_ACCURACY (0) for sub-μas accuracy or NOVAS_REDUCED_ACCURACY (1) - * for sub-mas accuracy. - * @param pl_mask Bitwise (1 << planet-number) mask indicating which planets to use, - * and have apparent positions and velocities defined. - * @param pl_pos [AU] Apparent geometric position of planets rel. to observer, antedated - * for light travel time to observer. - * @param pl_vel [AU/day] Apparent velocities of planets rel. to observer, antedated - * for light travel time to observer. + * @param planets Apparent planet data containing positions and velocities for the major + * gravitating bodies in the solar-system. * @param[out] out [AU] Nominal position vector of observed object, with respect to origin at * observer (or the geocenter), referred to ICRS axes, without gravitational * deflection, components in AU. It can be the same vector as the input, but not @@ -4315,23 +4405,19 @@ int grav_planets(const double *pos_src, const double *pos_obs, int pl_mask, cons * @since 1.1 * @author Attila Kovacs */ -int grav_undo_planets(const double *pos_app, const double *pos_obs, enum novas_accuracy accuracy, int pl_mask, const double pl_pos[][3], - const double pl_vel[][3], double *out) { +int grav_undo_planets(const double *pos_app, const double *pos_obs, const novas_planet_bundle *planets, double *out) { static const char *fn = "grav_undo_planets"; + const double tol = 1e-13; double pos_def[3] = { }, pos0[3] = { }; double l; - double tol = accuracy == NOVAS_FULL_ACCURACY ? 1e-13 : 1e-10; int i; if(!pos_app || !pos_obs) return novas_error(-1, EINVAL, fn, "NULL input 3-vector: pos_app=%p, pos_obs=%p", pos_app, pos_obs); - if(!pl_pos || !pl_vel) - return novas_error(-1, EINVAL, fn, "NULL input planet posvel: pl_pos=%p, pl_vel=%p", pl_pos, pl_vel); - - if(pl_pos == pl_vel) - return novas_error(-1, EINVAL, fn, "Indinstinct arrays: planet pos = vel"); + if(!planets) + return novas_error(-1, EINVAL, fn, "NULL input planet data"); if(!out) return novas_error(-1, EINVAL, fn, "NULL output 3-vector: out=%p", out); @@ -4345,10 +4431,10 @@ int grav_undo_planets(const double *pos_app, const double *pos_obs, enum novas_a memcpy(pos0, pos_app, sizeof(pos0)); - for(i = 0; i < INV_MAX_ITER; i++) { + for(i = 0; i < novas_inv_max_iter; i++) { int j; - prop_error(fn, grav_planets(pos0, pos_obs, pl_mask, pl_pos, pl_vel, pos_def), 0); + prop_error(fn, grav_planets(pos0, pos_obs, planets, pos_def), 0); if(novas_vdist(pos_def, pos_app) / l < tol) { memcpy(out, pos0, sizeof(pos0)); @@ -4359,8 +4445,7 @@ int grav_undo_planets(const double *pos_app, const double *pos_obs, enum novas_a pos0[j] -= pos_def[j] - pos_app[j]; } - novas_set_errno(ECANCELED, fn, "failed to converge"); - return -1; + return novas_error(-1, ECANCELED, fn, "failed to converge"); } /** @@ -4379,17 +4464,16 @@ int grav_undo_planets(const double *pos_app, const double *pos_obs, enum novas_a * @param pos_obs [AU] Position 3-vector of observer (or the geocenter), with respect to * origin at solar system barycenter, referred to ICRS axes, * components in AU. - * @param[in,out] pl_mask Bitwise (1 << planet-number) mask indicating which planets to use - * (input), and which of these have apparent positions and velocities - * returned (output). - * @param[out] pl_pos [AU] Apparent geometric position of planets rel. to observer, antedated - * for light travel time to observer. - * @param[out] pl_vel [AU/day] Apparent velocities of planets rel. to observer, antedated - * for light travel time to observer. + * @param pl_mask Bitwise `(1 << planet-number)` mask indicating which planets to request + * data for. See enum novas_planet for the enumeration of planet numbers. + * @param[out] planets Pointer to apparent planet data to populate. + * have positions and velocities calculated successfully. See enum + * novas_planet for the enumeration of planet numbers. * @return 0 if successful, -1 if any of the pointer arguments is NULL * or if the output vector is the same as pos_obs, or the error from * ephemeris(). * + * @sa enum novas_planet * @sa grav_planets() * @sa grav_undo_planets() * @sa set_planet_provider() @@ -4398,30 +4482,23 @@ int grav_undo_planets(const double *pos_app, const double *pos_obs, enum novas_a * @since 1.1 * @author Attila Kovacs */ -int obs_planets(double jd_tdb, enum novas_accuracy accuracy, const double *pos_obs, int *pl_mask, double pl_pos[][3], double pl_vel[][3]) { +int obs_planets(double jd_tdb, enum novas_accuracy accuracy, const double *pos_obs, int pl_mask, novas_planet_bundle *planets) { static const char *fn = "obs_planets"; static object body[NOVAS_PLANETS]; static int initialized; enum novas_debug_mode dbmode = novas_get_debug_mode(); - int i, mask, error = 0; + int i, error = 0; - if(!pl_mask) - return novas_error(-1, EINVAL, fn, "NULL planet mask parameter"); + if(!planets) + return novas_error(-1, EINVAL, fn, "NULL planet data"); - mask = *pl_mask; - *pl_mask = 0; + planets->mask = 0; if(!pos_obs) return novas_error(-1, EINVAL, fn, "NULL observer position parameter"); - if(!pl_pos || !pl_vel) - return novas_error(-1, EINVAL, fn, "NULL planet pos/vel array parameter(s): pl_pos=%p, pl_vel=%p", pl_pos, pl_vel); - - if(pl_pos == pl_vel) - return novas_error(-1, EINVAL, fn, "Indinstinct arrays: planet pos = vel"); - // Set up the structures of type 'object' containing the body information. if(!initialized) { for(i = 0; i < NOVAS_PLANETS; i++) @@ -4439,30 +4516,30 @@ int obs_planets(double jd_tdb, enum novas_accuracy accuracy, const double *pos_o double tl; int stat; - if((mask & bit) == 0) + if((pl_mask & bit) == 0) continue; // Calculate positions and velocities antedated for light time. - stat = light_time2(jd_tdb, &body[i], pos_obs, 0.0, accuracy, &pl_pos[i][0], &pl_vel[i][0], &tl); + stat = light_time2(jd_tdb, &body[i], pos_obs, 0.0, accuracy, &planets->pos[i][0], &planets->vel[i][0], &tl); if(stat) { if(!error) error = stat > 10 ? stat - 10 : -1; continue; } - *pl_mask |= bit; + planets->mask |= bit; } // Re-enable debug traces novas_debug(dbmode); // If could not calculate deflection due to Sun, return with error - if((*pl_mask & (1 << NOVAS_SUN)) == 0) + if((planets->mask & (1 << NOVAS_SUN)) == 0) prop_error("grav_init_planet:sun", error, 0); // If could not get positions for another gravitating body then // return error only if in extra debug mode... - if(*pl_mask != mask && novas_get_debug_mode() == NOVAS_DEBUG_EXTRA) + if(planets->mask != pl_mask && novas_get_debug_mode() == NOVAS_DEBUG_EXTRA) prop_error(fn, error, 0); return 0; @@ -4523,14 +4600,14 @@ short grav_def(double jd_tdb, enum novas_observer_place unused, enum novas_accur double *out) { static const char *fn = "grav_def"; - double pl_pos[NOVAS_PLANETS][3] = { { } }, pl_vel[NOVAS_PLANETS][3] = { { } }; + novas_planet_bundle planets = {}; int pl_mask = (accuracy == NOVAS_FULL_ACCURACY) ? grav_bodies_full_accuracy : grav_bodies_reduced_accuracy; if(!pos_src || !out) return novas_error(-1, EINVAL, fn, "NULL source position 3-vector: pos_src=%p, out=%p", pos_src, out); - prop_error(fn, obs_planets(jd_tdb, accuracy, pos_obs, &pl_mask, pl_pos, pl_vel), 0); - prop_error(fn, grav_planets(pos_src, pos_obs, pl_mask, pl_pos, pl_vel, out), 0); + prop_error(fn, obs_planets(jd_tdb, accuracy, pos_obs, pl_mask, &planets), 0); + prop_error(fn, grav_planets(pos_src, pos_obs, &planets, out), 0); return 0; } @@ -4581,14 +4658,14 @@ short grav_def(double jd_tdb, enum novas_observer_place unused, enum novas_accur int grav_undef(double jd_tdb, enum novas_accuracy accuracy, const double *pos_app, const double *pos_obs, double *out) { static const char *fn = "grav_undef"; - double pl_pos[NOVAS_PLANETS][3], pl_vel[NOVAS_PLANETS][3]; + novas_planet_bundle planets = {}; int pl_mask = (accuracy == NOVAS_FULL_ACCURACY) ? grav_bodies_full_accuracy : grav_bodies_reduced_accuracy; if(!pos_app || !out) return novas_error(-1, EINVAL, fn, "NULL source position 3-vector: pos_app=%p, out=%p", pos_app, out); - prop_error(fn, obs_planets(jd_tdb, accuracy, pos_obs, &pl_mask, pl_pos, pl_vel), 0); - prop_error(fn, grav_undo_planets(pos_app, pos_obs, accuracy, pl_mask, pl_pos, pl_vel, out), 0); + prop_error(fn, obs_planets(jd_tdb, accuracy, pos_obs, pl_mask, &planets), 0); + prop_error(fn, grav_undo_planets(pos_app, pos_obs, &planets, out), 0); return 0; } @@ -5288,7 +5365,7 @@ double planet_lon(double t, enum novas_planet planet) { lon = 5.311886286677 + 3.813303563778 * t; break; default: - errno = EINVAL; + novas_set_errno(EINVAL, "planet_lon", "invalid planet number: %d", planet); return NAN; } @@ -5688,9 +5765,11 @@ short cio_ra(double jd_tt, enum novas_accuracy accuracy, double *ra_cio) { * @param jd_tt [day] Terrestrial Time (TT) based Julian date * @param accuracy NOVAS_FULL_ACCURACY (0) or NOVAS_REDUCED_ACCURACY (1) * @param ra [h] The CIRS right ascension coordinate, measured from the CIO. - * @return [h] the apparent R.A. coordinate measured from the true equinox of date [0:24]. + * @return [h] the apparent R.A. coordinate measured from the true equinox of date [0:24], + * or NAN if the accuracy is invalid, or if there wan an error from cio_ra(). * - * @see app_to_cirs_ra() + * @sa app_to_cirs_ra() + * @sa cirs_to_tod() * * @since 1.0.1 * @author Attila Kovacs @@ -5699,7 +5778,11 @@ double cirs_to_app_ra(double jd_tt, enum novas_accuracy accuracy, double ra) { double ra_cio; // [h] R.A. of the CIO (from the true equinox) we'll calculate // Obtain the R.A. [h] of the CIO at the given date - cio_ra(jd_tt, NOVAS_FULL_ACCURACY, &ra_cio); + int stat = cio_ra(jd_tt, accuracy, &ra_cio); + if(stat) { + novas_trace("cirs_to_app_ra", stat, 0); + return NAN; + } // Convert CIRS R.A. to true apparent R.A., keeping the result in the [0:24] h range ra = remainder(ra + ra_cio, 24.0); @@ -5716,9 +5799,11 @@ double cirs_to_app_ra(double jd_tt, enum novas_accuracy accuracy, double ra) { * @param jd_tt [day] Terrestrial Time (TT) based Julian date * @param accuracy NOVAS_FULL_ACCURACY (0) or NOVAS_REDUCED_ACCURACY (1) * @param ra [h] the apparent R.A. coordinate measured from the true equinox of date. - * @return [h] The CIRS right ascension coordinate, measured from the CIO [0:24]. + * @return [h] The CIRS right ascension coordinate, measured from the CIO [0:24], + * or NAN if the accuracy is invalid, or if there wan an error from cio_ra(). * - * @see cirs_to_app_ra() + * @sa cirs_to_app_ra() + * @sa tod_to_cirs() * * @since 1.0.1 * @author Attila Kovacs @@ -5727,7 +5812,11 @@ double app_to_cirs_ra(double jd_tt, enum novas_accuracy accuracy, double ra) { double ra_cio; // [h] R.A. of the CIO (from the true equinox) we'll calculate // Obtain the R.A. [h] of the CIO at the given date - cio_ra(jd_tt, NOVAS_FULL_ACCURACY, &ra_cio); + int stat = cio_ra(jd_tt, accuracy, &ra_cio); + if(stat) { + novas_trace("app_to_cirs_ra", stat, 0); + return NAN; + } // Convert CIRS R.A. to true apparent R.A., keeping the result in the [0:24] h range ra = remainder(ra - ra_cio, 24.0); @@ -5739,8 +5828,10 @@ double app_to_cirs_ra(double jd_tt, enum novas_accuracy accuracy, double ra) { /** * Sets the CIO interpolaton data file to use to interpolate CIO locations vs the GCRS. - * The necessary binary data file may be obtained via the cio_file.c utility - * provided in this distribution under tools/. + * You can specify either the original `CIO_RA.TXT` file included in the distribution + * (preferred since v1.1), or else a platform-specific binary data file compiled from it + * (the old way), which may be obtained via the cio_file.c utility provided + * in this distribution under tools/. * * @param filename Path (preferably absolute path) to binary data file generated * by the cio_file.c utility from the CIO_RA.TXT @@ -6309,20 +6400,26 @@ short ephemeris(const double *jd_tdb, const object *body, enum novas_origin orig case NOVAS_EPHEM_OBJECT: { enum novas_origin eph_origin = NOVAS_HELIOCENTER; - error = -1; if(readeph2_call) { // If there is a newstyle epehemeris access routine set, we will prefer it. error = readeph2_call(body->name, body->number, jd_tdb[0], jd_tdb[1], &eph_origin, posvel, &posvel[3]); } -# ifdef DEFAULT_READEPH else { +# ifdef DEFAULT_READEPH // Use whatever readeph() was compiled or the equivalent user-defined call double *res = readeph(body->number, body->name, jd_tdb[0] + jd_tdb[1], &error); - if(res == NULL) return 3; - memcpy(posvel, res, sizeof(posvel)); - free(res); + if(res == NULL) { + error = 3; + errno = ENOSYS; + } + else { + memcpy(posvel, res, sizeof(posvel)); + free(res); + } +# else + return novas_error(-1, errno, "ephemeris:ephem_object", "No ephemeris provider was defined. Call set_ephem_provider() prior."); +# endif } -# endif prop_error("ephemeris:ephem_object", error, 20); @@ -6704,7 +6801,7 @@ double refract_astro(const on_surface *location, enum novas_refraction_model opt double refr = 0.0; int i; - for(i = 0; i < INV_MAX_ITER; i++) { + for(i = 0; i < novas_inv_max_iter; i++) { double zd_obs = zd_astro - refr; refr = refract(location, option, zd_obs); if(fabs(refr - (zd_astro - zd_obs)) < 3.0e-5) diff --git a/src/novascon.c b/src/novascon.c index 7150b5ab..9f881629 100644 --- a/src/novascon.c +++ b/src/novascon.c @@ -3,17 +3,23 @@ * * @author G. Kaplan and A. Kovacs * - * SuperNOVAS implementation for numerical constants that were used internally in novas.c. In SuperNOVAS it is no longer needed - * by novas.c, and you probably don't want to use it either with your application code. + * SuperNOVAS implementation for numerical constants that were used internally in novas.c. In + * SuperNOVAS it is no longer needed by novas.c, and you probably don't want to use it either + * with your application code. * - * @deprecated Use your own version for the selection of the constant you need, expressed in whatever units your application - * desires. We should not force you to adopt the internally used convention of NOVAS, not to mention the high - * chance of namespace conflicts with the super-simplistic naming scheme here. You are better off without this. + * @deprecated Use your own version for the selection of the constant you need, expressed in + * whatever units your application desires. We should not force you to adopt the + * internally used convention of NOVAS, not to mention the high chance of namespace + * conflicts with the super-simplistic naming scheme here. You are better off + * without this. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications */ // We shall never use these in the internal API. We define precompiler constants instead... diff --git a/src/nutation.c b/src/nutation.c index 4fb8376f..1cba7214 100644 --- a/src/nutation.c +++ b/src/nutation.c @@ -8,10 +8,13 @@ * IAU 2000A and IAU2000B series as well as a NOVAS-specific truncated low-precision * version we call NU2000K. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications */ #include diff --git a/src/readeph0.c b/src/readeph0.c index 73fd1c61..d3bb6aee 100644 --- a/src/readeph0.c +++ b/src/readeph0.c @@ -1,15 +1,18 @@ -/* - Naval Observatory Vector Astrometry Software (NOVAS) - C Edition, Version 3.1 - - readeph0.c: Dummy readeph for use when minor planet ephermeris is unavailable - - U. S. Naval Observatory - Astronomical Applications Dept. - Washington, DC - http://www.usno.navy.mil/USNO/astronomical-applications +/** + * @author G. Kaplan and A. Kovacs + * + * Dummy readeph() implementation for SuperNOVAS for use when minor planet ephermeris is + * unavailable. + * + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications + * */ - #include #include @@ -30,6 +33,7 @@ double * readeph_dummy(int mp, const char *name, double jd_tdb, int *error) { if(isnanf(jd_tdb)) { set_error(-1, EINVAL, fn, "NaN jd_tdb"); + *error = -1; return NULL; } diff --git a/src/refract.c b/src/refract.c index 75cf83ab..9aae5804 100644 --- a/src/refract.c +++ b/src/refract.c @@ -4,7 +4,11 @@ * @date Created on Jun 27, 2024 * @author Attila Kovacs * - * A collection of refraction models to used with novas_app_to_hor() + * A collection of refraction models and utilities to use with novas_app_to_hor() or + * novas_hor_to_app(). + * + * @sa novas_app_to_hor() + * @sa novas_hor_to_app(). */ /// \cond PRIVATE @@ -19,7 +23,6 @@ #include #include "novas.h" - static double novas_refraction(enum novas_refraction_model model, const on_surface *loc, enum novas_refraction_type type, double el) { if(!loc) { novas_error(-1, EINVAL, "novas_refraction", "NULL on surface observer location"); @@ -60,7 +63,7 @@ double novas_inv_refract(RefractionModel model, double jd_tt, const on_surface * const int dir = (type == NOVAS_REFRACT_OBSERVED ? 1 : -1); int i; - for(i = 0; i < INV_MAX_ITER; i++) { + for(i = 0; i < novas_inv_max_iter; i++) { double el1 = el0 + dir * refr; refr = model(jd_tt, loc, type, el1); @@ -134,7 +137,9 @@ double novas_optical_refraction(double jd_tt, const on_surface *loc, enum novas_ * * * @param jd_tt [day] Terrestrial Time (TT) based Julian data of observation - * @param loc Pointer to structure defining the observer's location on earth, and local weather + * @param loc Pointer to structure defining the observer's location on earth, and local weather. + * Make sure all weather values, including humidity (added in v1.1), are fully + * populated. * @param type Whether the input elevation is observed or astrometric: REFRACT_OBSERVED (-1) or * REFRACT_ASTROMETRIC (0). * @param el [deg] source elevation of the specified type. diff --git a/src/solsys-ephem.c b/src/solsys-ephem.c index 06eeca19..7ed951d4 100644 --- a/src/solsys-ephem.c +++ b/src/solsys-ephem.c @@ -1,8 +1,8 @@ /** * @file * - * SuperNOVAS major planet ephemeris handler via the same generic ephemeris reader that is configured by set_ephem_provider() prior - * to calling this routine. + * SuperNOVAS major planet ephemeris handler via the same generic ephemeris reader that is + * configured by set_ephem_provider() prior to calling this routine. * * @date Created on Jan 29, 2024 * @author Attila Kovacs diff --git a/src/solsys1.c b/src/solsys1.c index d7d4f176..4a0a9ab3 100644 --- a/src/solsys1.c +++ b/src/solsys1.c @@ -3,15 +3,19 @@ * * @author G. Kaplan and A. Kovacs * - * SuperNOVAS major planet ephemeris lookup implementation using JPL 1997 ephemeris data, to be used together - * with eph_manager.c. A more generic solution is to implement a novas_ephem_provider (e.g. relying on the current - * version of the CSPICE library) and set it as the default ephemeris handler via set_ephem_provider(), and then - * use solsys-ephem.c instead to use the same implementation for major planets. + * SuperNOVAS major planet ephemeris lookup implementation using JPL 1997 ephemeris data, to be + * used together with eph_manager.c. A more generic solution is to implement a + * novas_ephem_provider (e.g. relying on the current version of the CSPICE library) and set it as + * the default ephemeris handler via set_ephem_provider(), and then use solsys-ephem.c instead to + * use the same implementation for major planets. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications * * @sa solsys-ephem.c * @sa solsys1.c @@ -179,10 +183,9 @@ short solarsystem(double jd_tdb, short body, short origin, double *position, dou } short solarsystem_hp(const double jd_tdb[2], short body, short origin, double *position, double *velocity) { - if(!jd_tdb) { - errno = EINVAL; - return -1; - } + if(!jd_tdb) + return novas_error(-1, EINVAL, "solarsystem_hp", "NULL jd_tdb 2-component input array"); + prop_error("solarsystem_hp", planet_eph_manager_hp(jd_tdb, body, origin, position, velocity), 0); return 0; } diff --git a/src/solsys2.c b/src/solsys2.c index de7d5225..3eea275c 100644 --- a/src/solsys2.c +++ b/src/solsys2.c @@ -13,10 +13,13 @@ * e.g. to the SPICE library, via novas_ephem_provider instead, which you can then activate * at runtime with set_planet_provider(). * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications * * @sa solsys-ephem.c * @sa solsys1.c diff --git a/src/solsys3.c b/src/solsys3.c index b2d7396d..6c3eddcb 100644 --- a/src/solsys3.c +++ b/src/solsys3.c @@ -6,10 +6,13 @@ * SuperNOVAS plane calculator functions for the Earth and Sun only, with an orbital model * based on the DE405 ephemerides by JPL. * - * Based on the NOVAS C Edition, Version 3.1, U. S. Naval Observatory - * Astronomical Applications Dept. - * Washington, DC - * http://www.usno.navy.mil/USNO/astronomical-applications + * Based on the NOVAS C Edition, Version 3.1: + * + * U. S. Naval Observatory
+ * Astronomical Applications Dept.
+ * Washington, DC
+ * + * http://www.usno.navy.mil/USNO/astronomical-applications * * @sa solsys-ephem.c */ diff --git a/src/timescale.c b/src/timescale.c index 1945a8f6..e609a55f 100644 --- a/src/timescale.c +++ b/src/timescale.c @@ -5,8 +5,10 @@ * @author Attila Kovacs * @since 1.1 * - * A set of routines to make handling of astronomical timescale and conversions between them - * easier. + * A set of SuperNOVAS routines to make handling of astronomical timescales and conversions + * among them easier. + * + * @sa frames.c */ /// \cond PRIVATE @@ -44,9 +46,9 @@ /** * Sets an astronomical time to the fractional Julian Date value, defined in the specified - * timescale. The time set this way is accurate to a few μs (microsecond) due to the inherent + * timescale. The time set this way is accurate to a few μs (microseconds) due to the inherent * precision of the double-precision argument. For higher precision applications you may use - * `novas_set_split_time()` instead, which has an inherent accuracy at the picoseconds level. + * `novas_set_split_time()` instead, which has an inherent accuracy at the picosecond level. * * @param timescale The astronomical time scale in which the Julian Date is given * @param jd [day] Julian day value in the specified timescale @@ -71,8 +73,8 @@ int novas_set_time(enum novas_timescale timescale, double jd, int leap, double d /** * Sets an astronomical time to the split Julian Date value, defined in the specified timescale. * The split into the integer and fractional parts can be done in any convenient way. The highest - * precision is reached if the fractional part is on the order of ≤= 1 day. In that case, the - * time may be specified to picosecond accuracy, if needed. + * precision is reached if the fractional part is ≤ 1 day. In that case, the time may be + * specified to picosecond accuracy, if needed. * * The accuracy of Barycentric Time measures (TDB and TCB) relative to other time measures is * limited by the precision of `tbd2tt()` implementation, to around 10 μs. @@ -203,7 +205,7 @@ int novas_offset_time(const novas_timespec *time, double seconds, novas_timespec * double-precision result. For higher precision applications you may use `novas_get_split_time()` * instead, which has an inherent accuracy at the picosecond level. * - * @param time Pointer to the astronimical time specification data structure. + * @param time Pointer to the astronomical time specification data structure. * @param timescale The astronomical time scale in which the returned Julian Date is to be * provided * @return [day] The Julian date in the requested timescale. @@ -222,9 +224,9 @@ double novas_get_time(const novas_timespec *time, enum novas_timescale timescale /** * Returns the fractional Julian date of an astronomical time in the specified timescale, as an - * integer and fractional part. The two-component split of the time allows for abolute precisions + * integer and fractional part. The two-component split of the time allows for absolute precisions * at the picosecond level, as opposed to `novas_set_time()`, whose precision is limited to a - * ew microseconds. Ultimately, the accutacy of the time will depend on how it was set. + * few microseconds typically. * * The accuracy of Barycentric Time measures (TDB and TCB) relative to other time measures is * limited by the precision of the `tbd2tt()` implemenation, to around 10 μs. @@ -242,7 +244,7 @@ double novas_get_time(const novas_timespec *time, enum novas_timescale timescale * https://gssc.esa.int/navipedia/index.php/Transformations_between_Time_Systems * * - * @param time Pointer to the astronimical time specification data structure. + * @param time Pointer to the astronomical time specification data structure. * @param timescale The astronomical time scale in which the returned Julian Date is to be * provided * @param[out] ijd [day] The integer part of the Julian date in the requested timescale. It may @@ -324,8 +326,8 @@ double novas_get_split_time(const novas_timespec *time, enum novas_timescale tim * * @sa novas_set_time() * @sa novas_offset_time() - * @sa novas_tcb_diff() - * @sa novas_tcg_diff() + * @sa novas_diff_tcb() + * @sa novas_diff_tcg() * * @since 1.1 * @author Attila Kovacs @@ -340,45 +342,47 @@ double novas_diff_time(const novas_timespec *t1, const novas_timespec *t2) { } /** - * Returns the Geocentric Coordinate Time (TCG) based time difference (t1 - t2) in days between - * two astronomical time specifications. TCG progresses slightly faster, by a relative rate about - * 1.6×10-8 higher, than time on Earth due to the lack of gravitational time - * dilation by the Earth or Sun. + * Returns the Barycentric Coordinate Time (TCB) based time difference (t1 - t2) in days between + * two astronomical time specifications. TCB progresses slightly faster than time on Earth, at a + * rate about 1.6×10-8 higher, due to the lack of gravitational time dilation by + * the Earth or Sun. * * @param t1 First time * @param t2 Second time * @return [day] Precise TCB time difference (t1-t2), or NAN if one of the inputs was * NULL (errno will be set to EINVAL) * - * @sa novas_tgc_diff() + * @sa novas_diff_tcg() * @sa novas_diff_time() * * @since 1.1 * @author Attila Kovacs */ -double novas_tcb_diff(const novas_timespec *t1, const novas_timespec *t2) { +double novas_diff_tcb(const novas_timespec *t1, const novas_timespec *t2) { return novas_diff_time(t1, t2) * (1.0 + TC_LB); } /** - * Returns the Geocentric Coordinate Time (TCG) based time difference (t1 - t2) in days between two - * astronomical time specifications. TCG progresses slightly faster, by a relative rate about - * 7×10-10 higher, than time on Earth due to the lack of gravitational time - * dilation by Earth. + * Returns the Geocentric Coordinate Time (TCG) based time difference (t1 - t2) in days between + * two astronomical time specifications. TCG progresses slightly faster than time on Earth, at a + * rate about 7×10-10 higher, due to the lack of gravitational time dilation by + * Earth. TCG is an appropriate time measure for a spacecraft that is in the proximity of the + * orbit of Earth, but far enough from Earth such that the relativistic effects of Earth's gravity + * can be ignored. * * @param t1 First time * @param t2 Second time * @return [day] Precise TCG time difference (t1-t2), or NAN if one of the inputs was * NULL (errno will be set to EINVAL) * - * @sa novas_tcb_diff() + * @sa novas_diff_tcb() * @sa novas_diff_time() * * @since 1.1 * @author Attila Kovacs */ -double novas_tcg_diff(const novas_timespec *t1, const novas_timespec *t2) { +double novas_diff_tcg(const novas_timespec *t1, const novas_timespec *t2) { return novas_diff_time(t1, t2) * (1.0 + TC_LG); } @@ -426,9 +430,9 @@ int novas_set_unix_time(time_t unix_time, long nanos, int leap, double dut1, nov /** * Returns the UNIX time for an astronomical time instant. * - * @param time The astronomical time scale in which the returned Julian Date is to be provided - * @param nanos [ns] UTC sub-second component. It may be NULL if not required. - * @return [s] The integer UNIX time + * @param time Pointer to the astronomical time specification data structure. + * @param[out] nanos [ns] UTC sub-second component. It may be NULL if not required. + * @return [s] The integer UNIX time * * @sa novas_set_unix_time() * @sa novas_get_time() @@ -446,7 +450,7 @@ time_t novas_get_unix_time(const novas_timespec *time, long *nanos) { seconds = UNIX_J2000 + (ijd - IJD_J2000) * IDAY + isod; if(nanos) { - *nanos = round(1e9 * (sod - isod)); + *nanos = floor(1e9 * (sod - isod) + 0.5); if(*nanos == E9) { seconds++; *nanos = 0; diff --git a/test/Makefile b/test/Makefile index 5471fb9e..9c92ddcb 100644 --- a/test/Makefile +++ b/test/Makefile @@ -39,6 +39,7 @@ coverage: novas.c.gcov nutation.c.gcov solsys3.c.gcov timescale.c.gcov refract.c data: mkdir data +.PHONY: bad-data bad-data: ../cio_ra.bin bad-cio-data touch bad-cio-data/empty dd if=../cio_ra.bin of=bad-cio-data/bad-1.bin bs=27 count=1 @@ -89,6 +90,21 @@ clean-data: clean: clean-test clean-cov clean-data rm -f *.o + +.PHONY: help +help: + @echo + @echo "Syntax: make [target]" + @echo + @echo "The following targets are available:" + @echo + @echo " run (default) Compiles and runs regression tests." + @echo " coverage Extracts test coverage data." + @echo " clean Removes intermediate products." + @echo " distclean Deletes all generated files." + @echo + + vpath %.c ../src vpath %.c src diff --git a/test/src/test-errors.c b/test/src/test-errors.c index 88585af5..eee5ec7e 100644 --- a/test/src/test-errors.c +++ b/test/src/test-errors.c @@ -167,6 +167,28 @@ static int test_refract() { return n; } +static int test_refract_astro() { + on_surface surf = {}; + int n = 0; + + if(check_nan("refract_astro:converge", refract_astro(&surf, NOVAS_STANDARD_ATMOSPHERE, 85.0))) { + if(check("refract_astro:converge:errno", ECANCELED, errno)) n++; + } + + return n; +} + +static int test_inv_refract() { + on_surface surf = {}; + int n = 0; + + if(check_nan("inv_refract:converge", novas_inv_refract(novas_optical_refraction, NOVAS_JD_J2000, &surf, NOVAS_REFRACT_OBSERVED, 5.0))) { + if(check("inv_refract:converge:errno", ECANCELED, errno)) n++; + } + + return n; +} + static int test_limb_angle() { double pos[3] = { 0.01 }, pn[3] = { -0.01 }, pz[3] = {}, a, b; int n = 0; @@ -273,6 +295,20 @@ static int test_cirs_to_gcrs() { return n; } +static int test_cirs_to_app_ra() { + int n = 0; + if(check_nan("cirs_to_app_ra:accuracy:-1", cirs_to_app_ra(NOVAS_JD_J2000, -1, 0.0))) n++; + if(check_nan("cirs_to_app_ra:accuracy:2", cirs_to_app_ra(NOVAS_JD_J2000, 2, 0.0))) n++; + return n; +} + +static int test_app_to_cirs_ra() { + int n = 0; + if(check_nan("app_to_cirs_ra:accuracy:-1", app_to_cirs_ra(NOVAS_JD_J2000, -1, 0.0))) n++; + if(check_nan("app_to_cirs_ra:accuracy:2", app_to_cirs_ra(NOVAS_JD_J2000, 2, 0.0))) n++; + return n; +} + static int test_set_planet_provider() { if(check("set_planet_provider", -1, set_planet_provider(NULL))) return 1; return 0; @@ -321,12 +357,16 @@ static int test_radec_planet() { } static int test_mean_star() { - double x; + double x, y; int n = 0; - if(check("mean_star:ira", -1, mean_star(0.0, 0.0, 0.0, NOVAS_FULL_ACCURACY, NULL, &x))) n++; + if(check("mean_star:ira", -1, mean_star(0.0, 0.0, 0.0, NOVAS_FULL_ACCURACY, NULL, &y))) n++; if(check("mean_star:idec", -1, mean_star(0.0, 0.0, 0.0, NOVAS_FULL_ACCURACY, &x, NULL))) n++; + if(check("mean_star:converge", 1, mean_star(NOVAS_JD_J2000, 0.0, 0.0, NOVAS_REDUCED_ACCURACY, &x, &y))) { + if(check("mean_star:converge:errno", ECANCELED, errno)) n++; + } + return n; } @@ -849,45 +889,43 @@ static int test_grav_undef() { } static int test_grav_init_planets() { - double p[3] = {2.0}, pb[NOVAS_PLANETS][3] = {{}}, vb[NOVAS_PLANETS][3] = {{}}; - int pl_mask, n = 0; + novas_planet_bundle planets = {}; + double p[3] = {2.0}; + int n = 0; - if(check("grav_init_planets:pos_obs", -1, obs_planets(NOVAS_JD_J2000, NOVAS_FULL_ACCURACY, NULL, &pl_mask, pb, vb))) n++; - if(check("grav_init_planets:pl_pos", -1, obs_planets(NOVAS_JD_J2000, NOVAS_FULL_ACCURACY, p, &pl_mask, NULL, vb))) n++; - if(check("grav_init_planets:pl_vel", -1, obs_planets(NOVAS_JD_J2000, NOVAS_FULL_ACCURACY, p, &pl_mask, pb, NULL))) n++; - if(check("grav_init_planets:pl_pos+pl_vel", -1, obs_planets(NOVAS_JD_J2000, NOVAS_FULL_ACCURACY, p, &pl_mask, NULL, NULL))) n++; - if(check("grav_init_planets:pl_pos=pl_vel", -1, obs_planets(NOVAS_JD_J2000, NOVAS_FULL_ACCURACY, p, &pl_mask, pb, pb))) n++; - if(check("grav_init_planets:pl_mask", -1, obs_planets(NOVAS_JD_J2000, NOVAS_FULL_ACCURACY, p, NULL, pb, vb))) n++; + if(check("grav_init_planets:pos_obs", -1, obs_planets(NOVAS_JD_J2000, NOVAS_FULL_ACCURACY, NULL, 0, &planets))) n++; + if(check("grav_init_planets:planets", -1, obs_planets(NOVAS_JD_J2000, NOVAS_FULL_ACCURACY, p, 0, NULL))) n++; return n; } static int test_grav_planets() { - double p[3] = {2.0}, po[3] = {0.0, 1.0}, pb[NOVAS_PLANETS][3] = {{}}, vb[NOVAS_PLANETS][3] = {{}}, out[3] = {}; - int pl_mask = 0, n = 0; + novas_planet_bundle planets = {}; + double p[3] = {2.0}, po[3] = {0.0, 1.0}, out[3] = {}; + int n = 0; - if(check("grav_planets:pos_src", -1, grav_planets(NULL, po, pl_mask, pb, vb, out))) n++; - if(check("grav_planets:pos_obs", -1, grav_planets(p, NULL, pl_mask, pb, vb, out))) n++; - if(check("grav_planets:pl_pos", -1, grav_planets(p, po, pl_mask, NULL, vb, out))) n++; - if(check("grav_planets:pl_vel", -1, grav_planets(p, po, pl_mask, pb, NULL, out))) n++; - if(check("grav_planets:pl_pos+pl_vel", -1, grav_planets(p, po, pl_mask, NULL, NULL, out))) n++; - if(check("grav_planets:pl_pos=pl_vel", -1, grav_planets(p, po, pl_mask, pb, pb, out))) n++; - if(check("grav_planets:pos_src", -1, grav_planets(p, po, pl_mask, pb, vb, NULL))) n++; + if(check("grav_planets:pos_src", -1, grav_planets(NULL, po, &planets, out))) n++; + if(check("grav_planets:pos_obs", -1, grav_planets(p, NULL, &planets, out))) n++; + if(check("grav_planets:planets", -1, grav_planets(p, po, NULL, out))) n++; + if(check("grav_planets:pos_src", -1, grav_planets(p, po, &planets, NULL))) n++; return n; } static int test_grav_undo_planets() { - double p[3] = {2.0}, po[3] = {0.0, 1.0}, pb[NOVAS_PLANETS][3] = {{}}, vb[NOVAS_PLANETS][3] = {{}}, out[3] = {}; - int pl_mask = 0, n = 0; + novas_planet_bundle planets = {}; + double p[3] = {2.0}, po[3] = {0.0, 1.0}, out[3] = {}; + int n = 0; + + if(check("grav_undo_planets:pos_app", -1, grav_undo_planets(NULL, po, &planets, out))) n++; + if(check("grav_undo_planets:pos_obs", -1, grav_undo_planets(p, NULL, &planets, out))) n++; + if(check("grav_undo_planets:planets", -1, grav_undo_planets(p, po, NULL, out))) n++; + if(check("grav_undo_planets:pos_src", -1, grav_undo_planets(p, po, &planets, NULL))) n++; - if(check("grav_undo_planets:pos_app", -1, grav_undo_planets(NULL, po, NOVAS_REDUCED_ACCURACY, pl_mask, pb, vb, out))) n++; - if(check("grav_undo_planets:pos_obs", -1, grav_undo_planets(p, NULL, NOVAS_REDUCED_ACCURACY, pl_mask, pb, vb, out))) n++; - if(check("grav_undo_planets:pl_pos", -1, grav_undo_planets(p, po, NOVAS_REDUCED_ACCURACY, pl_mask, NULL, vb, out))) n++; - if(check("grav_undo_planets:pl_vel", -1, grav_undo_planets(p, po, NOVAS_REDUCED_ACCURACY, pl_mask, pb, NULL, out))) n++; - if(check("grav_undo_planets:pl_pos+pl_vel", -1, grav_undo_planets(p, po, NOVAS_REDUCED_ACCURACY, pl_mask, NULL, NULL, out))) n++; - if(check("grav_undo_planets:pl_pos=pl_vel", -1, grav_undo_planets(p, po, NOVAS_REDUCED_ACCURACY, pl_mask, pb, pb, out))) n++; - if(check("grav_undo_planets:pos_src", -1, grav_undo_planets(p, po, NOVAS_REDUCED_ACCURACY, pl_mask, pb, vb, NULL))) n++; + planets.mask = 1 << NOVAS_SUN; + if(check("grav_undo_planets:converge", -1, grav_undo_planets(p, po, &planets, out))) { + if(check("grav_undo_planets:converge:errno", ECANCELED, errno)) n++; + } return n; } @@ -1102,6 +1140,10 @@ static int test_geom_posvel() { frame.accuracy = 2; if(check("geom_posvel:frame:accuracy:2", -1, novas_geom_posvel(&o, &frame, NOVAS_ICRS, pos, vel))) n++; + frame.accuracy = NOVAS_REDUCED_ACCURACY; + make_ephem_object("blah", 111111, &o); + if(check("geom_posvel:ephem_object", -1, novas_geom_posvel(&o, &frame, NOVAS_ICRS, pos, vel))) n++; + return n; } @@ -1153,7 +1195,9 @@ static int test_app_to_geom() { if(check("app_to_geom:frame:init", -1, novas_app_to_geom(&frame, NOVAS_ICRS, 1.0, 2.0, 10.0, pos))) n++; novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &ts, 0.0, 0.0, &frame); - if(check("app_to_geom:frame:ok", 0, novas_app_to_geom(&frame, NOVAS_ICRS, 1.0, 2.0, 10.0, pos))) n++; + if(check("app_to_geom:frame:converge", -1, novas_app_to_geom(&frame, NOVAS_ICRS, 1.0, 2.0, 10.0, pos))) { + if(check("app_to_geom:frame:converge:errno", ECANCELED, errno)) n++; + } if(check("app_to_geom:pos", -1, novas_app_to_geom(&frame, NOVAS_ICRS, 1.0, 2.0, 10.0, NULL))) n++; if(check("app_to_geom:sys:-1", -1, novas_app_to_geom(&frame, -1, 1.0, 2.0, 10.0, pos))) n++; @@ -1302,8 +1346,12 @@ static int test_inv_transform() { } int main() { + extern int novas_inv_max_iter; int n = 0; + // For testing convergence errors. + novas_inv_max_iter = 0; + if(test_make_on_surface()) n++; if(test_make_in_space()) n++; if(test_make_observer()) n++; @@ -1318,6 +1366,8 @@ int main() { if(test_transform_hip()) n++; if(test_refract()) n++; + if(test_refract_astro()) n++; + if(test_inv_refract()) n++; if(test_limb_angle()) n++; if(test_ephemeris()) n++; @@ -1326,6 +1376,8 @@ int main() { if(test_tod_to_j2000()) n++; if(test_gcrs_to_cirs()) n++; if(test_cirs_to_gcrs()) n++; + if(test_cirs_to_app_ra()) n++; + if(test_app_to_cirs_ra()) n++; if(test_set_planet_provider()) n++; if(test_set_planet_provider_hp()) n++; diff --git a/test/src/test-super.c b/test/src/test-super.c index 71404c9a..dbdadade 100644 --- a/test/src/test-super.c +++ b/test/src/test-super.c @@ -782,6 +782,23 @@ static int test_radec_planet() { } +static int test_cirs_tod() { + double pos1[3] = {}, pos2[3] = {}; + double ra0, dec0, ra1, dec1; + + if(vector2radec(pos0, &ra0, &dec0) != 0) return 0; + + if(!is_ok("cirs_tod:cirs_to_tod", cirs_to_tod(tdb, NOVAS_FULL_ACCURACY, pos0, pos1))) return 1; + + vector2radec(pos1, &ra1, &dec1); + if(!is_equal("cirs_tod:cirs_to_tod:check", cirs_to_app_ra(tdb, NOVAS_FULL_ACCURACY, ra0), ra1, 1e-10)) return 1; + + if(!is_ok("cirs_tod:tod_to_cirs", tod_to_cirs(tdb, NOVAS_FULL_ACCURACY, pos1, pos2))) return 1; + if(!is_ok("cirs_tod:tod_to_cirs:check", check_equal_pos(pos2, pos0, 1e-13 * vlen(pos0)))) return 1; + + return 0; +} + static int test_observers() { double ps[3] = { 100.0, 30.0, 10.0 }, vs[3] = { 10.0 }; int n = 0; @@ -792,6 +809,8 @@ static int test_observers() { if(test_equ_ecl()) n++; if(test_equ_gal()) n++; + if(test_cirs_tod()) n++; + make_observer_at_geocenter(&obs); n += test_source(); @@ -912,6 +931,8 @@ static int test_cirs_app_ra() { return 0; } + + static int test_set_time() { novas_timespec tt, tt1, tai, gps, TDB, tcb, tcg, utc, ut1; int leap = 32; @@ -1690,13 +1711,13 @@ static int test_diff_time() { if(!is_equal("diff_time:check", novas_diff_time(&t1, &t), 0.5, 1e-9)) return 1; if(!is_equal("diff_time:check:rev", novas_diff_time(&t, &t1), -0.5, 1e-9)) return 1; - dt = novas_tcb_diff(&t, &t1) - (1.0 + LB) * novas_diff_time(&t, &t1); + dt = novas_diff_tcb(&t, &t1) - (1.0 + LB) * novas_diff_time(&t, &t1); if(!is_ok("diff_time:check:tcb", fabs(dt) >= 1e-9)) { printf("!!! missed TCB by %.9f\n", dt); return 1; } - dt = novas_tcg_diff(&t, &t1) - (1.0 + LG) * novas_diff_time(&t, &t1); + dt = novas_diff_tcg(&t, &t1) - (1.0 + LG) * novas_diff_time(&t, &t1); if(!is_ok("diff_time:check:tcg", fabs(dt) >= 1e-9)) { printf("!!! missed TCG by %.9f\n", dt); return 1;