Add routine that reproduces HORIZONS rise/set times, as alternative to Naval Observatory rise/set times

[ghost]

When elevation changes rise-set times have to change also due to atmospheric dip. Checked with JPL Horizons.

[brandon-rhodes]

Please provide the time and date of a specific rising or setting from HORIZONS, and a Skyfield script that produces a different prediction. Thanks!

Note that Skyfield is based on Naval Observatory code, not JPL's SPICE library, so the difference might simply be an automatic adjustment that the JPL makes that the USNO does not. So the resolution of this issue might be a documentation improvement in Skyfield explaining the difference in approaches and recommending to users which one they should use.

[ghost]

The following gives good results because dip = 0. from skyfield.api import load, Topos from skyfield import almanac body = load('de421.bsp') ts = load.timescale() corinth = Topos('37.94 N','22.93 E', elevation_m=0) # <== 0 meters is OK t0 = ts.utc(2020, 3, 15) t1 = ts.utc(2020, 3, 18) t, y = almanac.find_discrete(t0, t1, almanac.sunrise_sunset(body, corinth)) for ti, yi in zip(t, y): print(ti.utc_iso(), "Rise" if yi else "Set")

output skyfield--------------------------------------------JPL Horizons 2020-03-15T04:39:04Z Rise---------------2020-Mar-15 04:40
2020-03-15T16:35:44Z Set--------------- 2020-Mar-15 16:36 2020-03-16T04:37:33Z Rise--------------- 2020-Mar-16 04:38 2020-03-16T16:36:41Z Set--------------- 2020-Mar-16 16:37 2020-03-17T04:36:02Z Rise--------------- 2020-Mar-17 04:37 2020-03-17T16:37:37Z Set--------------- 2020-Mar-17 16:38

corinth = Topos('37.94 N','22.93 E', elevation_m=500) # not OK !!!

output skyfield--------------------------------------------JPL Horizons 2020-03-15T04:39:04Z Rise---------------2020-Mar-15 04:36 2020-03-15T16:35:44Z Set--------------- 2020-Mar-15 16:40 2020-03-16T04:37:33Z Rise--------------- 2020-Mar-16 04:34 2020-03-16T16:36:41Z Set--------------- 2020-Mar-16 16:41 2020-03-17T04:36:02Z Rise--------------- 2020-Mar-17 04:33 2020-03-17T16:37:37Z Set--------------- 2020-Mar-17 16:42

corinth = Topos('37.94 N','22.93 E', elevation_m=1500) # even worst !!!

output skyfield--------------------------------------------JPL Horizons 2020-03-15T04:39:04Z Rise-------------2020-Mar-15 04:33 2020-03-15T16:35:44Z Set--------------- 2020-Mar-15 16:43 2020-03-16T04:37:33Z Rise--------------- 2020-Mar-16 04:31 2020-03-16T16:36:41Z Set--------------- 2020-Mar-16 16:44 2020-03-17T04:36:02Z Rise--------------- 2020-Mar-17 04:30 2020-03-17T16:37:37Z Set----------------- 2020-Mar-17 16:45

[brandon-rhodes]

What setting are you using in HORIZONS for elevation? I have not found one. I see an "Altitude", but that is different than elevation. Elevation means "the height of the whole landscape above sea level", and it's presumed that the horizon moves up to the same height: sunrise on, say, the Colorado Plateau involves both the ground and the horizon being at, say, 5,000 ft above sea level. Whereas "Altitude" means "an object's height above the surrounding terrain, including above the terrain at the horizon."

Would additional documentation in the Almanac section on risings and settings have helped prevent confusion between the two concepts? We could add a small guide to computing where the horizon is if one's altitude is nonzero.

[ghost]

These are the settings @ https://ssd.jpl.nasa.gov/horizons.cgi

Ephemeris Type [change] :	OBSERVER
Target Body [change] :	Sun [Sol] [10]
Observer Location [change] :	Topocentric ( 22°55'48.0''E, 37°56'24.0''N, 500 m )
Time Span [change] :	Start=2020-03-15, Stop=2020-03-18, Step=1 m
Table Settings [change] :	QUANTITIES=4; angle format=DEG; refraction model=REFRACTED; RTS flag=TVH; object page=NO
Display/Output [change] :	default (formatted HTML)

Skyfield should change the rise-set times when elevation (or altitude at JPL) changes.

[ghost]

Let me put it differently. How can i compute the different rise-set times at locations A and B with skyfield (since elevation is not the case)?

[brandon-rhodes]

I will try to put together some sample code!

It is possible that my earlier explanation is not what is going on with HORIZONS; I will have to try coding alternatives until we find one that matches the HORIZONS data. Only that will tell us the difference between the ideas "The difference is atmospheric refraction, which changes with elevation" and "The difference is that HORIZONS thinks you are on a hill."

In either case, the key is going to be creating a routine like sunrise_sunset() in almanac.py but that uses a special value instead of the standard -0.8333 specified by the Naval Observatory. You could try cutting and pasting that sunrise_sunset() routine into your own code and trying different values there to see if you can match HORIZONS in the two cases you cite above?

[ghost]

This could help https://digitalcommons.mtu.edu/cgi/viewcontent.cgi?article=1776&context=etdr pages, 28-29

To calculate the dip due to altitude is not so difficult. Just add some coefficient*sqrt(altitude) to 0.833deg, but it is not crystal clear to me issues like refraction bellow horizon (0deg).

Any way the must read paper is this one http://articles.adsabs.harvard.edu/cgi-bin/nph-iarticle_query?1997AN....318..305W&defaultprint=YES&filetype=.pdf

[brandon-rhodes]

You can also search the Skyfield code for the word "refraction" — I think it has one refraction formula already? I think it's an option when calling altaz().

[ghost]

From your code almanac.py (lines 155-157)

Skyfield uses the same definition as the United States Naval Observatory: the Sun is up when its center is 0.8333 degrees below the horizon,

which horizon?

First think for simplicity [case location A] that there is no refraction no atmosphere at all. what happens then? We see not from 0deg to 90deg but from 0deg to lets say 92.4deg, where 2.4deg = 0.02931*sqrt(h) [Wittmann paper, previous 2nd link] is the horizon dip. Maybe it's that simple.

The sun rises not from horizon at angle 90deg (from zenith), but from horizon at angle 92.4deg!

In my opinion you could keep the same refraction since there is no big difference to rise set times (how big could the refraction difference be at 90deg or 92deg), but the code must include the change of the horizon. That makes a significant difference.

[brandon-rhodes]

In the altazimuth system, "the horizon" is defined as the great circle where the altitude is zero:

https://en.wikipedia.org/wiki/Horizontal_coordinate_system

the code must include the change of the horizon. That makes a significant difference.

But wouldn’t that break Skyfield by making its output disagree with the official sunrise and sunset proclaimed by the Naval Observatory’s Astronomical Almanac?

Have you tried adjusting the refraction or horizon to see if you can reproduce exactly the sunrise and sunset times of HORIZONS? Only once we have code that produces exactly the same time as HORIZONS can we be confident that we know exactly how their approach is different from the official one.

[ghost]

Yes if you include the horizon dip, rise-set time must be ok as far as i remember when a couple of months ago i did the calculation with my fortran code and the spice library. The spice library provides only position vectors and i had to include the horizon dip to take similar results with JPL Horizons.

But wouldn’t that break Skyfield by making its output disagree with the official sunrise and sunset proclaimed by the Naval Observatory’s Astronomical Almanac?

What if it is clear that this is not true? Everybody knows that the sun rises earlier if you are in a tall place. I would suggest you to put a new flag-variable lets say REAL rise-set time = ON or something like that. There are people out there who wish to know when the sun rises and they dont really care about the formal definitions of institutions or wikipedia.

[ghost]

P.S. If you will not make it i will stuck with my fortran numbers once again (which is not a bad thing because i like fortran). Anyway i would indeed like to give you my congratulations for your wonderful software.

[brandon-rhodes]

Everybody knows that the sun rises earlier if you are in a tall place.

Skyfield currently has no way to represent situation "A" that you illustrate — the Topos class can represent only a location with an elevation that it shares with the surrounding terrain, like Flagstaff, Arizona, or Denver, Colorado, both of which are locations more than a mile high but whose horizon terrain is also a mile high, cancelling the effect. If Skyfield did auto-adjust the definition of "sunrise" to compensate for such locations, I think it would make sunrise later, because the air would be thinner and not refract as much? Do we know whether HORIZONS makes that adjustment?

Handling situation "A", where sunrise is observed from the top of a tower, or an airplane, or an Earth satellite, could be handled in one of two ways. Either the Earth could be considered a simple sphere and basic trigonometry could provide an approximate negative altitude angle for the position of the horizon, or a better calculation could be done considering Earth to be an oblate spheroid, perhaps that of WGS84. Do you happen to know which one HORIZONS is using?

[ghost]

Here in Greece and in Europe in general there are plenty of sharp mountains (case A) and as far as i know this is not the case in the U.S (although America's landscapes are amazing!).

Anyway JPL uses a simple model for refraction (it computes refraction just for standard atmospheric conditions 10oC and 1010hPa) and yes

Either the Earth could be considered a simple sphere and basic trigonometry could provide an approximate negative altitude angle for the position of the horizon

i think this is the solution.

[brandon-rhodes]

Thanks for the answers!

It may be several weeks until I can make any progress on this issue. Once I have time to get started, it will take time to figure out exactly which formulae will produce exact agreement with HORIZONS, so that Skyfield and HORIZONS give exactly the same second for sunrise and sunset. If you happen to run across a formula that produces exact to-the-second agreement for any of the HORIZONS output you provided above, please share it here, as it would make adding the feature to Skyfield very fast.

[ghost]

1) The thinner air above the observer does not contribute to refraction (according to several papers). 2) As far as i remember i took rise-set times from JPL for several heights with and without refraction and i noticed (through interpolation) the sqrt dependence of the refracted dip.

It is indeed a complicated issue, i'll let you know if i find something interesting. Thanks for the conversation.