0

At Sunrise sun is supposed to be at the horizon (ie 0 degree altitude at east). but if I check ecliptic longitude of:

  1. the sun at sunrise and
  2. from_altaz() of zero degrees east

I get different values for both. But is it not supposed to be the same?

from skyfield import almanac, api
from pytz import timezone
from datetime import datetime, timedelta


ts = api.load.timescale(builtin=True)
eph = api.load('de421.bsp')
earth = eph['EARTH']
sun = eph['SUN']
tz = timezone('UTC')

# 2019-12-23 01:23:58.273000+00:00
testTime = datetime(year=2019, month=12, day=23, hour=1,
                    minute=23, second=58)
testTime = tz.localize(testTime)

lat = '33.775867N'
lon = '84.39733E'
observer = api.Topos(lat, lon)

t0 = ts.utc(testTime)
t1 = ts.utc(testTime + timedelta(days=1))
t, y = almanac.find_discrete(t0, t1, almanac.sunrise_sunset(eph, observer))
for i, j in zip(t, y):
    if j:
        print('Sunrise : ', i.utc_datetime().astimezone(tz), j)

observer = earth + api.Topos(lat, lon)

zeroDegreeEast = observer.at(t0).from_altaz(alt_degrees=0, az_degrees=90)
_, lonAtHorizon, _ = zeroDegreeEast.ecliptic_latlon(epoch=t0)

sunAtSunRise = observer.at(t0).observe(sun).apparent()
_, lonSun, _ = sunAtSunRise.ecliptic_latlon(epoch=t0)
print('TestDate: ', testTime.isoformat())
print(lonAtHorizon)
print(lonSun)
Sunrise :  2019-12-23 01:23:58.273000+00:00 True
TestDate:  2019-12-23T01:23:58+00:00
288deg 05' 48.3"
270deg 53' 48.0"

what am I missing?

Is there anything I could do to get the same values?

AviD1511
  • 133
  • 7

1 Answers1

1

There are two additional effects that your script does not consider.

  1. By definition, sunrise is when the center of the sun is 0.8333 degrees below the horizon, because it's half a degree wide and additionally its image is refracted upward and is visible ahead of schedule because of the atmosphere.
  2. The sun never rises directly to the east, because of the seasons.

You can ask for the Sun's position on the horizon at the moment of sunrise that you've computed to account for the second effect:

alt, az, distance = ((observer + eph['earth']).at(t[0])
                     .observe(eph['sun']).apparent().altaz())
print(alt, az)
print(az.degrees)

Which prints:

-00deg 49' 59.9" 117deg 57' 17.6"
117.95487601490409

or about 27° south of directly east. You can then account for the first effect by adjusting the altaz you are computing:

zeroDegreeEast = observer.at(t0).from_altaz(alt_degrees=-0.8333, az_degrees=117.954)

The result should very closely match what you expect.

Brandon Rhodes
  • 83,755
  • 16
  • 106
  • 147