As mentioned in several issues (508, 464, 453), astroplan's current 1-pass, grid-based rise/set/transit methods and related others like noon() have problems with accuracy at coarse grid, and with speed at fine grid. A two-pass/optimizer approach has been suggested (issues 441 and 440) as enhancement.
I have confirmed these reports with an exhaustive test (52,707 calls) of astroplan 0.8's Observer.moon_rise_time() as example method. For an observer at Apache Point Observatory (which latitude should never see a 28-hour period without a moonrise, moonset, or transit):
100 calls (of 52,707) gave no result (masked array was returned),
140 missed the 'next' rise by a day or more,
139 missed the 'nearest' rise by a day or more, and
curiously, none missed the 'previous' rise by a day or more although there were some masked-array (missed) results.
This test was conducted on 'next', 'nearest', and 'previous' calls on every Time in grid spanning the calendar year 2021 on a grid spacing of 30 minutes, each Time dithered reproducibly (random seed set) on a uniform distribution of +- 10 minutes.
To remedy this, I have now written a two-pass moon_rise_time() demo function. In the same exhaustive test, it found the rise time in every case of 52,707. It also improved the accuracy [per moon_altaz()] (average estimated error 30 microseconds vs 0.42 seconds for astroplan), and cut the per-call run time by half (50 ms vs 108 ms). [One estimates the error as (- metric function @ time / metric function derivative @ time)].
In the proposed code, only two items distinguish all the rise/set/transit methods:
A metric function, which need only be smooth, differentiable, and must cross zero at event time in a direction that distinguishes rise from set (or transit from antitransit). E.g., for moon_rise_time(), the metric function is simply (moon altitude - horizon angle).
The zero-crossing's required direction (e.g., positive delta-altitude for rise, positive delta-azimuth for transit).
Thus the grid-and-refinement engine takes up almost all the code; an individual rise/set/transit need only define the metric function and direction of crossing and hand off to the engine. E.g., for moon_rise_time() and the above metric function, the direction is positive.
So each proposed rise/set/transit method requires only about 5 lines of code to define these two items; the engine code is absolutely the same for all of them.
The new, proposed algorithm has two passes:
Call Observer.moon_altaz() on very coarse grid of 2-hour spacing over 28 hours centered on ('nearest'), starting on ('next'), or ending on ('previous') the user's start time. It finds the bracket containing the metric function's zero-crossing, finds the best third adjacent grid point based on low-metric function value, performs a quadratic solution on the three grid points to find the two roots, chooses the correct root, and passes this approximate event time to the next pass.
Call Observer.moon_altaz() on a 3-point grid centered on the best first-pass time and spaced by about one minute. repeat the quadratic solution to find the correct root. The method immediately returns this time.
So, only two very small calls to Observer.moon_altaz() suffice, whence the speed of the proposed demo code. This approach succeeds because the of the smoothness and locally near-quadratic shape of a properly chosen metric function.
[By the way: for sun and moon rise/set/transit, a 24-hour grid duration as I think I see in astroplan 0.8 is not enough--you will definitely miss many actual events, and will misidentify the nearest event for others. The moon cycle averages about 24.8 hours, and the sun cycle is prone to edge cases at 24-hour grid search. This is probably the root cause of most of astroplan's current rise/transit/set failure to return a useful event time, as mentioned above.]
Thus,
I propose that I construct a pull request to rewrite astroplan's Observer's rise/set/transit/antitransit methods.
Questions before I would start:
Is there interest?
Should I make a first pull request of only moon_rise_time() (for your evaluation), or should I do all the methods at one time? Doing them all in one go should take little time to write, but unfortunately would take considerable time to exhaustively test each of them before making the PR.
Should I also test against package ephem? Seems like a good idea.
How should we handle cases (usually for circumpolar Observer locations or large horizon angles) where there really are no rises or sets? Shall we retain the current Warning and masked Time approach?
Is it worth capturing cases where the target rises and sets within 2 hours, which would sometimes evade the proposed grid search? This is more common for Fixed Targets; for the sun and moon it should happen only for Observers at circumpolar latitudes or for significantly non-zero horizon angles. Adding this feature would not significantly slow the algorithm, but doing it properly may lengthen the code substantially. Testing should not be much harder. Perhaps this would be reserved for a later PR.
You may have questions of me as well. I stand ready to try answering them.
As mentioned in several issues (508, 464, 453),
astroplan's current 1-pass, grid-based rise/set/transit methods and related others like noon() have problems with accuracy at coarse grid, and with speed at fine grid. A two-pass/optimizer approach has been suggested (issues 441 and 440) as enhancement.I have confirmed these reports with an exhaustive test (52,707 calls) of astroplan 0.8's Observer.moon_rise_time() as example method. For an observer at Apache Point Observatory (which latitude should never see a 28-hour period without a moonrise, moonset, or transit):
This test was conducted on 'next', 'nearest', and 'previous' calls on every
Timein grid spanning the calendar year 2021 on a grid spacing of 30 minutes, eachTimedithered reproducibly (random seed set) on a uniform distribution of +- 10 minutes.To remedy this, I have now written a two-pass moon_rise_time() demo function. In the same exhaustive test, it found the rise time in every case of 52,707. It also improved the accuracy [per moon_altaz()] (average estimated error 30 microseconds vs 0.42 seconds for astroplan), and cut the per-call run time by half (50 ms vs 108 ms). [One estimates the error as (- metric function @ time / metric function derivative @ time)].
I've put up a minimal demonstration repo at https://github.com/edose/moon_rise_time_demo .
In the proposed code, only two items distinguish all the rise/set/transit methods:
So each proposed rise/set/transit method requires only about 5 lines of code to define these two items; the engine code is absolutely the same for all of them.
The new, proposed algorithm has two passes:
So, only two very small calls to Observer.moon_altaz() suffice, whence the speed of the proposed demo code. This approach succeeds because the of the smoothness and locally near-quadratic shape of a properly chosen metric function.
[By the way: for sun and moon rise/set/transit, a 24-hour grid duration as I think I see in astroplan 0.8 is not enough--you will definitely miss many actual events, and will misidentify the nearest event for others. The moon cycle averages about 24.8 hours, and the sun cycle is prone to edge cases at 24-hour grid search. This is probably the root cause of most of astroplan's current rise/transit/set failure to return a useful event time, as mentioned above.]
Thus,
I propose that I construct a pull request to rewrite astroplan's Observer's rise/set/transit/antitransit methods.
Questions before I would start:
ephem? Seems like a good idea.You may have questions of me as well. I stand ready to try answering them.
Thank you, Eric Dose