I tried more distant objects like Jupiter ... Neptune.
This revealed that at increasing distances, the convergence
threshold in inverse_terra needed to increased also.
So now I use 1 AU as a baseline, and scale up linearly
for more distant objects.
Asking the latitude and longitude directly beneath
the Sun causes inverse_terra not to converge, because the
convergence increment `W` never got below 1.48e-8, but the
convergence limit was 1.0e-8. I increased the limit to 2.0e-8
in all programming language versions.
I'm hoping that is a big enough tolerance for all cases now,
but I will do more testing to see if further fixes are required
for even more distant bodies than the Sun.
Prevents crashes due to dereferencing NULL time pointers.
Passing in NULL for a `time` pointer will no longer cause
a crash in an Astronomy Engine function.
Wherever possible, a NULL time pointer will result in a
status code `ASTRO_INVALID_PARAMETER`.
`Astronomy_Horizon` has no way to report a status code,
so a null pointer causes it to return all NAN values.
Perhaps it should return a status code (considering for separate commit).
Thanks to [Steven Booth](https://github.com/sbooth) for suggesting this!
It turns out that GetSystemTimeAsFileTime only returns
time with millisecond resolution.
In order to get microsecond resolution in Astronomy_CurrentTime(),
I had to switch to GetSystemTimePreciseAsFileTime for Windows.
Example output from unit test:
C Test_AstroTime: PASS - realtime increment = 3.143e-07 seconds after 1 iterations.
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.
This patch adds support for microsecond time resolution on the UNIX and
WIN32 platforms.
Implementation details:
Previously, the `Astronomy_CurrentTime()` function used `time(NULL)` to
get the current time, thus providing 1-second resolution. This patch
uses `gettimeofday()` for UNIX and `GetSystemTimeAsFileTime()` for WIN32
platforms. If neither are supported, it will fall back to `time()` and
issue a #warning. (Note that windows.h defines ARRAYSIZE, which may or
may not be compatible. Thus, the internal define `ARRAYSIZE` is renamed
to `ASTRO_ARRAYSIZE`.)
The UNIX code was tested in Linux, and in arm-eabi-none under newlib.
The WIN32 code was tested using MinGW/64 under WINE.
Use case:
One-second resolution is enough in most cases. However, there are cases
where higher resolutions are desirable. For example:
We are using the astronomy.c library to control a mechanical rotor to
track celestial objects (planets, stars, and satellites). The rotor
controller uses a PID controller with 100 tick/sec updates, and tracks
the velocity of the azimuth and elevation angles from the previous tick.
With 1-second resolution, the PID controller jerks and oscillates once
per second as it adjusts to the new position. With at least
10-millisecond resolution (100/sec), it can calculate the per-tick
velocity change and track smoothly with far less jitter.
More information about the project:
Source using astronomy.c:
https://github.com/KJ7NLL/space-ham
Lego-controlled az/el:
https://youtu.be/vrlw4QPKMRY
Lego-controlled telescope focus:
https://youtu.be/p-5dOQG95xg
APID+SMC control algorithm:
https://doi.org/10.1016/j.precisioneng.2022.01.006
Signed-off-by: Eric Wheeler, KJ7LNW <astronomy-git@z.ewheeler.org>
Tested-by: Zeke Wheeler, KJ7NLL <kj7nll@gmail.com>
Slightly different cppcheck dev 2.11 behaviors have added
another warning that I don't care about. I don't want to
have to convert callback pointers to const, then cast them
to const.
However, it did find a couple of useful cases I fixed in
astronomy.c where GravSim parameters could be made const.
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.
The C functions for calculating calendar dates used the
type `long` to perform calculations that require 64-bit
integers. However, in some C compilers, `long` is still
32 bits. This caused a failure in Windows for extreme
year values. So I now use the type `int64_t` to explicitly
require a 64-bit integer.
Fixed problems converting AstroTime to calendar dates and back.
Also expose struct CalendarDateTime to outside callers,
for convenience dealing with Gregorian calendar dates.
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.
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.
I'm taking gcc 12.2 for a test drive today.
It reports a few warnings that slipped through earlier versions.
None of the warnings concern me for actual code safety,
but I went ahead and resolved them to keep the build clean.
Also provide a hook for a CPP environment variable to override
the C++ compiler to use, instead of forcing g++.
Define ECT = True Ecliptic of Date in the documentation.
I will soon convert the Ecliptic() functions to return ECT instead of
ECL, but I will retain ECL support via rotation matrix functions.
Somehow gcc didn't warn me that a `void` function
was trying to return a `void` value. It makes sense
from a functional language "unit type" perspective,
but it wasn't intentional, and it is weird. Fixed it.
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.
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.
