search

pytroll / pyorbital

Orbital and astronomy computations in python
http://pyorbital.readthedocs.org/
GNU General Public License v3.0
226 stars 77 forks source link

Use Scipy brentq method instead of secant method to perform root-finding #35

Closed rdbisme closed 6 years ago

rdbisme commented 6 years ago

Most of the work is by @AndreaMalgarini

PS: I also added a fix to avoid deprecation warning for numpy indices indexing by objects that are not a tuple.

coveralls commented 6 years ago

Coverage Status

Coverage increased (+0.06%) to 83.401% when pulling 9e1c68b2e5cf705498a5f6cf254ca3250628537d on rubendibattista:hotfix/get_next_passes_brentq into 2bd8087cff802f040165c3483775359d8a61a896 on pytroll:master.

rdbisme commented 6 years ago

Should we maybe add a unit test with the problematic TLEs?

mraspaud commented 6 years ago

@rubendibattista @AndreaMalgarini thanks a lot for this ! Unit tests would be very good indeed.

mraspaud commented 6 years ago

@johnmccombs would you mind trying if this fixes your errors too ?

johnmccombs commented 6 years ago

@johnmccombs would you mind trying if this fixes your errors too ?

Hmmm, till getting a bad result. To reproduce the issue, comment out 'tlefile.fetch(TLE_FILE_NAME)' on line 19 of my sample code (issue #22), so it doesn't update the TLE file, then use the TLE file attached iridium-tle.txt

Line 1=rise, 2=approx transit, 3=PyOrbital transit, 4=set

1.   2018-03-12 14:48:35.500958 182.38926546757895 39.99999995021594
2.   2018-03-12 14:50:34.500958 248.12961834962076 89.97802252575404
3.   2018-03-12 13:59:54.556282 280.9413732517446 -85.20348517125001  
4.   2018-03-12 14:52:32.707255 1.6747595142409606 40.000000068334586

iridium-tle.txt

AndreaMalgarini commented 6 years ago

@johnmccombs I had a look at your issue since I had some free time. This is what I used:

TLE: (stored in to_test.txt)

IRIDIUM 7 [+]           
1 24793U 97020B   18137.55660552  .00000073  00000-0  19102-4 0  9994
2 24793  86.3989 179.2463 0002393  90.7421 269.4049 14.34253938100896

With the same code you used, which I report here.

from pyorbital import tlefile
from pyorbital.orbital import Orbital
from datetime import datetime

HORIZON = 40

# altitude (km) above seal-level and lat lon of the observer
ALTITUDE_ASL = 0.5
LON = 170.556
LAT = -43.368

orb = Orbital('IRIDIUM 7 [+]', 'to_test.txt')
d = datetime(2018, 3, 12, 14, 0, 15)
passage = orb.get_next_passes(d, 1, LON, LAT, ALTITUDE_ASL, horizon=HORIZON)

So basically what I found out is that by tuning the initial guess, d, you can either get a good or a bad result. For example with the value of d set as in the example you get the correct behaviour. If you change just the microseconds it diverges, for example try with d = datetime(2018, 3, 12, 14, 0, 0).

The issue lies in the get_max_parab method and it is affected by the same logic the get_root_secant had, since it depends from the initial bounds and sometimes it might diverge. A quick solution I found is to provide it the bounds of risetime and falltime, without considering the middle to reduce the search path. So doing something like that: highest = utc_time + timedelta(minutes=get_max_parab(elevation_inv, risemins, fallmins, tol=tol / 60.0))

However I am quite sure this might do more bad than good. An alternative would be either to use a more robust method or to put checks within get_max_parab to prevent it to go too far from the initial search space, and if this happens, just reset and try with a more relaxed guess.

Hope this is of help, as soon as I have more time I will try to dig deeper in this matter (since I am quite interested in it) to see if there is something to do from a mathematical point of view.

johnmccombs commented 6 years ago

Hope this is of help, as soon as I have more time I will try to dig deeper in this matter (since I am quite interested in it) to see if there is something to do from a mathematical point of view.

This is pretty much the same conclusion I came to. I did a brute-force search to solve my immediate problem and haven't had time to get back to it. Would be great if you have the time to fix it properly.

rdbisme commented 6 years ago

For what concerns this PR (maybe discuss the other issue on its thread), I'll try to setup a unit test that does not validate with current master and then validates with this PR. Give me a bit of time, I'm quite charded in this period :)

rdbisme commented 6 years ago

Done. The added unit test fails on master but validates here.

pep8speaks commented 6 years ago

Hello @rubendibattista! Thanks for updating the PR.

mraspaud commented 6 years ago

@rubendibattista @AndreaMalgarini thanks for fixing this!