8000 Use the IAU 2000B nutation series for reduced (mas-level) accuracy by attipaci · Pull Request #157 · Smithsonian/SuperNOVAS · GitHub
[go: up one dir, main page]
More Web Proxy on the site http://driver.im/
Skip to content

Use the IAU 2000B nutation series for reduced (mas-level) accuracy #157

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 2 commits into from
Apr 10, 2025
Merged

Use the IAU 2000B nutation series for reduced (mas-level) accuracy #157

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5ef83fb
Use the IAU 2000B nutation series for reduced (mas-level) accuracy.
attipaci Apr 10, 2025
b9480ab
Always calculate all terms of fund_args()
attipaci Apr 10, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
11 changes: 10 additions & 1 deletion CHANGELOG.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -29,7 +29,7 @@ calculations.

- #156: `obs_posvel()` called `geo_posvel()` with TDB instead of TT. It less than a 2 ms difference, so quite
inconsequential.

### Added

- #113: New `novas_frame_lst()` convenience function to readily return the Local (apparent) Sidereal Time for a given
Expand Down Expand Up @@ -152,6 +152,13 @@ calculations.

- #156: Tighten `precession()` time checking.

- #157: Switch to `iau2000b()` as the low-precision nutation series. It is good to about 1 arcsec, which is the
promised precision in reduced accuracy, but a lot faster than the `nu2000k()` series, which in turn is barely
faster than the full `iau2000a()` nutation series.

- #157: `fund_args()` to always calculate full series. There is little performance to gain from the truncation we
used in the vicinity of J2000. No real change in the results.

- In reduced accuracy mode apply gravitational deflection for the Sun only. In prior versions, deflection corrections
were applied for Earth too. However, these are below the mas-level accuracy promised in reduced accuracy mode, and
without it, the calculations for `place()` and `novas_sky_pos()` are significantly faster.
Expand All @@ -165,6 +172,8 @@ calculations.
- Use `SIZE_OF_OBJ_NAME` and `SIZE_OF_CAT_NAME` instead of `sizeof(obj->starname)` and `sizeof(obj->catalog)`
internally for improved portability.

- Various arithmetic simplifications to give a minor performance boost.

- Updated `README.md` for v1.3 and benchmarks, including comparisons to __astropy__.

- Default `THREAD_LOCAL` definition extended to C23 `thread_local` keyword also.
Expand Down
11 changes: 7 additions & 4 deletions README.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -149,6 +149,9 @@ provided by SuperNOVAS over the upstream NOVAS C 3.1 code:

- Fixes potential string overflows and eliminates associated compiler warnings.

- Nutation series used truncated expressions for the fundamental arguments, resulting in mas-level errors in just
20--30 years away from the reference epoch of J2000.

- [__v1.1__] Fixes division by zero bug in `d_light()` 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.
Expand Down Expand Up @@ -1075,10 +1078,10 @@ aberration and gravitational deflection corrections from the observer's point of
| | full | 2708014 |
| `place()`, same time, same observer | reduced | 898609 |
| | full | 833300 |
| `novas_sky_pos()`, individual | reduced | 69568 |
| | full | 30728 |
| `place()`, individual | reduced | 61844 |
| | full | 28586 |
| `novas_sky_pos()`, individual | reduced | 173211 |
| | full | 35509 |
| `place()`, individual | reduced | 135618 |
| | full | 32293 |

For reference, we also provide the reduced accuracy benchmarks from NOVAS C 3.1.

Expand Down
54 changes: 28 additions & 26 deletions benchmark/benchmark-place.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -143,88 +143,90 @@ int main(int argc, const char *argv[]) {
clock_gettime(CLOCK_REALTIME, &unix_time);
for(i = 0; i < N; i++) calc_pos(&stars[i], &obs_frame);
clock_gettime(CLOCK_REALTIME, &end);
printf(" - reduced accuracy, same frame: %12.1f positions/sec\n",
printf(" - novas_sky_pos(), same frame, red. acc.: %12.1f positions/sec\n",
N / (end.tv_sec - unix_time.tv_sec + 1e-9 * (end.tv_nsec - unix_time.tv_nsec)));

// -------------------------------------------------------------------------
// Benchmark reduced accuracy, place(), same time
// Benchmark full accuracy, same frame
obs_frame.accuracy = NOVAS_FULL_ACCURACY;
clock_gettime(CLOCK_REALTIME, &unix_time);
for(i = 0; i < N; i++) calc_place(&stars[i], &obs_frame);
for(i = 0; i < N; i++) calc_pos(&stars[i], &obs_frame);
clock_gettime(CLOCK_REALTIME, &end);
printf(" - reduced accuracy place(), same time: %12.1f positions/sec\n",
printf(" - novas_sky_pos, same frame, full acc.: %12.1f positions/sec\n",
N / (end.tv_sec - unix_time.tv_sec + 1e-9 * (end.tv_nsec - unix_time.tv_nsec)));



// -------------------------------------------------------------------------
// Benchmark full accuracy, same frame
obs_frame.accuracy = NOVAS_FULL_ACCURACY;
// Benchmark place() reduced accuracy, same frame
clock_gettime(CLOCK_REALTIME, &unix_time);
for(i = 0; i < N; i++) calc_pos(&stars[i], &obs_frame);
for(i = 0; i < N; i++) calc_place(&stars[i], &obs_frame);
clock_gettime(CLOCK_REALTIME, &end);
printf(" - full accuracy, same frame: %12.1f positions/sec\n",
printf(" - place(), same frame, red. acc.: %12.1f positions/sec\n",
N / (end.tv_sec - unix_time.tv_sec + 1e-9 * (end.tv_nsec - unix_time.tv_nsec)));

// -------------------------------------------------------------------------
// Benchmark full accuracy, place(), same time
// Benchmark place() full accuracy, same frame
obs_frame.accuracy = NOVAS_FULL_ACCURACY;
clock_gettime(CLOCK_REALTIME, &unix_time);
for(i = 0; i < N; i++) calc_place(&stars[i], &obs_frame);
clock_gettime(CLOCK_REALTIME, &end);
printf(" - full accuracy place() same frame: %12.1f positions/sec\n",
printf(" - place(), same frame, full acc.: %12.1f positions/sec\n",
N / (end.tv_sec - unix_time.tv_sec + 1e-9 * (end.tv_nsec - unix_time.tv_nsec)));


// individual frames are expected to be significantly slower, so
// benchmark over fewer iterations
N /= 10;

// -------------------------------------------------------------------------
// Benchmark reduced accuracy different frames
// Benchmark reduced accuracy, individual fames
clock_gettime(CLOCK_REALTIME, &unix_time);
for(i = 0; i < N; i++) {
novas_set_unix_time(unix_time.tv_sec, unix_time.tv_nsec, LEAP_SECONDS, DUT1, &obs_time);
novas_make_frame(NOVAS_REDUCED_ACCURACY, &obs, &obs_time, POLAR_DX, POLAR_DY, &obs_frame);
calc_pos(&stars[i], &obs_frame);
}
clock_gettime(CLOCK_REALTIME, &end);
printf(" - reduced accuracy, own frames: %12.1f positions/sec\n",
printf(" - novas_sky_pos, individual, red. acc.: %12.1f positions/sec\n",
N / (end.tv_sec - unix_time.tv_sec + 1e-9 * (end.tv_nsec - unix_time.tv_nsec)));


// -------------------------------------------------------------------------
// Benchmark reduced accuracy, place(), different times()
obs_frame.accuracy = NOVAS_REDUCED_ACCURACY;
// Benchmark full accuracy, individual frames
clock_gettime(CLOCK_REALTIME, &unix_time);
for(i = 0; i < N; i++) {
obs_frame.time.ijd_tt += (i % 2) ? 1 : -1; // alternate dates.
calc_place(&stars[i], &obs_frame);
novas_set_unix_time(unix_time.tv_sec, unix_time.tv_nsec, LEAP_SECONDS, DUT1, &obs_time);
novas_make_frame(NOVAS_FULL_ACCURACY, &obs, &obs_time, POLAR_DX, POLAR_DY, &obs_frame);
calc_pos(&stars[i], &obs_frame);
}
clock_gettime(CLOCK_REALTIME, &end);
printf(" - reduced accuracy, place() own frames: %12.1f positions/sec\n",
printf(" - novas_sky_pos, individual, full acc.: %12.1f positions/sec\n",
N / (end.tv_sec - unix_time.tv_sec + 1e-9 * (end.tv_nsec - unix_time.tv_nsec)));



// -------------------------------------------------------------------------
// Benchmark full accuracy different frames
// Benchmark place() reduced accuracy, individual frames
obs_frame.accuracy = NOVAS_REDUCED_ACCURACY;
clock_gettime(CLOCK_REALTIME, &unix_time);
for(i = 0; i < N; i++) {
novas_set_unix_time(unix_time.tv_sec, unix_time.tv_nsec, LEAP_SECONDS, DUT1, &obs_time);
novas_make_frame(NOVAS_FULL_ACCURACY, &obs, &obs_time, POLAR_DX, POLAR_DY, &obs_frame);
calc_pos(&stars[i], &obs_frame);
obs_frame.time.ijd_tt += (i % 2) ? 1 : -1; // alternate dates.
calc_place(&stars[i], &obs_frame);
}
clock_gettime(CLOCK_REALTIME, &end);
printf(" - full accuracy, own frames: %12.1f positions/sec\n",
printf(" - place(), individual, red. acc.: %12.1f positions/sec\n",
10000 N / (end.tv_sec - unix_time.tv_sec + 1e-9 * (end.tv_nsec - unix_time.tv_nsec)));


// -------------------------------------------------------------------------
// Benchmark full accuracy, place(), different times
// Benchmark place full accuracy, individual frames
obs_frame.accuracy = NOVAS_FULL_ACCURACY;
clock_gettime(CLOCK_REALTIME, &unix_time);
for(i = 0; i < N; i++) {
obs_frame.time.ijd_tt += (i % 2) ? 1 : -1;
calc_place(&stars[i], &obs_frame);
}
clock_gettime(CLOCK_REALTIME, &end);
printf(" - full accuracy, place() own time: %12.1f positions/sec\n",
printf(" - place(), individual, full acc.: %12.1f positions/sec\n",
N / (end.tv_sec - unix_time.tv_sec + 1e-9 * (end.tv_nsec - unix_time.tv_nsec)));


Expand Down
Binary file modified resources/SuperNOVAS-benchmark.png
View file
Open in desktop
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
22 changes: 5 additions & 17 deletions src/equinox.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -329,23 +329,11 @@ int fund_args(double t, novas_delaunay_args *restrict a) {
if(!a)
return novas_error(-1, EINVAL, "fund_args", "NULL output pointer");

// higher order terms (for 0.1 uas precision) only if |t| > 0.0001
if(fabs(t) > 1e-4) {
const double t2 = t * t;
a->l = t2 * (31.8792 + t * (0.051635 + t * (-0.00024470)));
a->l1 = t2 * (-0.5532 + t * (0.000136 + t * (-0.00001149)));
a->F = t2 * (-12.7512 + t * (-0.001037 + t * (0.00000417)));
a->D = t2 * (-6.3706 + t * (0.006593 + t * (-0.00003169)));
a->Omega = t2 * (7.4722 + t * (0.007702 + t * (-0.00005939)));
}
else
memset(a, 0, sizeof(*a));

a->l += 485868.249036 + t * 1717915923.2178;
a->l1 += 1287104.793048 + t * 129596581.0481;
a->F += 335779.526232 + t * 1739527262.8478;
a->D += 1072260.703692 + t * 1602961601.2090;
a->Omega += 450160.398036 - t * 6962890.5431;
a->l = 485868.249036 + t * (1717915923.2178 + t * (31.8792 + t * (0.051635 + t * (-0.00024470))));
a->l1 = 1287104.793048 + t * (129596581.0481 + t * (-0.5532 + t * (0.000136 + t * (-0.00001149))));
a->F = 335779.526232 + t * (1739527262.8478 + t * (-12.7512 + t * (-0.001037 + t * (0.00000417))));
a->D = 1072260.703692 + t * (1602961601.2090 + t * (-6.3706 + t * (0.006593 + t * (-0.00003169))));
a->Omega = 450160.398036 + t * (-6962890.5431 + t * (7.4722 + t * (0.007702 + t * (-0.00005939))));

a->l = novas_norm_ang(a->l * ARCSEC);
a->l1 = novas_norm_ang(a->l1 * ARCSEC);
Expand Down
2 changes: 1 addition & 1 deletion src/plugin.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -34,7 +34,7 @@ static novas_ephem_provider readeph2_call = NULL;

/// Function to use for reduced-precision calculations. (The full IAU 2000A model is used
/// always for high-precision calculations)
static novas_nutation_provider nutate_lp = nu2000k;
static novas_nutation_provider nutate_lp = iau2000b;

/**
* Sets the function to use for obtaining position / velocity information for minor planets,
Expand Down
27 changes: 11 additions & 16 deletions src/timescale.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,20 +69,21 @@ int strncasecmp(const char *s1, const char *s2, size_t n);
* Computes the Terrestrial Time (TT) or Terrestrial Dynamical Time (TDT) Julian date
* corresponding to a Barycentric Dynamical Time (TDB) Julian date.
*
* Expression used in this function is a truncated form of a longer and more precise
* series given in the first reference. The result is good to about 10 microseconds.
* The simple analytical expression containing just the leading orbital term from Moyer
* 1981 is fast to calculate and goot to about 30 us. This is sufficient toobtain
* meter-level absolute positional accuracy for Solar-system sources. (Relative
* positional accuracy is
*
* @deprecated Use the less computationally intensive an more accurate tt2tdb()
* routine instead.
*
* REFERENCES:
* <ol>
* <li>Fairhead, L. & Bretagnon, P. (1990) Astron. & Astrophys. 229, 240.</li>
* <li>Kaplan, G. (2005), US Naval Observatory Circular 179.</li>
* <li><a href="https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/time.html#The%20Relationship%20between%20TT%20and%20TDB">
* https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/FORTRAN/req/time.html</a></li>
* <li>Moyer, T.D. (1981), Celestial mechanics, Volume 23, Issue 1, pp. 57-68<.li>
* <li><a href="https://gssc.esa.int/navipedia/index.php/Transformations_between_Time_Systems">
* https://gssc.esa.int/navipedia/index.php/Transformations_between_Time_Systems</a></li>
* <li>Fairhead, L. & Bretagnon, P. (1990) Astron. & Astrophys. 229, 240.</li>
* <li>Kaplan, G. (2005), US Naval Observatory Circular 179.</li>
* </ol>
*
* @param jd_tdb [day] Barycentric Dynamic Time (TDB) based Julian date
Expand All @@ -99,17 +100,11 @@ int tdb2tt(double jd_tdb, double *restrict jd_tt, double *restrict secdiff) {
static THREAD_LOCAL double d;

if(!novas_time_equals(jd_tdb, last_tdb)) {
const double t = (jd_tdb - JD_J2000) / JULIAN_CENTURY_DAYS;

// Expression given in USNO Circular 179, eq. 2.6.
// AK: This agrees poorly with the canonical expression further below...
//d = 0.001657 * sin(628.3076 * t + 6.2401) + 0.000022 * sin(575.3385 * t + 4.2970) + 0.000014 * sin(1256.6152 * t + 6.1969)
//+ 0.000005 * sin(606.9777 * t + 4.0212) + 0.000005 * sin(52.9691 * t + 0.4444) + 0.000002 * sin(21.3299 * t + 5.5431)
//+ 0.000010 * t * sin(628.3076 * t + 4.2490);

// The simpler formula with a precision of ~30 us.
const double g = 6.239996 + 630.0221385924 * t;
d = 0.001657 * sin(g + 0.01671 * sin(g));
const double t = (jd_tdb - NOVAS_JD_J2000) / JULIAN_CENTURY_DAYS;
d = 0.001657 * sin(628.3076 * t + 6.2401) + 0.000022 * sin(575.3385 * t + 4.2970) + 0.000014 * sin(1256.6152 * t + 6.1969)
+ 0.000005 * sin(606.9777 * t + 4.0212) + 0.000005 * sin(52.9691 * t + 0.4444) + 0.000002 * sin(21.3299 * t + 5.5431)
+ 0.000010 * t * sin(628.3076 * t + 4.2490);

last_tdb = jd_tdb;
}
Expand Down
60 changes: 30 additions & 30 deletions test/novasc3.1/aberration.out
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,62 +59,62 @@

10000.0 16-20 S2 O2 A0: -0.46964757761 -0.81389470906 -0.34206220986 0 0 0

-10000.0 22+20 S2 O0 A1: 0.81383625 -0.46976251 0.34204348 0 0 0
-10000.0 22+20 S2 O0 A1: 0.8138362 -0.4697625 0.3420435 0 0 0

-10000.0 22+20 S2 O1 A1: 0.81387407 -0.46967986 0.34206698 0 0 0
-10000.0 22+20 S2 O1 A1: 0.8138741 -0.4696799 0.3420670 0 0 0

-10000.0 22+20 S2 O2 A1: 0.81388607 -0.46966596 0.34205753 0 0 0
-10000.0 22+20 S2 O2 A1: 0.8138861 -0.4696660 0.3420575 0 0 0

-10000.0 16-20 S2 O0 A1: -0.46983397 -0.81380408 -0.34202186 0 0 0
-10000.0 16-20 S2 O0 A1: -0.4698340 -0.8138041 -0.3420219 0 0 0

-10000.0 16-20 S2 O1 A1: -0.46982216 -0.81381036 -0.34202315 0 0 0
-10000.0 16-20 S2 O1 A1: -0.4698222 -0.8138104 -0.3420231 0 0 0

-10000.0 16-20 S2 O2 A1: -0.46979563 -0.81382324 -0.34202894 0 0 0
-10000.0 16-20 S2 O2 A1: -0.4697956 -0.8138232 -0.3420289 0 0 0

0.0 22+20 S2 O0 A1: 0.81375973 -0.46989854 0.34203868 0 0 0
0.0 22+20 S2 O0 A1: 0.8137597 -0.4698985 0.3420387 0 0 0

0.0 22+20 S2 O1 A1: 0.81372223 -0.46995031 0.34205679 0 0 0
0.0 22+20 S2 O1 A1: 0.8137222 -0.4699503 0.3420568 0 0 0

0.0 22+20 S2 O2 A1: 0.81373305 -0.46993802 0.34204793 0 0 0
0.0 22+20 S2 O2 A1: 0.8137330 -0.4699380 0.3420479 0 0 0

0.0 16-20 S2 O0 A1: -0.46991614 -0.81376333 -0.34200593 0 0 0
0.0 16-20 S2 O0 A1: -0.4699161 -0.8137633 -0.3420059 0 0 0

0.0 16-20 S2 O1 A1: -0.46998479 -0.81372958 -0.34199192 0 0 0
0.0 16-20 S2 O1 A1: -0.4699848 -0.8137296 -0.3419919 0 0 0

0.0 16-20 S2 O2 A1: -0.46995998 -0.81374174 -0.34199708 0 0 0
0.0 16-20 S2 O2 A1: -0.4699600 -0.8137417 -0.3419971 0 0 0

10000.0 22+20 S2 O0 A1: 0.81381282 -0.46985608 0.34197069 0 0 0
10000.0 22+20 S2 O0 A1: 0.8138128 -0.4698561 0.3419707 0 0 0

10000.0 22+20 S2 O1 A1: 0.81382804 -0.46986538 0.34192168 0 0 0
10000.0 22+20 S2 O1 A1: 0.8138280 -0.4698654 0.3419217 0 0 0

10000.0 22+20 S2 O2 A1: 0.81383922 -0.46985310 0.34191196 0 0 0
10000.0 22+20 S2 O2 A1: 0.8138392 -0.4698531 0.3419120 0 0 0

10000.0 16-20 S2 O0 A1: -0.46975995 -0.81383982 -0.34203850 0 0 0
10000.0 16-20 S2 O0 A1: -0.4697599 -0.8138398 -0.3420385 0 0 0

10000.0 16-20 S2 O1 A1: -0.46967475 -0.81388122 -0.34205700 0 0 0
10000.0 16-20 S2 O1 A1: -0.4696748 -0.8138812 -0.3420570 0 0 0

10000.0 16-20 S2 O2 A1: -0.46964758 -0.81389471 -0.34206221 0 0 0
10000.0 16-20 S2 O2 A1: -0.4696476 -0.8138947 -0.3420622 0 0 0

10000.0 22+20 S2 O0 A1: 0.81381282 -0.46985608 0.34197069 0 0 0
10000.0 22+20 S2 O0 A1: 0.8138128 -0.4698561 0.3419707 0 0 0

10000.0 22+20 S2 O1 A1: 0.81382804 -0.46986538 0.34192168 0 0 0
10000.0 22+20 S2 O1 A1: 0.8138280 -0.4698654 0.3419217 0 0 0

10000.0 22+20 S2 O2 A1: 0.81383922 -0.46985310 0.34191196 0 0 0
10000.0 22+20 S2 O2 A1: 0.8138392 -0.4698531 0.3419120 0 0 0

10000.0 16-20 S2 O0 A1: -0.46975995 -0.81383982 -0.34203850 0 0 0
10000.0 16-20 S2 O0 A1: -0.4697599 -0.8138398 -0.3420385 0 0 0

10000.0 16-20 S2 O1 A1: -0.46967475 -0.81388122 -0.34205700 0 0 0
10000.0 16-20 S2 O1 A1: -0.4696748 -0.8138812 -0.3420570 0 0 0

10000.0 16-20 S2 O2 A1: -0.46964758 -0.81389471 -0.34206221 0 0 0
10000.0 16-20 S2 O2 A1: -0.4696476 -0.8138947 -0.3420622 0 0 0

10000.0 22+20 S2 O0 A1: 0.81381283 -0.46985607 0.34197069 0 0 0
10000.0 22+20 S2 O0 A1: 0.8138128 -0.4698561 0.3419707 0 0 0

10000.0 22+20 S2 O1 A1: 0.81382801 -0.46986543 0.34192170 0 0 0
10000.0 22+20 S2 O1 A1: 0.8138280 -0.4698654 0.3419217 0 0 0

10000.0 22+20 S2 O2 A1: 0.81383923 -0.46985307 0.34191196 0 0 0
10000.0 22+20 S2 O2 A1: 0.8138392 -0.4698531 0.3419120 0 0 0

10000.0 16-20 S2 O0 A1: -0.46975994 -0.81383982 -0.34203850 0 0 0
10000.0 16-20 S2 O0 A1: -0.4697599 -0.8138398 -0.3420385 0 0 0

10000.0 16-20 S2 O1 A1: -0.46967478 -0.81388122 -0.34205697 0 0 0
10000.0 16-20 S2 O1 A1: -0.4696748 -0.8138812 -0.3420570 0 0 0

10000.0 16-20 S2 O2 A1: -0.46964758 -0.81389471 -0.34206221 0 0 0
10000.0 16-20 S2 O2 A1: -0.4696476 -0.8138947 -0.3420622 0 0 0
Loading
0