I was able to eyeball the slope from a graph of deduced
phase curve by analyzing JPL Horizons output.
Now my Pluto magnitude values are well within agreement.
I found a paper by James L Hilton (USNO) that provides
formulas for the phase curves of Mercury and Venus that match
the JPL Horizons tool within 0.012 mag.
https://iopscience.iop.org/article/10.1086/430212
JPL Horizons data does not include ring aspect for Saturn,
so don't bother plotting the Saturn data.
I already deleted the Saturn data so it's not there anymore.
The JPL Horizons data for Saturn says at the bottom of the file
that it does not include variations due to ring aspect.
Therefore, it is not suitable for my needs.
I realized if I'm going to create regression models of planets,
it will be helpful to have known correct values of heliocentric
distance of the planet also. Will need to figure out what to do
for Saturn's rings.
Not fully validated, but I did tweak Montenbruck/Pfleger
formulas to match JPL Horizons output for 2018-04-27 02:00 UTC.
Still need to implement formulas for Moon and Saturn.
"Phase angle" means the angle between the Sun and the Earth
as seen by a third body.
This function calculates the angle between the Sun and a body
as seen by the Earth, with a range that goes all the way to
360 degrees, allowing finding all 4 quarters of the Moon's cycle.
When searching for oppositions and conjunctions of Mercury and Mars,
dynamically adjust to their eccentric orbits by tweaking the
effective synodic period based on how far we missed the mark
on each iteration.
This brings the average for Mercury from 19 down to 6.4.
All the other planets got at least a little better.
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.
Compiler does not understand that there is no need for
a switch default because the enclosing switch statement
has already guarded against that. Changed to ternary
operator to make compiler warning go away.
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.
I learned by debugging the code that the Chebyshev approximation
does not exactly reproduce the endpoints like I thought. It does so
for the x values indicated by Alpha[1][...].
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.
