Rise/set/transit time calculations and convenience methods - #14

Just in case your notifications are set differently @astroplanners/astroplanners-contributors : take a look at this sometime this week!

After discussion with @eteq today, we noted that I should relate the N parameter in _calc_riseset() to a precision of the rise/set time in units of time, and make a kwarg for sunrise/sunset that will allow users to adjust the precision (in exchange for speed).

@bmorris3 – Great!

Could you please have a look a the coverage report for the Observer class and add tests for the few missing statements that do computations (e.g. here) (tests that cover the invalid input errors aren't really necessary.)

I built the docs locally and the "Methods Summary" section for the observer class is empty. I like having the overview of available methods ... it's much quicker than scrolling through the long descriptions. Anyone know why the "Methods Summary" isn't filled?

For the docstrings ... most don't have a Returns section because it's kind of obvious what they return. But it might still be useful to mention that a Time object is returned just to be clear ... maybe in the first line of the docstring if you think a Returns section is too verbose?

@bmorris3 – What do you think is the best strategy with narrative docs: http://astroplan.readthedocs.org/en/latest/#general-documentation ?
Add / improve piece by piece with each PR or only implement the classes and write in a separate PR later?

@bmorris3 – I've left a bunch of nitpicky comments about the docstrings. If you disagree with some, OK, I don't feel very strongly about them, it's just how I would have written it.

I built the docs locally and the "Methods Summary" section for the observer class is empty.

@cdeil How are you building the docs locally? When I run python setup.py build_sphinx -lo, I get the methods summary. Were you maybe looking at Observation rather than Observer?

I'm getting an empty "Methods summary" section of the Observer class:

Weird that it's there for you ... I also use python setup.py build_sphinx -l.
Can someone else check, please?

This is what I see:

After the Hangout with @cdeil , we have some things to discuss (in particular, @cdeil and @eteq might have some things to hash out, and any input from @ejeschke and @adrn could help too!).

Precision

These methods work and will produce high precision results where necessary, but can we decide on in which use cases precision will matter the most? In our conversation, the rise/set/transit times that we thought were most important to solve with precision and with consideration for the atmosphere were for the Sun, but in almost all other cases, the rise and set times can be calculated with low precision and without the atmosphere, since almost no observers will be looking at a target on the horizon (apart from solar observers).

@cdeil suggests that I implement a simple equivalent to my fancy algorithm that gets the altitude of targets over time from simple trigonometry rather than astropy's AltAz frame, since this is the biggest bottleneck in the computation. Plus, I've already implemented Meeus' trigonometric methods before in a function that would be easy to implement in the framework already developed in this PR. However, that method relies on calculation of LST, so @eteq – is the calculation for local sidereal time any cheaper than generating AltAz frames, or is it equivalent? My past experience suggests it's faster.

Would anyone (@eteq) be opposed to imposing that when the user asks for rise/set of a target other than the sun, we default to using the trigonometric approximations? Or even for the sun for that matter?

Vectorization

@cdeil brought up the good point that we should perhaps be more concerned with properly implementing vectorization in these algorithms that with their precision for now. After this PR gets merged, maybe I should focus on implementing use cases (like computing rise/set on hundreds of targets simultaneously) rather than making super-awesome calculators.

@cdeil @eteq: Thinking on this some more, we could always use the trigonometric rise/set equations when pressure=0. That would be a simple setting for users to figure out, and it remains self-consistent.

@cdeil - I agree with you about this not being a high priority, but that's why I think it's better to go with all AltAz for now: we don't need any checking to make sure it's consistent with the "trig" answer if we don't use the trig answer at all. So I was suggesting we stick with the AltAz approach for the present and treat adding a trig answer (or better yet, figuring out where to optimize AltAz), is better done later once we can see if it's really a critical bottleneck.

Similarly, my suggestion about making sure we cover a grid of points was not meant to be exhaustive like coordinates-benchmark, but rather a relatively quick check to ensure something isn't very wrong. Like select 25 random points on the earth and 25 points in the sky or something like that. Just doing 0/0 lat/lon and 2000-1-1 is not sufficient here because most of the failure modes I know of are worse the further you go from there. Not something that should take a lot of time, just a few minutes of adding a for loop in pyephem to get the "good" answers for a few cases.

Also, @bmorris3, one related point of clarification: the approach PyEphem uses isn't actually fixed to 1980, if I understood the tweet correctly: rather it's an algorithm published in 1980. They differ from the astropy (2006) algorithms only at arcsecond levels at the very most (or centuries in the past/future), so it shouldn't be noticable at all here.

I don't agree with the "doesn't work for the eskimos and penguins" bit, though: the case I had in mind is if I'm at, say, Subaru, and ask for the rise/set time of Polaris. That's a nonsense question because it's circumpolar, and its important for observing planning to know that. and for the rise/set methods to appropriately do something in that case. Similarly, If I ask for the rise/set time of the Large Magellanic Cloud, it should also tell me that the LMC never rises at Subaru. But I have no objection to adding that later as a separate PR. But we do need to decide how we are going to tell the user about these use cases - that's a significant API choice.

@adrn @ejeschke @eteq – Do you have time for a final review of this PR? (I don't.)
It would be great if we can merge this tomorrow or Friday.

Sorry I missed the hangout yesterday!

@astroplanners/astroplanners-contributors I'll look this over today -- sorry, yesterday was a bit crazy.

@adrn If you have any questions, find me on Slack! I look forward to getting this merged and moving on.

@bmorris3 – What about @eteq's last comment? What currently happens when you ask for rise times for objects that never rise or set? Do you think this behaviour is appropriate? Or do you want to improve it here or make a new issue to discuss how to handle this? Is there a test covering this?

@cdeil I'm working on custom exceptions for these cases. There is a test checking that a ValueError is thrown for calculating the rise/set time of Polaris.

@cdeil I will implement these in a separate PR once this one is merged. I'm managing a whole bunch of overlapping branches of code, and it would be really convenient if rise/set was merged so I can build off of it.

@adrn – I'm +1 to merge this today and let @bmorris3 continue in other PRs ... I'll leave the final review and merging to you.

@bmorris3 – I didn't look yet and your exceptions to handle non-rising and non-setting objects. But raising an exception means it's not possible to process arrays of targets where only some rise, right? What numpy does there is raise a warning (which can be suppressed), but process the array and return np.nan or np.inf ... example:

>>> import numpy as np
>>> np.sqrt([4, -4, 9, -9])
__main__:1: RuntimeWarning: invalid value encountered in sqrt
array([  2.,  nan,   3.,  nan])

I just briefly looked at your tests ... doesn't look like arrays or targets are supported (or at least not tested) at the moment? Is it possible / easy to implement? Do you plan to do it in a future PR?
(one thing that can be useful is to start a new issue or PR with a task list and add the things you don't want to address in this PR, to make sure they aren't forgotten)

@cdeil The next really important thing to work on is to allow for easy vectorization of all of the methods implemented in this PR. The purpose of moving forward from this PR is to get agreement on the approach for single targets so we can move onto the more important cases with multiple targets.

I have a few minor inline comments above, but a larger issue is that I think there should be many more tests. The tests right now are largely functional tests of end-to-end calculations and comparison with pyephem. It would be good to add a bunch of unit tests (as @cdeil hinted at in a previous comment) -- e.g., testing your "internal" functions like _horiz_cross on many different possible inputs, and checking that it fails where you expect. Now, normally I would say that this is required before the PR is accepted, but I do agree that we need to find a balance here for moving forward / writing tests / merging code that other stuff depends on...soooooo I guessssss I'd be willing to merge this (after all comments are addressed) and have a separate PR expanding the tests.

BTW there were weird test failures (I think travis issues) so I restarted the tests.

@adrn It looks like travis issues are still happening.

Merged, but let's all keep in mind:

edge cases (objects never rise / set)
#14 (comment)
#18
#19