skyfielders / python-skyfield

Elegant astronomy for Python
MIT License
1.41k stars 212 forks source link

find_discrete results apparently depend on the Skyfield version #872

Closed aendie closed 1 year ago

aendie commented 1 year ago

I have a big issue with almanac.find_discrete which I need to describe in detail. Of course I am looking into my particular issue and I am unaware of other issues that might have been fixed in the meantime. My example relates to one particular condition: MARS rise and set times at latitude 64.0°N. That sounds very specific. I have not found any issues with Mercury, Venus, Jupiter or Saturn. I have no issues with Mars at lower latitudes. So here's how the problem originally surfaced....

I am making graphs of planet visibility and Mars looks okay at latitude 63.0N:

Mars 2023 latitude 63 0N

There is a small mistake (a shading omission) in the graph above but that's for me to fix. I then tried latitude 64.0N and everything broke. Looking at the graphs without shading I saw a problem immediately:

Mars 2023 latitude 64 0N

The Mars SET between Jan 1 and Mar 2 has a gap in it. So does Mars RISE in the same timeframe. It's between Feb 24 and Feb 26 inclusive. In order to investigate this I developed a test program (included below) which outputs Mars' rise and set times in chronological order per day. I also include a calculated 'hbh' field which prints the hours the planet is below the horizon. This field supports my argument why I believe the data is incorrect.

Next I had to investigate if this issue is reproducible in Windows (on Intel) as well as Ubuntu (on AMD). It is! Eventually I found it is tightly connected with the Skyfield version. I also tried the de440s_plus_MarsPC.bsp ephemeris besides using de421.bsp... there was no difference: both had the same issue.

Looking at the output from my test program with Skyfield 1.46 you see:

Linux - Skyfield 1 46

Up to Feb 23 (inclusive) the 'hbh' span slow reduces to Mars being 0.96 hours below the horizon on Feb 23. Then it jumps suddenly to zero for the next three days. And from Feb 27 to Mar 2 it's well-behaved again going down from 0.65 to 0.28 before Mars correctly enters the range Mar 3 to Apr 4 where Mars is above the horizon all day. The issue is with Feb 24 to 26 inclusive. It even suggests that the correct data probably gradually steps 'hbh' down from 0.96 to 0.65 hours.

Next I began to look for the results from other Skyfield versions .... all screenshots are from Ubuntu Desktop 22.04.2 LTS. Skyfield 1.35 gives these results:

Linux - Skyfield 1 35

Note that there is no "zero hours above horizon" range gap during February but I find the final Feb 23 hbh (0.96 hours) and the initial Apr 5 hbh (0.65 hours) rather high on either side of the time when Mars never sets. Skyfield 1.40 gives me the same results. Skyfield 1.41 on the other hand gives me beautiful data:

Linux - Skyfield 1 41

And one sees that the Feb 23 to Apr 5 range in previous versions excluded several days of significant data. The final Mar 2 hbh is just 0.28 hours and Apr 3 just 0.35 hours. This has reduced the "days above horizon" from 40 days to 31 days! And the data just looks and feels good. Skyfield 1.42 and 1.43 give the same good results. Skyfield 1.44 is where things break:

Linux - Skyfield 1 44

This introduces two days (Feb 24 and 25) with no rise or set (Mars still above the horizon) which is illogical. Skyfield 1.45 is the same. Skyfield 1.46 goes a step worse with three days (Feb 24, 25, 26) with no rise or set (Mars still above the horizon).

Linux - Skyfield 1 46

Please pardon the exhaustive analysis, but it might prove useful in fixing this problem.

Finally I include here a ZIP file that contains my test program that gave me the above results. Please forgive my coding - it was put together in a hurry. I thank all who look into this deeper and I look forward to future versions of Skyfield! Sincere Regards, Andrew

Mars.zip

brandon-rhodes commented 1 year ago

Alas, the logic used by find_discrete() routine has not changed between 1.35 and 1.46, so we must look elsewhere for the reason you are getting different results. Your search technique looks like it would be a bit fragile, so my first guess is that updates to the ∆T table are what changes its results. And it turns out that's indeed the cause! Running Skyfield from the repository, I get no gaps when I run your script immediately after:

git co d5e627c

but get the gaps when I try the next revision:

git co 914f87d

The commit in question updated the ∆T table right before the 1.44 release:

914f87d

So, how is your search routine fragile, so that slight differences in its starting conditions can throw it off? I think it's because of this setting:

    is_obj_up_at.rough_period = 0.5  # twice a day

That setting does not, in fact, mean ‘this happens twice a day’ — but, rather, ‘check the value of this variable every 4.0 * rough_period hours’ — a confusing meaning, that has resulted in rough_period getting deprecated back in Skyfield 1.23; the documentation now recommends a much more straightforward value step_days which lets the step size between samples be specified directly.

So by setting rough_period to 0.5, you are getting a sample taken every 2 hours. Which means that when Mars is up for at least two hours, you are guaranteed to find out, since one of your samples will land inside and one outside the state where Mars is up.

But if Mars is up for less than 2 hours, you are taking a gamble. Maybe you'll be lucky and, as Skyfield takes samples 2 hours apart, it will happen to see Mars change state from up to down. But when you're not lucky, then, as you see in your examples, you can miss Mars being up. So I'd recommend two changes:

  1. Abandon rough_period and use step_days.
  2. If you want to guarantee seeing Mars up for even, say, 15 minutes a day, then your step size needs to be 15 minutes — or, since it's specified in days, 1.0 / 24.0 / 60.0 * 15.0. Once Mars is up even more briefly, you'll again be gambling with whether you miss it.

Another possibility would be to grab the whole year, and run find_maxima to find all the times when Mars is highest. That way, it doesn't matter whether Mars is up at all, since it will still experience a maximum and a minimum altitude about 12 hours apart, even if both of them are below the horizon, so there isn't a danger of missing those events. Then, filter the resulting list and look for rises and sets around the maxima that actually bring Mars above the horizon.

There's an even better approach towards risings-and-settings than that, that I started coding up once, but I can't remember in what state I left the routine. I'll take a look tomorrow and let you know what I find! But these ideas should at least get your investigation moving again.

aendie commented 1 year ago

Many many thanks for your immediate response. I really appreciate that. Apologies for the poor title I originally gave this - I have changed it into a more realistic title (I hope).

Your ideas are very valuable and I fully agree with you that the problem must lie elsewhere - independent of the Skyfield version.

Now, I have zero experience with almanac.risings_and_settings. I implemented it in my program more "by guesswork" and was quite surprised to get dates that agreed absolutely correctly with the best results I got above with Skyfield 1.41, 1.42 and 1.43.

Now this is where I begin to feel incredibly stupid...

Using almanac.risings_and_settings leaves even larger gaps in the expected data! Am I doing something wrong? Here are my results using Skyfield 1.46. Note that all data with set and rise values AND all data from April 13th onwards agrees exactly with the best results I got above using my day5length() function above.

D:\_DEVelopment\Astronomical coding\Z_STUFF\z_Mars>py Mars.py
default latitude: 64.0°N
to pick an alternative latitude use a command line argument, e.g. >py Mars.py -lat 63.0N

  Enter year as 'YYYY':

        hbh = 'hours below horizon (per day)'
Jan 1   hbh =  2.09  set: 2023-01-01, 08:41:25  rise: 2023-01-01, 10:46:54
Jan 2   hbh =  0.00
Jan 3   hbh =  0.00
Jan 4   hbh =  0.00
Jan 5   hbh =  0.00
Jan 6   hbh =  0.00
Jan 7   hbh =  0.00
Jan 8   hbh =  0.00
Jan 9   hbh =  0.00
Jan 10  hbh =  2.22  set: 2023-01-10, 07:58:33  rise: 2023-01-10, 10:11:43
Jan 11  hbh =  2.23  set: 2023-01-11, 07:54:20  rise: 2023-01-11, 10:07:51
Jan 12  hbh =  2.23  set: 2023-01-12, 07:50:12  rise: 2023-01-12, 10:04:00
Jan 13  hbh =  2.23  set: 2023-01-13, 07:46:11  rise: 2023-01-13, 10:00:09
Jan 14  hbh =  2.23  set: 2023-01-14, 07:42:17  rise: 2023-01-14, 09:56:18
Jan 15  hbh =  2.23  set: 2023-01-15, 07:38:28  rise: 2023-01-15, 09:52:28
Jan 16  hbh =  2.23  set: 2023-01-16, 07:34:46  rise: 2023-01-16, 09:48:38
Jan 17  hbh =  0.00
Jan 18  hbh =  0.00
    ... "hours below horizon" zero until ...
Apr 8   hbh =  0.00
Apr 9   hbh =  0.00
Apr 10  hbh =  1.16  set: 2023-04-10, 04:47:20  rise: 2023-04-10, 05:57:03
Apr 11  hbh =  1.25  set: 2023-04-11, 04:42:58  rise: 2023-04-11, 05:58:07
Apr 12  hbh =  0.00
Apr 13  hbh =  1.43  set: 2023-04-13, 04:34:26  rise: 2023-04-13, 06:00:04
Apr 14  hbh =  1.51  set: 2023-04-14, 04:30:14  rise: 2023-04-14, 06:00:59

The data above with no RISE or SET, up to April 13th 2023, is absolutely incorrect. For reference, it should be the following (which I created using my original program with Skyfield 1.43):

D:\_DEVelopment\Astronomical coding\Z_STUFF\z_Mars>py Mars_v0.py
default latitude: 64.0°N
to pick an alternative latitude use a command line argument, e.g. >py Mars.py -lat 63.0N

  Enter year as 'YYYY':

        hbh = 'hours below horizon (per day)'
Jan 1   hbh =  2.09  set: 2023-01-01, 08:41:25  rise: 2023-01-01, 10:46:54
Jan 2   hbh =  2.11  set: 2023-01-02, 08:36:14  rise: 2023-01-02, 10:42:56
Jan 3   hbh =  2.13  set: 2023-01-03, 08:31:09  rise: 2023-01-03, 10:39:00
Jan 4   hbh =  2.15  set: 2023-01-04, 08:26:10  rise: 2023-01-04, 10:35:04
Jan 5   hbh =  2.16  set: 2023-01-05, 08:21:18  rise: 2023-01-05, 10:31:09
Jan 6   hbh =  2.18  set: 2023-01-06, 08:16:32  rise: 2023-01-06, 10:27:15
Jan 7   hbh =  2.19  set: 2023-01-07, 08:11:53  rise: 2023-01-07, 10:23:21
Jan 8   hbh =  2.20  set: 2023-01-08, 08:07:20  rise: 2023-01-08, 10:19:28
Jan 9   hbh =  2.21  set: 2023-01-09, 08:02:53  rise: 2023-01-09, 10:15:35
Jan 10  hbh =  2.22  set: 2023-01-10, 07:58:33  rise: 2023-01-10, 10:11:43
Jan 11  hbh =  2.23  set: 2023-01-11, 07:54:20  rise: 2023-01-11, 10:07:51
Jan 12  hbh =  2.23  set: 2023-01-12, 07:50:12  rise: 2023-01-12, 10:04:00
Jan 13  hbh =  2.23  set: 2023-01-13, 07:46:11  rise: 2023-01-13, 10:00:09
Jan 14  hbh =  2.23  set: 2023-01-14, 07:42:17  rise: 2023-01-14, 09:56:18
Jan 15  hbh =  2.23  set: 2023-01-15, 07:38:28  rise: 2023-01-15, 09:52:28
Jan 16  hbh =  2.23  set: 2023-01-16, 07:34:46  rise: 2023-01-16, 09:48:38
Jan 17  hbh =  2.23  set: 2023-01-17, 07:31:11  rise: 2023-01-17, 09:44:48
Jan 18  hbh =  2.22  set: 2023-01-18, 07:27:41  rise: 2023-01-18, 09:40:58
Jan 19  hbh =  2.21  set: 2023-01-19, 07:24:18  rise: 2023-01-19, 09:37:08
Jan 20  hbh =  2.21  set: 2023-01-20, 07:21:00  rise: 2023-01-20, 09:33:18
Jan 21  hbh =  2.19  set: 2023-01-21, 07:17:49  rise: 2023-01-21, 09:29:28
Jan 22  hbh =  2.18  set: 2023-01-22, 07:14:44  rise: 2023-01-22, 09:25:38
Jan 23  hbh =  2.17  set: 2023-01-23, 07:11:44  rise: 2023-01-23, 09:21:48
Jan 24  hbh =  2.15  set: 2023-01-24, 07:08:51  rise: 2023-01-24, 09:17:58
Jan 25  hbh =  2.13  set: 2023-01-25, 07:06:03  rise: 2023-01-25, 09:14:08
Jan 26  hbh =  2.12  set: 2023-01-26, 07:03:21  rise: 2023-01-26, 09:10:17
Jan 27  hbh =  2.10  set: 2023-01-27, 07:00:44  rise: 2023-01-27, 09:06:27
Jan 28  hbh =  2.07  set: 2023-01-28, 06:58:13  rise: 2023-01-28, 09:02:36
Jan 29  hbh =  2.05  set: 2023-01-29, 06:55:47  rise: 2023-01-29, 08:58:45
Jan 30  hbh =  2.02  set: 2023-01-30, 06:53:27  rise: 2023-01-30, 08:54:54
Jan 31  hbh =  2.00  set: 2023-01-31, 06:51:12  rise: 2023-01-31, 08:51:03
Feb 1   hbh =  1.97  set: 2023-02-01, 06:49:02  rise: 2023-02-01, 08:47:12
Feb 2   hbh =  1.94  set: 2023-02-02, 06:46:57  rise: 2023-02-02, 08:43:20
Feb 3   hbh =  1.91  set: 2023-02-03, 06:44:56  rise: 2023-02-03, 08:39:27
Feb 4   hbh =  1.88  set: 2023-02-04, 06:43:01  rise: 2023-02-04, 08:35:35
Feb 5   hbh =  1.84  set: 2023-02-05, 06:41:11  rise: 2023-02-05, 08:31:42
Feb 6   hbh =  1.81  set: 2023-02-06, 06:39:25  rise: 2023-02-06, 08:27:48
Feb 7   hbh =  1.77  set: 2023-02-07, 06:37:44  rise: 2023-02-07, 08:23:55
Feb 8   hbh =  1.73  set: 2023-02-08, 06:36:08  rise: 2023-02-08, 08:20:00
Feb 9   hbh =  1.69  set: 2023-02-09, 06:34:36  rise: 2023-02-09, 08:16:05
Feb 10  hbh =  1.65  set: 2023-02-10, 06:33:08  rise: 2023-02-10, 08:12:10
Feb 11  hbh =  1.61  set: 2023-02-11, 06:31:46  rise: 2023-02-11, 08:08:13
Feb 12  hbh =  1.56  set: 2023-02-12, 06:30:28  rise: 2023-02-12, 08:04:16
Feb 13  hbh =  1.52  set: 2023-02-13, 06:29:14  rise: 2023-02-13, 08:00:18
Feb 14  hbh =  1.47  set: 2023-02-14, 06:28:05  rise: 2023-02-14, 07:56:19
Feb 15  hbh =  1.42  set: 2023-02-15, 06:27:01  rise: 2023-02-15, 07:52:19
Feb 16  hbh =  1.37  set: 2023-02-16, 06:26:02  rise: 2023-02-16, 07:48:18
Feb 17  hbh =  1.32  set: 2023-02-17, 06:25:07  rise: 2023-02-17, 07:44:15
Feb 18  hbh =  1.26  set: 2023-02-18, 06:24:18  rise: 2023-02-18, 07:40:10
Feb 19  hbh =  1.21  set: 2023-02-19, 06:23:34  rise: 2023-02-19, 07:36:03
Feb 20  hbh =  1.15  set: 2023-02-20, 06:22:55  rise: 2023-02-20, 07:31:54
Feb 21  hbh =  1.09  set: 2023-02-21, 06:22:22  rise: 2023-02-21, 07:27:42
Feb 22  hbh =  1.03  set: 2023-02-22, 06:21:56  rise: 2023-02-22, 07:23:27
Feb 23  hbh =  0.96  set: 2023-02-23, 06:21:37  rise: 2023-02-23, 07:19:08
Feb 24  hbh =  0.89  set: 2023-02-24, 06:21:25  rise: 2023-02-24, 07:14:44
Feb 25  hbh =  0.81  set: 2023-02-25, 06:21:23  rise: 2023-02-25, 07:10:14
Feb 26  hbh =  0.73  set: 2023-02-26, 06:21:32  rise: 2023-02-26, 07:05:35
Feb 27  hbh =  0.65  set: 2023-02-27, 06:21:56  rise: 2023-02-27, 07:00:44
Feb 28  hbh =  0.55  set: 2023-02-28, 06:22:39  rise: 2023-02-28, 06:55:36
Mar 1   hbh =  0.43  set: 2023-03-01, 06:23:55  rise: 2023-03-01, 06:49:58
Mar 2   hbh =  0.28  set: 2023-03-02, 06:26:18  rise: 2023-03-02, 06:43:15
Mar 3   hbh =  0.00
Mar 4   hbh =  0.00
    ... "hours below horizon" zero until ...
Apr 1   hbh =  0.00
Apr 2   hbh =  0.00
Apr 3   hbh =  0.35  set: 2023-04-03, 05:23:34  rise: 2023-04-03, 05:44:23
Apr 4   hbh =  0.51  set: 2023-04-04, 05:16:49  rise: 2023-04-04, 05:47:43
Apr 5   hbh =  0.65  set: 2023-04-05, 05:11:08  rise: 2023-04-05, 05:50:00
Apr 6   hbh =  0.76  set: 2023-04-06, 05:05:56  rise: 2023-04-06, 05:51:49
Apr 7   hbh =  0.87  set: 2023-04-07, 05:01:02  rise: 2023-04-07, 05:53:21
Apr 8   hbh =  0.97  set: 2023-04-08, 04:56:20  rise: 2023-04-08, 05:54:43
Apr 9   hbh =  1.07  set: 2023-04-09, 04:51:47  rise: 2023-04-09, 05:55:56
Apr 10  hbh =  1.16  set: 2023-04-10, 04:47:20  rise: 2023-04-10, 05:57:03
Apr 11  hbh =  1.25  set: 2023-04-11, 04:42:58  rise: 2023-04-11, 05:58:07
Apr 12  hbh =  1.34  set: 2023-04-12, 04:38:40  rise: 2023-04-12, 05:59:07
Apr 13  hbh =  1.43  set: 2023-04-13, 04:34:26  rise: 2023-04-13, 06:00:04
Apr 14  hbh =  1.51  set: 2023-04-14, 04:30:14  rise: 2023-04-14, 06:00:59

P.S. PROBLEM NOW SOLVED I did forget to include f.step_days = 0.1 before the call to almanac.find_discrete. The data improves but leaves Jan 20 to 22 with no rise/set. Using f.step_days = 0.05 improves the data again slightly. Using f.step_days = 0.01 finally gives me the correct data! Thanks for your help!!

I include here the corrected test program (with minor fixes) in case anyone else finds it useful:

Mars_v3.zip

brandon-rhodes commented 1 year ago

@aendie — Okay, based on an experiment I did back at the end of 2021, I've spent the past few days adding a much faster rising-setting finder to Skyfield. Here's the documentation I've just added today about how to use it for sunrises and sunsets:

https://github.com/skyfielders/python-skyfield/blob/1ac697dca28c0feb23b9947360d827d9e191caeb/skyfield/documentation/almanac.rst#sunrise-and-sunset

But as you can see, it's not specific to the Sun, because it takes the target body as an argument. Though currently, in cases that aren't the Sun, you need to provide a 5th argument which is where you want the horizon to be placed; taking what I think is the USNO's standard estimate for refraction at the horizon, we get something like:

t, y = almanac.find_risings(bluffton, mars, t0, t1, -0.5667)
print(t.utc_iso(' '))
print(y)

I would be very interested to know:

If you want to try it out, you can simply install the development version of Skyfield:

pip install -U https://github.com/skyfielders/python-skyfield/archive/master.zip

If you get the chance to try it out, let me know what you find — it would be helpful to know if it works for someone else before the new feature's release!

brandon-rhodes commented 1 year ago

I am going to go ahead and close this, since we've tracked down the underlying problem, but would still be interested to hear any comments about the new rising / setting routines, so please feel free to make more comments here!