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.
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.
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.
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.
For consistency, I now calculate the Sun's apparent position
correcting for aberration, for lunar eclipses and transits.
This didn't make much difference for accuracy.
Before correcting for Sun aberration:
$ ./ctest transit
C TransitFile(eclipse/mercury.txt ): PASS - verified 94, max minutes = 10.709, max sep arcmin = 0.2120
C TransitFile(eclipse/venus.txt ): PASS - verified 30, max minutes = 9.108, max sep arcmin = 0.6772
$ ./ctest -v lunar_fraction
C LunarFractionCase(2010-06-26) fraction error = 0.00900991
C LunarFractionCase(2012-06-04) fraction error = 0.00491535
C LunarFractionCase(2013-04-25) fraction error = 0.00143039
C LunarFractionCase(2017-08-07) fraction error = 0.00682336
C LunarFractionCase(2019-07-16) fraction error = 0.00572727
C LunarFractionCase(2021-11-19) fraction error = 0.00350680
C LunarFractionCase(2023-10-28) fraction error = 0.00370518
C LunarFractionCase(2024-09-18) fraction error = 0.00429906
C LunarFractionCase(2026-08-28) fraction error = 0.00322697
C LunarFractionCase(2028-01-12) fraction error = 0.00405870
C LunarFractionCase(2028-07-06) fraction error = 0.00857840
C LunarFractionCase(2030-06-15) fraction error = 0.00557106
C LunarFractionTest: PASS
After correcting for Sun aberration:
$ ./ctest transit
C TransitFile(eclipse/mercury.txt ): PASS - verified 94, max minutes = 10.709, max sep arcmin = 0.2120
C TransitFile(eclipse/venus.txt ): PASS - verified 30, max minutes = 9.108, max sep arcmin = 0.6772
$ ./ctest -v lunar_fraction
C LunarFractionCase(2010-06-26) fraction error = 0.00762932
C LunarFractionCase(2012-06-04) fraction error = 0.00606322
C LunarFractionCase(2013-04-25) fraction error = 0.00111560
C LunarFractionCase(2017-08-07) fraction error = 0.00571542
C LunarFractionCase(2019-07-16) fraction error = 0.00713913
C LunarFractionCase(2021-11-19) fraction error = 0.00298979
C LunarFractionCase(2023-10-28) fraction error = 0.00448445
C LunarFractionCase(2024-09-18) fraction error = 0.00367044
C LunarFractionCase(2026-08-28) fraction error = 0.00405559
C LunarFractionCase(2028-01-12) fraction error = 0.00347340
C LunarFractionCase(2028-07-06) fraction error = 0.00729982
C LunarFractionCase(2030-06-15) fraction error = 0.00680776
C LunarFractionTest: PASS
C function Astronomy_SearchLocalSolarEclipse now reports
obscuration for the peak of a solar eclipse seen at a given location.
Some of the test data from aa.usno.navy.mil looks suspiciously inaccurate.
I am finding better agreement with Fred Espenak's eclipsewise.com and
Xavier M. Jubier's xjubier.free.fr online calculators.
The larger obscuration differences are from aa.usno.navy.mil:
don@spearmint:~/github/astronomy/generate $ ./ctest -v solar_fraction
C GlobalAnnularCase(2023-10-14) obscuration error = 0.00009036, expected = 0.90638000, calculated = 0.90647036
C GlobalAnnularCase(2024-10-02) obscuration error = 0.00005246, expected = 0.86975000, calculated = 0.86980246
C GlobalAnnularCase(2027-02-06) obscuration error = 0.00007237, expected = 0.86139000, calculated = 0.86146237
C GlobalAnnularCase(2028-01-26) obscuration error = 0.00003656, expected = 0.84787000, calculated = 0.84790656
C GlobalAnnularCase(2030-06-01) obscuration error = 0.00008605, expected = 0.89163000, calculated = 0.89171605
C LocalSolarCase(2023-10-14) obscuration diff = 0.00006323, expected = 0.90638000, calculated = 0.90644323
C LocalSolarCase(2023-10-14) obscuration diff = -0.00521043, expected = 0.57800000, calculated = 0.57278957
C LocalSolarCase(2024-04-08) obscuration diff = 0.00000000, expected = 1.00000000, calculated = 1.00000000
C LocalSolarCase(2024-04-08) obscuration diff = 0.00304558, expected = 0.34000000, calculated = 0.34304558
C LocalSolarCase(2024-10-02) obscuration diff = 0.00007858, expected = 0.86975000, calculated = 0.86982858
C LocalSolarCase(2024-10-02) obscuration diff = -0.00343797, expected = 0.43600000, calculated = 0.43256203
C LocalSolarCase(2030-06-01) obscuration diff = 0.00007259, expected = 0.89163000, calculated = 0.89170259
C LocalSolarCase(2030-06-01) obscuration diff = -0.00059871, expected = 0.67240000, calculated = 0.67180129
C SolarFractionTest: PASS
I am going to rework using xjubier.free.fr in a subsequent commit.
Fixed bug in internal C function Obscuration().
It was not calculating obscuration correctly when
the second disc is completely inside the first disc,
i.e. the annular solar eclipse case.
Added C/C++ calculation of obscuration for global solar
eclipses that are total or annular at their peak
location and time. Obscuration is not calculated
for partial solar eclipses by SearchGlobalSolarEclipse.
Soon SearchLocalSolarEclipse will calculate obscuration
for partial solar eclipses at the geographic location
provided by the caller.
Starting to add support for calculating the intensity
of lunar eclipses and solar eclipses in terms of "obscuration".
This commit adds calculation of obscuration for lunar eclipses
in C/C++. The structure returned by SearchLunarEclipse and
NextLunarEclipse now includes an `obscuration` field whose value
is in the range [0, 1], indicating what fraction of the Moon's
apparent disc is covered by the Earth's umbra at the eclipse's peak.
I ported the NOVAS C 3.1 functions julian_date and cal_date to Python,
and removed the dependence on the standard datetime class for calculating UT.
Now we can create Time objects for a much wider range of year values.
Simplified the julian_date formula in C and C#.
In the Python version, I had to account for a difference
in the way integer division works for negative numbers.
In Python, integer division always rounds down, not toward
zero like it does in C/C#. So I reworked the formulas to
avoid dividing a negative integer (month-14), dividing the
positive quantity (14-month) instead and toggling addition
of the term with subtraction of the term.
I use the reworked (14-month) version in C and C# for consistency.
Also, the formatting of the formula was wacky and didn't make sense,
so now it easier to read and understand.
The Python regex for parsing dates has been expanded to allow
years before 0 and after 9999.
Allow converting Python Time to string for years before 0 and after 9999.
The .NET type System.DateTime is limited to the years 0000..9999.
The Astronomy Engine type AstroTime was using System.DateTime to
convert a (year, month, day, hour, minute, second) tuple into a
fractional day value. This caused an exception for years outside
the supported range 0000..9999.
I ported the NOVAS C 3.1 functions julian_date and cal_date to C#,
and removed the dependence on System.DateTime.
Now we can create AstroTime for a much wider range of year values.
Allow converting AstroTime to string for years before 0 and after 9999.
The following C# functions now support searching
in forward or reverse chronological order:
Astronomy.SearchRiseSet
Astronomy.SearchAltitude
Astronomy.SearchHourAngle
Also fixed places where I forgot to update documentation
for the corresponding changes to the C code.
The following C functions now support searching
in forward or reverse chronological order:
Astronomy_SearchRiseSet
Astronomy_SearchAltitude
Astronomy_SearchHourAngleEx
The function Astronomy_SearchHourAngleEx replaces
Astronomy_SearchHourAngle, adding a new `direction` parameter.
A #define for Astronomy_SearchHourAngle preserves backward
compatibility for older code.
Implementation notes:
Astronomy_SearchRiseSet and Astronomy_SearchAltitude used
to call a private function InternalSearchAltitude.
That function has been split into two functions:
BackwardSearchAltitude and ForwardSearchAltitude,
for searching both directions in time.
Fixed a bug where it was possible to report a successful
altitude event that went outside the time limit specified
by `limitDays`.
The quadratic interpolator used by `Search` was returning
an unused output: `x`. Search does not need this dimensionless
value; it only cares about the time solution `t` and the slope
of the function at `t`. Removed `x` from the return value
to make slightly smaller/faster code.
Enhanced the C function Astronomy_SearchMoonPhase
to allow searching forward in time when the `limitDays`
argument is positive, or backward in time when `limitDays`
is negative.
Added unit test "moon_reverse" to verify this new feature.
The C version of Astronomy Engine does not work correctly
when gcc "fast math" optimizations are enabled.
The problem is that Astronomy Engine uses NAN values to
represent invalid/uninitialized floating point numbers.
The -Ofast option breaks the ability of the runtime to
check for NAN values, resulting in multiple failures
and incorrect behaviors at runtime.
Added a compile-time check for the __FAST_MATH__ preprocessor
symbol, which gcc defines to signal that the optimization
was enabled. If detected, this results in a compiler error
to make it obvious that something is wrong before invalid
code would be executed.
This is not an ideal fix for two reasons:
1. I don't know if this will detect similar problems for
other compilers than gcc.
2. If individual risky math optimizations are enabled, instead
of the combination of options included in -Ofast, the
__FAST_MATH__ preprocessor symbol will not be defined
and no compiler error will occur. I could not find a
way to detect individual risky optimizations.
However, this change is much better than nothing, and
hopefully it will prevent most cases of overly-aggressive
optimization.
Implemented two new C functions for generic correction
of light travel time:
1. Astronomy_CorrectLightTravel
This is a completely generic solver that is passed
a pointer to a function that returns a relative
position vector for a given time.
The generic solver keeps calling the function
whose address is passed, until the light travel
solution converges. This usually takes 3 or 4 iterations.
2. Astronomy_BackdatePosition
This is a more specific solver that uses Astronomy_CorrectLightTravel
to backdate the position of a target body as seen from an
observer body at a specified observation time.
It implements the aberration option needed by Astronomy_HelioVector.
Astronomy_HelioVector now uses Astronomy_BackdatePosition,
instead of having a redundant light travel solver,
so the code isn't much larger than it was before.
Corrected a mistake in the explanation of the
C function Astronomy_GravSimInit: the `bodyStates`
parameter is NOT barycentric -- it is relative to the
originBody parameter.
Python had improperly formatted documentation for
Time.FromTerrestrialTime parameter `tt`.
The Python markdown generator `pydown` did not
correctly handle links to compound symbols like
`#GravitySimulator.Update`. It also was trying
to link to `StateVector[]` instead of `StateVector`.
Removed unnecessary and unhelpful documentation
for C# internal class constructors. They do not appear
in the generated markdown documentation anyway.
Other minor wording revisions in the documentation.
Finished coding the Python version of the gravity simulator.
No unit tests have been written yet.
Cleaned up documentation in the other languages.
Made some functions static that did not need to be members.
Added the following functions:
GravSimTime
GravSimNumBodies
GravSimOrigin
GravSimSwap
The GravSimSwap function allows an instantaneous "undo"
of a simulation step, which required refactoring the internal
data structures of the simulator. I did this because I realized
I needed a way to undo exploration of time steps near a fixed
current time. This will make it easier to implement a
light travel time solver.
I have decided the ability to select different
collections of gravitating bodies causes far
more complexity in the code than it is worth.
So now the gravity simulator always calculates
the Sun and all planets except Pluto.
This greatly simplifies the core code, gets
a good balance between efficiency and accuracy,
and makes the test matrix much simpler.
The GravSimOriginState function was doing too
much work for bodies that weren't already calculated
and cached. Instead of always calling Astronomy_BaryState,
it is possible to do a lot less work for most bodies
by recycling the SSB calculation that has already
been done.
Other minor code cleanup -- mostly stuff that makes
the code easier to read and understand.
Starting implementation of a generalized gravity simulator.
Already calculating the movement of Ceres, but with less
accuracy than I had hoped. I don't know if the lack of modeling
pull of the other asteroids has a larger effect than I expected,
or there is just something wrong with the implementation.
