Can Almanac consider the elevation in Topos objects?

[dkazdan]

I would like to be able to find time of sunrise at an elevation above mean sea level (300 km up). I don't think the Almanac methods consider elevation, is there a way to do this?

Thanks...

[brandon-rhodes]

Good question!

Skyfield assumes that your local horizon is at the same elevation as you are — which is generally the correct assumption: someone watching sunrise in Denver is not at the top of a mile-high cliff looking down at a sea from which the Sun is rising. Instead, their local horizon, which is typically about ~3 miles away, is probably at nearly the exact same height they are.

Given your figure of 300 km, far higher than any elevation on Earth, let's instead use the name "altitude" — how can you measure rising and setting for an object at a very high altitude above the surface?

A first thought: if the object at 300 km altitude is an Earth satellite, you will soon be able to call a method to find out if it's in sunlight or not:

https://github.com/skyfielders/python-skyfield/pull/369

I'd expect that feature to land in a few days.

A second thought: you could use the Almanac routine risings_and_settings() but give it a more negative than usual value for horizon_degrees to account for being at a high altitude above the surface.

Could you share a few more details of the situation you're modeling? Maybe there should be a helper routine to help folks at high altitude compute the right amount of negative horizon.

[brandon-rhodes]

@dkazdan — I just saw your comment on #370!

Wow, that's an exciting project. I've only ever heard WWV from the ground — my late grandfather would tune in when showing me his ham radio, and we would set our watches to the tones. I also heard at least one leap second that way.

I have just drawn a sketch showing a satellite at 300 km altitude, with its horizon defined by a line passing through it that's tangent to the curve of the Earth. That produces a right triangle whose right angle is right there at the tangent point where the Earth's radius meets the line-of-sight from the satellite. Doing a bit of trig, I get:

from numpy import arccos
from skyfield.constants import ERAD, RAD2DEG
altitude_m = 300e3
horizon_dip = arccos(ERAD / (ERAD + altitude_m))
print(horizon_dip * RAD2DEG)

17.238916804212003

So the object at 300km altitude would "see" a horizon 17 degrees below what it would see if it were at the surface. Does that value help you move forward, or are there move moving pieces we should assemble?

[dkazdan]

Hmm, drawings like these, that my dear wife and I were just drawing? Please do look in on us at www.hamsci.org I'll write more this weekend.  73, David CC:Kristina Collins, KD8OXT, PhD candidate in this stuff. 

[brandon-rhodes]

@dkazdan — Alas, the comment that GitHub created from your reply failed to include the drawings! Feel free to email me at brandon@rhodesmill.org if you want to send them over a less noisy channel.

[dkazdan]

OK, thanks.  I'm sorry to be thick here--but do I now need to write a rootfinding routine to find when the sun is that far below the horizon from the corresponding Topos point? I see the appropriate function in the PyEphem but not in Skyfield.        David

[brandon-rhodes]

OK, thanks. I'm sorry to be thick here--but do I now need to write a rootfinding routine to find when the sun is that far below the horizon from the corresponding Topos point? I see the appropriate function in the PyEphem but not in Skyfield.

I think you can use the Almanac routines in Skyfield.

Would the F-layer above Ohio be illuminated by the Sun about 1 hour 15 minutes before a ground observer? I had not expected the gap to be that big, but:

# How far down is the horizon if you are floating in the F-layer?

from numpy import arccos
from skyfield.constants import ERAD, RAD2DEG
altitude_m = 300e3
horizon_degrees = - arccos(ERAD / (ERAD + altitude_m)) * RAD2DEG

# When does the sun rise for an earth-bound observer, and for an
# observer up in the F-layer?

from skyfield import api
from skyfield import almanac
ts = api.load.timescale()

eph = api.load('de421.bsp')
sun = eph['sun']
topos = api.Topos('40.8939 N', '83.8917 W')

EDT = -4
t0 = ts.utc(2020, 5, 4, -EDT)
t1 = ts.utc(2020, 5, 5, -EDT)

f = almanac.risings_and_settings(eph, sun, topos)
t, y = almanac.find_discrete(t0, t1, f)
for ti, yi in zip(t, y):
    if yi:
        print(ti.utc_jpl(), '  Sunrise')

f = almanac.risings_and_settings(eph, sun, topos, horizon_degrees)
t, y = almanac.find_discrete(t0, t1, f)
for ti, yi in zip(t, y):
    if yi:
        print(ti.utc_jpl(), '  F-layer sunrise')

produces:

A.D. 2020-May-04 10:30:58.6079 UT   Sunrise
A.D. 2020-May-04 08:47:31.6571 UT   F-layer sunrise

Does that answer look reasonable, given your understanding of the geometry? (Those UT times are 4:47am and 6:30am EDT. I'm not usually up at 4:47am! The F-layer is woken early.)

[dkazdan]

Thank you! Yes, that's exactly what I need, and that time difference is roughly what I expected.  Solving the altitude triangle for the distance to the horizon from 300 kilometers altitude, I got d=ERAD*tan(arccos(ERAD/(ERAD+elevation_m))), or about 1500 kilometers.  The time difference to sunrise will depend on the latitude, but TLAR: that looks about right.  Our immediate interest is around 40 degrees north latitude. If you're interested, we can describe at greater length what the implications of that are for our WWV measurements.  The short version is that the Doppler shift increase in the beacon's received frequency for receivers east of the beacon begins at about that time, and seeing the day to day variation has been a good observation for us.  This computation gives much predictive power to our models. Care to join us for our weekly Zoom conference this Thursday?  Our group from CWRU and others is in the morning; the larger HamSCI call is in the afternoon.

           73 (amateur radio operator's "best regards"),           David

[brandon-rhodes]

I’m glad we might have found you a workable approach! I’ll see what my schedule looks like and email you about Thursday.