test.py

test.py

import numpy
import sys
import datetime

def hours_per_day(latitude, day_of_year):
    # From https://www.quora.com/How-can-you-calculate-the-length-of-the-day-on-Earth-at-a-given-latitude-on-a-given-date-of-the-year
    latitude = numpy.radians(latitude)
    # compute declination using approx from
    # https://en.wikipedia.org/wiki/Position_of_the_Sun
    # XXX check that day_of_year is correcly 0 or 1 referenced
    delta = numpy.radians(-23.44 * numpy.cos(numpy.radians((360 / 365) * (day_of_year + 10))))
    cosh = -numpy.tan(latitude) * numpy.tan(delta)
    fill = (1 + numpy.sign(cosh)) * 90 # 0 or 180 depending on sign of cosh
    h = numpy.where(abs(cosh) > 1, fill, numpy.arccos(cosh))
    return 2 * numpy.degrees(h) / 15

print hours_per_day(55.7558, 110)