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.
I bootstrapped based on the pretty good optimizations that
codegen did for the Python version of the (now truncated)
IAU2000B nutation formula. I will do the same for the other
nutation formulas.
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.
Modified the code generator and source templates to allow
using fewer than 77 terms from the IAU2000B nutation model.
The number of terms causes calculations to be far slower than
I would like, plus most of these terms provide far less than
one arcsecond of difference in the output.
I want to experiment with truncated versions of the series.
https://apps.dtic.mil/sti/pdfs/AD1112517.pdf
Page 4 in the above document references a shorter series
"NOD version 2" that has only 13 terms instead of 77 as used in IAU2000B.
I found that code and confirmed the 13 terms are the first 13
consecutive terms from the original 77.
This commit does not change the number of terms calculated;
it only enables doing so by changing a single #define in
generate/codegen.c.
Once I change the nutation model, I'm sure there will be multiple
tweaks to unit tests to get everything working again.
Astronomy_SphereFromVector was not checking its vector
argument for having a bad status. Now it does.
Added comments that clarify exactly what nutation and precession
functions do.
The documentation for SearchRiseSet and SearchAltitude needed
clarification about refraction and the part of the body solved
for (center versus limb). The JavaScript version was especially
lacking compared to documentation for the other languages.
Also documented SearchAltitude's limitations; it does not
work at or near maximum/minimum altitude.
Mention that user-defined stars are allowed for
SearchRiseSet, SearchAltitude, and SearchHourAngle.
Fixed a couple places where the Kotlin documentation had
broken links to other functions.
I had a copy-n-paste typo in the `dec` parameters
for all of the DefineStar functions. Fixed it.
The TypeScript version of HelioState did not handle
user-defined stars. Added support there.
Because we instantly know the heliocentric
distance of a user-defined star, there is no
need to convert it into a vector and then take
the length of the vector.
All of the HelioDistance functions now return
the distance directly, as an optimization.
Also, I decided it didn't make sense to have a
default definition for user-defined stars.
If the caller doesn't define a star, it should
be treated as an invalid body.
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.
DefineStar now requires passing in the heliocentric
distance of the star expressed in light-years.
That way, I can directly support returning vectors
to a star from HelioVector, GeoVector, etc.
SearchRiseSet and SearchHourAngle now work with user-defined stars.
I'm starting to implement the ability to define
up to 8 distinct points in the sky as "stars"
that will be allowed as a `body` parameter to
some Astronomy Engine functions, to be determined.
Made sure all the altitude search functions
verify that the geographic latitude and target altitude
are valid numbers in the range [-90, +90].
Reworked the C version of the code to be clearer:
eliminated goofy ALTDIFF macro, split out max
altitude derivative into its own function MaxAltitudeSlope,
just like the other language implementations do.
Minor rewording of comments in MaxAltitudeSlope functions.
Python InvalidBodyError now includes the invalid body
in the diagnostic message.
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.
Finding apparent angular radius by dividing body radius
by distance is not quite as accurate as taking the
arcsine of the same ratio. Changed this for a slightly
more accurate answer.
Refactored so there is a single altitude function/context pair
instead of separate ones for SearchAltitude and SearchRiseSet.
This allows InternalSearchAltitude to fully understand the
problem to be solved, and also makes the code smaller.
InternalSearchAltitude will need to perform min/max altitude
calculations, which means it needs to know all the parameters
of the altitude search directly. They are no longer hidden
inside an opaque function pointer and context.
So now SearchAltitude and SearchRiseSet no longer create
contexts, nor pass in a function pointer. They just
pass the correct numeric parameters to the generic solver
InternalSearchAltitude, which packs the parameters into
the context for Search().
Unified the two context types into a single context type,
and the two callback functions into a single callback function.
Updated the version number so I can create a new
npm package to test the pull request from @matheo
that should allow TypeScript types to be exported correctly.
Updated CodeQL config to ignore source templates,
because they are not syntactically valid source code.
Ignore other stuff that is irrelevant to published
code quality.
Made various fixes based on helpful CodeQL analysis.
The rise/set search based on hour angles is complicated,
and doesn't handle oddities that happen close to the poles.
I'm starting to rework rise/set as a more brute force solution
that iterates through finite time steps.
I also added a series of Moon data for the arctic circle,
which includes some of the more painful special cases.
For example:
Moon 130 80 2034-05-16T13:21Z s
Moon 130 80 2034-05-16T13:51Z r
Here the Moon sets, then rises 30 minutes later.
So now I'm trying to figure out how to discover
arbitrarily brief intervals like this.
I want the time increments to scale intelligently
so that we don't waste time during long periods
of inactivity (body above or below the horizon continuously),
but without missing examples like the one above.
I use a variety of doxygen versions among the systems
that run the unit tests, including my personal development
systems and GitHub Actions.
Older versions of doxygen complain about options in Doxyfile
they don't know about.
Newer versions of doxygen complain about obsolete options in Doxyfile.
I'm trying to get rid of all these noisy warnings by pruning
all options in Doxyfile that are equivalent to their default
values. I only keep the options that I need to be different
from the default.
Hopefully this quiets down the warnings.
This is my second attempt to release eclipse obscuration.
I discovered there was a missing unit test for obscuration
for lunar eclipses in the Kotlin library. It has been added.
