This is a follow-up to work provided by:
Eric Wheeler, KJ7LNW <astronomy-git@z.ewheeler.org>
Before now, the C function Astronomy_GetCurrentTime returned
the current time from the system clock, but only with whole
second resolution. Now it supports microsecond resolution on
Linux/Unix, Mac OS, and Windows.
For unsupported platforms, a compiler error will occur
to indicate that microsecond resolution is not available.
However, it is possible to define one of the following two
preprocessor symbols to work around the compiler error:
1. ASTRONOMY_ENGINE_NO_CURRENT_TIME
Excludes the function Astronomy_CurrentTime from the build.
If your project does not need to obtain the current time,
or your hardware platform does not provide current date
and time in the first place, this is likely the better option.
2. ASTRONOMY_ENGINE_WHOLE_SECOND
If your project does need to use the current date and time
for astronomy calculations, and it can tolerate whole
second resolution, this option provides a version of
Astronomy_CurrentTime that uses a call to `time(NULL)`.
Notes:
- Added unit test to confirm at least millisecond resolution.
Because these tests have to run on GitHub Actions cloud platform,
and those systems can be heavily CPU-loaded, I want to be tolerant
of resolution and avoid false failures.
- Added detection of Mac platform.
- Added preprocessor options documented above.
- On Windows, use ULARGE_INTEGER and eliminated one integer division.
- Added comments and developer documentation.
- Converted tabs to spaces in astronomy.c, for consistent code format.
Because I plan on adding metersAboveGround as a parameter
that defaults to 0.0 in the other languages, and I want
the language implementations to be reasonably consistent,
I moved the metersAboveGround parameter to the end of
the parameter list for the C version of SearchRiseSetEx.
I realized I had to rework the RiseSetEx function so that
it accepts a height above ground level, rather than a generic
altitude angle correction, because atmospheric effects are
important both for the horizon dip angle and for the amount
of refraction at ground level.
The atmosphere calculations were interesting enough that
I made them public as a new function Astronomy_Atmosphere.
This returns the idealized temperature, pressure, and relative
density of air at the given elevation above/below sea level.
This is the first step toward calculating body rise/set times
for an observer that is significantly above the ground.
It figures out the angular correction of the horizon
using both parallax and refractive correction of a light
ray from the horizon to the observer's eye.
Enhanced the Time class to correctly calculate calendar
dates for the year range -999999 to +999999.
Made unit tests in C, C#, and Kotlin all exercise
the full year range, for February 28 and March 1 in each year,
to make sure we cover before and after each potential leap day.
With more rigorous testing, I discovered more bugs
in the C functions for converting calendar dates
to times and vice versa.
Astronomy_UtcFromTime():
When the year went before -4714, the value of the variable
`djd` went negative, causing the typecast `(long)djd` to
round toward zero instead of taking the true floor.
Changed this to `(long)floor(djd)`.
Astronomy_MakeTime():
Reworked the logic so that none of the integer divisions
involve negative values over the year range -999999..+999999.
Addressed limitations of the logic I copied from NOVAS cal_date().
Its formulas did not work for years much before -12000 due to
integer division going negative. I figured out how to make the
formulas work for plus or minus 1 million years from the present era.
We already had a function to search for the next time a body
reaches a certain hour angle. But we didn't have a function
to ask what the current hour angle of a body is.
This will resolve that problem, which will also answer
questions about true solar time: use the Sun as the body,
and add 12 to the hour angle, modulo 24.
It looks like the C unit tests can use __func__ to avoid
hardcoded function names everywhere. I know this is defined
in C++11 and C99, so it should be fine. If it passes unit
testing across all platforms I exercise in GitHub Actions,
I feel comfortable using it, since it will not affect
production code.
Now the C function Astronomy_Ecliptic returns ecliptic
coordinates in true equinox of date instead of the
J2000 mean equinox. I'm doing this because it is a
better fit for physical phenomena that ecliptic
coordinates are often used for. For example, lunar nodes,
eclipses, phase angles, and oppositions make more sense
with true latitude and longitude of date.
I will port these changes to the other languages also.
More work standardizing the nomenclature of the
orientation systems across all language documents.
Added C functions to calculate rotation matrices
for EQJ/ECT and ECT/EQJ.
Added EclipticGeoMoon as output to the temp/*_check.txt files as 'm' lines.
This ensures that all the languages calculate nearly identical values.
Optimized EclipticGeoMoon a little more by eliminating a redundant
call to mean_obliq.
Updated the C function Astronomy_EclipticGeoMoon to
correct ecliptic coordinates for nutation.
This means that the returned value is expressed in
true equinox of date instead of mean equinox of date.
This results in the moon_ecm test decreasing the max
longitude error from 24 arcseconds to 6 arcseconds.
EclipticGeoMoon is now about 40% slower, but it still
runs in about 0.4 microseconds per call.
While trying to convert ecliptic coordinates from mean
equinox of date to true equinox of date, I ran into excessive
overhead from the IAU2000B nutation model. The fact that it
uses 77 trigonometric terms made the calculations a lot slower.
https://apps.dtic.mil/sti/pdfs/AD1112517.pdf
Page 4 in the above document mentions a shorter series
“NOD version 2” that has 13 terms instead of 77 as used in IAU2000B.
I had not noticed NOD2 before, because it appears only in
the FORTRAN version of NOVAS 3.x, not the C version.
After reading the FORTRAN code, I realized NOD2 is the same
as IAU2000B, only it keeps the first 13 of 77 terms.
The terms are already arranged in descending order of
significance, so it is easy to truncate the series.
Based on this discovery, I realized I could achieve all of
the required accuracy needed for Astronomy Engine by
keeping only the first 5 terms of the nutation series.
This tremendously speeds up nutation calculations while
sacrificing only a couple of arcseconds of accuracy.
It also makes the minified JavaScript code smaller:
Before: 119500 bytes.
After: 116653 bytes.
So that's what I did here. Most of the work was updating
unit tests for accepting slightly different calculation
results.
The nutation formula change did trigger detection of a
lurking bug in the inverse_terra functions, which convert
a geocentric vector into latitude, longitude, and elevation
(i.e. an Observer object). The Newton's Method loop in
this function was not always converging, resulting in
an infinite loop. I fixed that by increasing the
convergence threshold and throwing an exception
if the loop iterates more than 10 times.
I also fixed a couple of bugs in the `demotest` scripts.
Added a performance test to measure how fast
a call to EclipticGeoMoon is. Results:
Over 3 trials, mean time = 24.064 seconds.
Number of calls = 14,609,700.
Time per call = 1.647 microseconds.
That is faster than I thought!
Also replaced conditional compilation with a runtime flag
to indicate whether a test should be skipped by the "all"
command. This allows me to run performance tests without
hacking a #define in the code, plus it ensures that even
though I don't run the performance tests every time I change
Astronomy Engine, they at least compile without errors.
This was another test where I was trying to figure out
an apparent anomaly that was actually based on my
flawed assumption about why obliquity changes over time.
I don't need the MoonLatitudes test any more.
Its purpose was to investigate what I thought was an
ecliptic latitude discrepancy. Now I understand that
changes in obliquity are caused by changes in the ecliptic
plane, not the equatorial plane.
I'm still trying to understand a discrepancy
between ecliptic latitudes calculated in two
different equinoxes: mean equinox of date and
mean J2000 equinox. I keep thinking the two
latitude values should be the same, because they
are measured against an essentially unchanging
ecliptic plane. I understand that ecliptic longitudes
should change as the Earth's equator precesses
by 47 arcseconds per century.
Added new test MoonEcliptic(), that compares
JPL Horizons calculations of the Moon's ecliptic
coordinates relative to equator of date versus
Astronomy_EclipticGeoMoon():
C MoonEcliptic: PASS: count=73050, max lat=1.825 arcsec, max lon=23.462 arcsec.
It is so weird that the maximum latitude discrepancy is less than 2 arcseconds.
Updated MoonVector test to measure maximum errors in addition to RMS errors.
Display distance errors in kilometers instead of dimensionless units.
Defined constant JD_2000 so all my tests that
convert Julian Dates from JPL Horizons can share it.
There is a discrepancy between ecliptic latitudes measured
using two different mean equinoxes: 1900 and 2000.
The discrepancy is very close to the change of obliquity
expected for a 100-year period: 47 arcseconds.
But ecliptic latitudes should remain constant as
the Earth's axis precesses!
MoonLatitudes now generates a CSV file for me to study the
ecliptic latitude discrepancy as a function of time.
I can see that it increases away from the year 2000,
and that it is periodic with the Moon's orbit.
This indicates an error in correcting for precession.
Generated new JPL Horizons test data for MoonVector that
covers the years 1900..2100 instead of 1980..2020.
Confirmed that its maximum error is still 1.16 arcsec.
That is error comparing EQJ vectors and ECL vectors.
The C unit test program ctest.c now includes a MoonLatitudes
function invoked by the name `moon_lat`.
This function demonstrates that there is a 47 arcsecond
discrepancy between ecliptic latitudes using the mean
equator of date and ecliptic latitudes using the mean
equator of J2000. This should not happen! I'm investigating.
Also, when running ctest without command-line arguments,
print a list of all available test names, so I don't have
to go open up the file ctest.c and read the code.
Added Python support for user-defined stars.
Defined new StateVector methods: Position and Velocity.
Defined division operator: Vector / float.
Bumped version number to 2.1.12.
SearchRiseSet and SearchHourAngle now work with user-defined stars.
Fixed a bug in InternalSearchAltitude where I failed to return an error
code when its call to MaxAltitudeSlope failed.
The C# version of the altitude searches now work correctly
near the Earth's poles. This is a port of the C version.
Cleaned up a little code in the C version.
Do not recurse deeper when any solution would
have a rise+set or a set+rise consecutively
within one second. This would risk confusing the
search, which has a 0.1 second tolerance for convergence.
For practical reasons, such events are dubious anyway,
being swamped by unpredictable atmospheric refraction.
Plus it's just faster to not have to recurse so deeply.
Add a safety valve that fails the search if recursion
gets deeper than the above tolerance would ever allow
to happen in the first place.
This change improves performance significantly:
C RiseSet: passed 5909 lines: time errors in minutes: rms=0.2904, max=1.1541, recur=15, altcount=156881
Added helper function AscentError for cleaner code in error cases.
Track the maximum recursion depth in FindAscent.
Count the number of times the altitude_diff function is called.
These are important performance statistics that I want to
measure, so I can try a few optimization ideas.
Current results:
don@spearmint:~/github/astronomy/generate $ ./generate source && ./ctbuild && ./ctest -v riseset
./ctbuild: C compiler = gcc
./ctbuild: Built 'ctest' program.
C RiseSet: Moon lat=-61.0 lon=103.0
C RiseSet: Moon lat=-45.0 lon=150.0
C RiseSet: Moon lat=29.0 lon=-81.0
C RiseSet: Moon lat=80.0 lon=130.0
C RiseSet: Moon lat=83.0 lon=105.0
C RiseSet: Moon lat=84.0 lon=125.0
C RiseSet: Moon lat=85.0 lon=135.0
C RiseSet: Moon lat=85.0 lon=140.0
C RiseSet: Sun lat=-90.0 lon=0.0
C RiseSet: Sun lat=-88.0 lon=45.0
C RiseSet: Sun lat=-60.0 lon=-150.0
C RiseSet: Sun lat=15.0 lon=75.0
C RiseSet: Sun lat=29.0 lon=-81.0
C RiseSet: Sun lat=60.0 lon=0.0
C RiseSet: Sun lat=89.0 lon=30.0
C RiseSet: Sun lat=90.0 lon=0.0
C RiseSet: passed 5909 lines: time errors in minutes: rms=0.2904, max=1.1541, recur=19, altcount=171438
This is a whole new algorithm that efficiently finds
all rise/set events, even near the poles.
It uses a recursive bisection search that limits
recursion depth by knowing the maximum possible
|da/dt| = change in altitude with respect to time.
