Confirmed that Mercury is taking 19 iterations/call on average,
and Mars is taking about 10 iterations/call.
The other planets average 6 iterations/call.
Still need to verify that there is a consistent interval between
consecutive events.
Mercury is taking way too long to converge.
Mars is kind of slow too.
Need to improve the efficiency of SearchRelativeLongitude!
Using relative heliocentric ecliptic longitude of the Earth
and the other planet. Home in on when both planets have the
same longitude (the difference is 0).
For some reason, all my calculations are about 8 minutes earlier
than predictions from the test data. I suspect this is because
of light travel time from the Sun (equivalently, aberration).
Added new function Astronomy.SunPosition().
It is supposed to return ecliptic coordinates of date for the Sun
as seen from the center of the Earth.
The values look reasonable but I need to test them.
Will use the Sun's longitude in the return value from SunPosition()
to determine solstices and equinoxes.
All callers of sidereal_time ended up needing it for apparent time,
not mean time. So I simplified the code so it no longer has extra
stuff for calculating GMST.
Astronomy.GeoVector no longer iterates to try to correct light
travel time for the apparent position of the Sun. The Sun's
heliocentric coordinates are always (0,0,0), so there is no need
to do that.
Search limit can be adjusted in options passed into Search().
After 20 iterations, we should have divided the search
region by a factor of more than a million. If quadratic
interpolation can't finish the job at that point, something
is really wrong.
Allow caller to pass in pre-evaluated endpoints to begin the search.
This eliminates 2 function calls per search, reducing the
average from 8 calls/search down to 6 calls/search.
I think this is about as good of performance as I'm going to get.
The smaller the slope magnitude |df/dt| is, the larger
the uncertainty in dt. That means we are better off using
an estimated value for the slope each time than underestimating
the time error like we were doing.
This also simplifies the code, and does not make it very
much slower.
Now pass in max slope of function to be searched, expressed
in units/day. By seeing how far the function is from zero,
we can deduce whether we are within the specified time tolerance
of finding the event.
Use a simplified refraction model in the rise/set search so that
the function is better fit by parabolas. Assume constant refraction
instead of variable refraction, because it only matters near the horizon
anyway. Use a canned value of +34 minutes, which creates close fit with
test data.
I am interested in optimizing the Search function.
Right now it is a very simple binary search that keeps breaking
an interval in half to narrow in on the time of where the supplied
function ascends through zero. I know this can be made much better,
and this is important because the function calls are very expensive
in some cases.
So this commit adds the beginning of some simple metrics tracking
where unit test code can retrieve the number of times Search
sampled the function it is trying to find the ascending root for.
The refraction formula went nuts near altitude angle -5.11 degrees.
We were taking the tangent of a value that zoomed toward infinity
near that value, causing essentially random numbers without any
upper bound to their size. Just like JPL Horizons, truncate any
angle more than 1 degree below the horizon, only I have a linear
taper down to 0 refraction as the altitude angle approaches -90.
I did not want any chance of creating an altitude less than -90.
Removed unit tests for the Sun at latitude -80 degrees.
It is too easy for my code to behave differently from another
calculator, because tiny changes in atmospheric modeling can
cause disagreement about whether there even is a sunset/sunrise.
This is because for observers so close to the pole, the Sun
sometimes barely dips below the horizon and then comes back
up for less than an hour.
This is the first time it has passed the unit test,
although the unit test is just exercising whether the predictions
occur in the right order. I will need to add check for how accurate
the predictions are.
This way people don't have to figure out how to iterate
through moon quarters. Use SearchMoonQuarter to start iteration,
NextMoonQuarter to iterate through as many more as desired.
Can now search for the next new moon, first quarter,
full moon, or third quarter.
Verified against US Navy Observatory data.
Predictions are confirmed to within 2 minutes of time
for years between 1800 and 2100.
The floating point constants emitted for the Pluto Chebyshev
model did not quite match between Linux and Windows.
I suspect the real problem is not the operating system
but that I'm using different versions of Node on both:
Windows: v10.15.3
Linux: v8.15.1
Now I print only 12 place after the decimal instead of 18.
This makes no difference in the unit test output,
and reduces the JS code size significantly.
Expanded the Chebyshev model for Pluto.
Resampled VSOP models to have required accuracy over wider date range.
Decreased astro_check.js sampling rate to allow tests to run faster.
Now the JavaScript code uses UT and TT values expressed
in days since 2000, instead of Julian Dates.
This makes the numeric values much smaller and thus
should yield less floating point error when time solvers
are added later.
Using historic, recent, and predicted values of TT-UT instead of
UTC leap seconds. With linear interpolation, there are no longer
discrete jumps in the calculated TT values. Hopefully, this will
make event solvers (rise, set, etc) more well-behaved.
Astronomy.GeoVector now corrects for light travel time from
the observed object. This reduced worst case angular error
from 1.16 arcmin to 0.89 arcmin (0.27 arcmin improvement).
I found some online resources that helped me track down the
formula for the refraction model used in the JPL Horizons
online tool. Now the JavaScript library allows 4 different
refraction options in Astronomy.Horizon():
false : no refraction
'novas' : use the NOVAS C 3.1 algorithm.
'jplhor' : JPL Horizons algorithm, clamped beyond 1 degree below horizon.
'sae' : same as 'jplhor', only without clamping.
Now passes the jplcheck unit test without filtering out objects below the horizon!
Always compile the C code when executing the script './run'.
