search

skyfielders / python-skyfield

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

How to compute culmination on polar night/day ? #243

Closed int8 closed 4 years ago

int8 commented 5 years ago

Hi there, I'm looking for a function that will point to the brightest/darkest time of a polar night / polar day - how do I get it with skyfield?

looking for culmination from here:

https://www.suncalc.org/#/69.9604,21.1157,5/2019.07.03/17:43/1/3

in ephem that is next_transit() function <= please take a look here https://chrisramsay.co.uk/posts/2017/03/fun-with-the-sun-and-pyephem/#Rising,-Transit-&-Setting

brandon-rhodes commented 5 years ago

I don't yet have a general-purpose routine that I'm happy with for finding a maxima or minima, because if you look very very closely at the floating point numbers generated for times very close together, there's bits of jaggedness that means that a moment like solar noon in fact looks like several little local minima, even though really there's a smooth curve behind the numbers.

For example, here's an attempt to use what I have so far to compute solar noon for a location 1° from the pole:


from skyfield.api import load, Topos

ts = load.timescale()

planets = load('de421.bsp')
earth = planets['earth']
sun = planets['sun']

t0 = ts.utc(2019, 3, 10)
t1 = ts.utc(2019, 3, 20)

place = Topos('89 N', '0 E') + earth

def f(t):
    alt, az, distance = place.at(t).observe(sun).apparent().altaz()
    return alt.degrees

f.rough_period = 1.0

from skyfield.almanac import _find_maxima
times, maxima = _find_maxima(t0, t1, f)

for t, altitude in zip(times, maxima):
    print('At {} the sun is {:.2f} degrees high'.format(t.utc_jpl(), altitude))

Its output is:

At A.D. 2019-Mar-10 12:24:40.6406 UT the sun is -3.11 degrees high
At A.D. 2019-Mar-10 12:24:40.6418 UT the sun is -3.11 degrees high
At A.D. 2019-Mar-10 12:24:40.6423 UT the sun is -3.11 degrees high
At A.D. 2019-Mar-11 12:24:26.7786 UT the sun is -2.71 degrees high
At A.D. 2019-Mar-11 12:24:26.7791 UT the sun is -2.71 degrees high
At A.D. 2019-Mar-13 12:23:57.3623 UT the sun is -1.93 degrees high
At A.D. 2019-Mar-14 12:23:41.8461 UT the sun is -1.53 degrees high
At A.D. 2019-Mar-14 12:23:41.8467 UT the sun is -1.53 degrees high
At A.D. 2019-Mar-15 12:23:25.8202 UT the sun is -1.14 degrees high
At A.D. 2019-Mar-15 12:23:25.8213 UT the sun is -1.14 degrees high
At A.D. 2019-Mar-16 12:23:09.3082 UT the sun is -0.74 degrees high
At A.D. 2019-Mar-16 12:23:09.3091 UT the sun is -0.74 degrees high
At A.D. 2019-Mar-17 12:22:52.3370 UT the sun is -0.35 degrees high
At A.D. 2019-Mar-18 12:22:34.9363 UT the sun is 0.05 degrees high
At A.D. 2019-Mar-18 12:22:34.9368 UT the sun is 0.05 degrees high

On some days like 2019-Mar-17 it sees a unique moment when the sun is highest. But on other days it gets several.

If you don't need microsecond precision you could do something like:

one_second = 1.0 / 24 / 60 / 60
times, maxima = _find_maxima(t0, t1, f, epsilon=one_second)

and it would become statistically far less likely that the routine would find two noons per day. But it would still be possible to encounter it.

Maybe it would make sense to post-filter the results to throw out all but the first of several points that are close together?

JoshPaterson commented 5 years ago

You can also try the function culminations from PR #211. Here's an example:

from skyfield.api import load, Topos
from skyfield.contrib.almanac2.almanac2 import culminations

ts = load.timescale()

planets = load('de421.bsp')
earth = planets['earth']
sun = planets['sun']

t0 = ts.utc(2019, 3, 10)
t1 = ts.utc(2019, 3, 20)

place = earth + Topos('89 N', '0 E')

all_culminations, kind = culminations(place, sun, t0, t1)
times = all_culminations[kind=='upper']

maxima = place.at(times).observe(sun).apparent().altaz()[0].degrees

for t, altitude in zip(times, maxima):
    print('At {} the sun is {:.2f} degrees high'.format(t.utc_jpl(), altitude))

The output is:

At A.D. 2019-Mar-10 12:24:40.6418 UT the sun is -3.11 degrees high
At A.D. 2019-Mar-11 12:24:26.7782 UT the sun is -2.71 degrees high
At A.D. 2019-Mar-12 12:24:12.3459 UT the sun is -2.32 degrees high
At A.D. 2019-Mar-13 12:23:57.3621 UT the sun is -1.93 degrees high
At A.D. 2019-Mar-14 12:23:41.8462 UT the sun is -1.53 degrees high
At A.D. 2019-Mar-15 12:23:25.8207 UT the sun is -1.14 degrees high
At A.D. 2019-Mar-16 12:23:09.3084 UT the sun is -0.74 degrees high
At A.D. 2019-Mar-17 12:22:52.3371 UT the sun is -0.35 degrees high
At A.D. 2019-Mar-18 12:22:34.9361 UT the sun is 0.05 degrees high
At A.D. 2019-Mar-19 12:22:17.1402 UT the sun is 0.44 degrees high

These times are within a half millisecond of the maximum altitude.

If you end up giving that a try I would welcome any feedback on the PR!

brandon-rhodes commented 4 years ago

@JoshPaterson — I do intend sometime this year to figure out what to do with the directed search algorithm you implemented (complete with tests!) in your pull request. Now that the more general maxima function is in good condition, I can dedicate a weekend sometime to trying out both schemes in all kinds of circumstances to see if the extra complexity of the directed solver is a big enough win over the more general finders to make me confident in maintaining it forevermore. For situations where a directed solver is perfect, like sunrise and sunset, I would have to weigh telling users to use the directed solver explicitly, versus having the various almanac functions have a hidden attribute that tells a top-level solving routine which of several lower-level ones to use. Anyway, thanks in the meantime for telling users about the alternative you’ve offered!

@int8 — You can now accomplish this using Skyfield's maxima finder, and it no longer has the problem of a jagged curve making it think that there were two culminations. Thus:

from skyfield.api import load, Topos
from skyfield.searchlib import find_maxima

ts = load.timescale()

planets = load('de421.bsp')
earth = planets['earth']
sun = planets['sun']

t0 = ts.utc(2019, 3, 10)
t1 = ts.utc(2019, 3, 20)

place = Topos('89 N', '0 E') + earth

def f(t):
    alt, az, distance = place.at(t).observe(sun).apparent().altaz()
    return alt.degrees

f.rough_period = 1.0

times, maxima = find_maxima(t0, t1, f, 1.0 / 24.0 / 3600.0, 12)

for t, altitude in zip(times, maxima):
    print('At {} the sun is {:.2f} degrees high'.format(t.utc_jpl(), altitude))

Produces:

At A.D. 2019-Mar-10 12:24:40.6535 UT the sun is -3.11 degrees high
At A.D. 2019-Mar-11 12:24:26.9977 UT the sun is -2.71 degrees high
At A.D. 2019-Mar-12 12:24:12.5942 UT the sun is -2.32 degrees high
At A.D. 2019-Mar-13 12:23:57.3778 UT the sun is -1.93 degrees high
At A.D. 2019-Mar-14 12:23:41.9013 UT the sun is -1.53 degrees high
At A.D. 2019-Mar-15 12:23:26.0021 UT the sun is -1.14 degrees high
At A.D. 2019-Mar-16 12:23:09.4851 UT the sun is -0.74 degrees high
At A.D. 2019-Mar-17 12:22:52.3179 UT the sun is -0.35 degrees high
At A.D. 2019-Mar-18 12:22:35.1181 UT the sun is 0.05 degrees high
At A.D. 2019-Mar-19 12:22:17.1706 UT the sun is 0.44 degrees high

Enjoy! :)