Local solar time function from UTC and longitude

Question

I'm wondering if there is a python function/module that calculates the local time after midnight (or local solar time) given the UTC time and longitude? It doesn't need to take into account daylight saving time.

Thanks in advance.

Answer 1

I took a look at Jean Meeus' Astronomical Algorithms. I think you might be asking for local hour angle, which can be expressed in time (0-24hr), degrees (0-360) or radians (0-2pi).

I'm guessing you can do this with ephem. But just for the heck of it, here's some python:

#!/usr/local/bin/python

import sys
from datetime import datetime, time, timedelta
from math import pi, sin, cos, atan2, asin

# helpful constant
DEG_TO_RAD = pi / 180

# hardcode difference between Dynamical Time and Universal Time
# delta_T = TD - UT
# This comes from IERS Bulletin A
# ftp://maia.usno.navy.mil/ser7/ser7.dat
DELTA = 35.0

def coords(yr, mon, day):
    # @input year (int)
    # @input month (int)
    # @input day (float)
    # @output right ascention, in radians (float)
    # @output declination, in radians (float)

    # get julian day (AA ch7)
    day += DELTA / 60 / 60 / 24 # use dynamical time
    if mon <= 2:
        yr -= 1
        mon += 12
    a = yr / 100
    b = 2 - a + a / 4
    jd = int(365.25 * (yr + 4716)) + int(30.6 * (mon + 1)) + day + b - 1524.5

    # get sidereal time at greenwich (AA ch12)
    t = (jd - 2451545.0) / 36525

    # Calculate mean equinox of date (degrees)
    l = 280.46646 + 36000.76983 * t + 0.0003032 * t**2
    while (l > 360):
        l -= 360
    while (l < 0):
        l += 360

    # Calculate mean anomoly of sun (degrees)
    m = 357.52911 + 35999.05029 * t - 0.0001537 * t**2

    # Calculate eccentricity of Earth's orbit
    e = 0.016708634 - 0.000042037 * t - 0.0000001267 * t**2

    # Calculate sun's equation of center (degrees)
    c = (1.914602 - 0.004817 * t - .000014 * t**2) * sin(m * DEG_TO_RAD) \
        + (0.019993 - .000101 * t) * sin(2 * m * DEG_TO_RAD) \
        + 0.000289 * sin(3 * m * DEG_TO_RAD)

    # Calculate the sun's radius vector (AU)
    o = l + c # sun's true longitude (degrees)
    v = m + c # sun's true anomoly (degrees)

    r = (1.000001018 * (1 - e**2)) / (1 + e * cos(v * DEG_TO_RAD))

    # Calculate right ascension & declination
    seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * 0.001813))
    e0 = 23 + (26 + (seconds / 60)) / 60

    ra = atan2(cos(e0 * DEG_TO_RAD) * sin(o * DEG_TO_RAD), cos(o * DEG_TO_RAD)) # (radians)
    decl = asin(sin(e0 * DEG_TO_RAD) * sin(o * DEG_TO_RAD)) # (radians)

    return ra, decl

def hour_angle(dt, longit):
    # @input UTC time (datetime)
    # @input longitude (float, negative west of Greenwich)
    # @output hour angle, in degrees (float)

    # get gregorian time including fractional day
    y = dt.year
    m = dt.month
    d = dt.day + ((dt.second / 60.0 + dt.minute) / 60 + dt.hour) / 24.0 

    # get right ascention
    ra, _ = coords(y, m, d)

    # get julian day (AA ch7)
    if m <= 2:
        y -= 1
        m += 12
    a = y / 100
    b = 2 - a + a / 4
    jd = int(365.25 * (y + 4716)) + int(30.6 * (m + 1)) + d + b - 1524.5

    # get sidereal time at greenwich (AA ch12)
    t = (jd - 2451545.0) / 36525
    theta = 280.46061837 + 360.98564736629 * (jd - 2451545) \
            + .000387933 * t**2 - t**3 / 38710000

    # hour angle (AA ch13)
    ha = (theta + longit - ra / DEG_TO_RAD) % 360

    return ha

def main():
    if len(sys.argv) != 4:
        print 'Usage: hour_angle.py [YYYY/MM/DD] [HH:MM:SS] [longitude]'
        sys.exit()
    else:
        dt = datetime.strptime(sys.argv[1] + ' ' + sys.argv[2], '%Y/%m/%d %H:%M:%S')
        longit = float(sys.argv[3])
    ha = hour_angle(dt, longit)
    # convert hour angle to timedelta from noon
    days = ha / 360
    if days > 0.5:
        days -= 0.5
    td = timedelta(days=days)
    # make solar time
    solar_time = datetime.combine(dt.date(), time(12)) + td
    print solar_time

if __name__ == '__main__':
    main()