Daylight Algorithm

daylight.py

def get_julian_day(date):
    a = (14 - date.month) // 12
    y = date.year + 4800 - a
    m = date.month + 12 * a - 3
    julian_date = date.day + (153 * m + 2) // 5 + 365 * y + y // 4 - y // 100 + y // 400 - 32045

    return julian_date - 2451545 + .0008

def daylight(date, lat, lon):
    """ Returns daylight hours for a single lat/lon point """

    julian_day = get_julian_day(date)
    solar_noon = julian_day - lon//360
    solar_anomaly = (357.5291 + 0.98560028*solar_noon) % 360
    equation_of_center = (
        1.9148*math.sin(math.radians(solar_anomaly)) +
        0.0200*math.sin(math.radians(2*solar_anomaly)) +
        0.0003*math.sin(math.radians(3*solar_anomaly))
    )
    ecliptic_longitude = (solar_anomaly + equation_of_center + 180 + 102.9372) % 360
    solar_transit = (
        2451545.5 + solar_noon + 0.0053*math.sin(math.radians(solar_anomaly)) -
        0.0069*math.sin(math.radians(2*ecliptic_longitude))
    )
    declination = math.asin(math.sin(math.radians(ecliptic_longitude))*math.sin(math.radians(23.44)))
    hour_angle = math.acos(
        (math.sin(math.radians(-.83)) - math.sin(math.radians(lat))*math.sin(declination)) /
        (math.cos(math.radians(lat))*math.cos(declination))
    )

    sunrise = solar_transit - math.degrees(hour_angle)/360
    sunset = solar_transit + math.degrees(hour_angle)/360

    return (sunset - sunrise) * 24
  
  def daylight_array(date, lat, lon):
    """ Returns daylight hours for arrays of lat/lon points """

    julian_day = get_julian_day(date)

    lat_arr = numpy.tile(lat.reshape(len(lat), 1), (1, len(lon)))
    lon_arr = numpy.tile(lon, (len(lat), 1))

    # solar_noon = julian_day - lon//360
    solar_noon = lon_arr
    del lon_arr
    solar_noon //= 360
    solar_noon -= julian_day
    solar_noon *= -1

    # solar_anomaly = (357.5291 + 0.98560028 * solar_noon) % 360
    solar_anomaly = solar_noon * 0.98560028
    solar_anomaly += 357.5291
    solar_anomaly %= 360

    # equation_of_center = (
    #     1.9148 * math.sin(math.radians(solar_anomaly)) +
    #     0.0200 * math.sin(math.radians(2 * solar_anomaly)) +
    #     0.0003 * math.sin(math.radians(3 * solar_anomaly))
    # )
    equation_of_center = numpy.radians(solar_anomaly)
    numpy.sin(equation_of_center, equation_of_center)
    equation_of_center *= 1.9148
    equation_of_center_2 = solar_anomaly * 2
    numpy.radians(equation_of_center_2, equation_of_center_2)
    numpy.sin(equation_of_center_2, equation_of_center_2)
    equation_of_center_2 *= 0.0200
    equation_of_center += equation_of_center_2
    equation_of_center_2 = solar_anomaly * 3
    numpy.radians(equation_of_center_2, equation_of_center_2)
    numpy.sin(equation_of_center_2, equation_of_center_2)
    equation_of_center_2 *= 0.0003
    equation_of_center += equation_of_center_2
    del equation_of_center_2

    # ecliptic_longitude = (solar_anomaly + equation_of_center + 180 + 102.9372) % 360
    ecliptic_longitude = equation_of_center
    del equation_of_center
    ecliptic_longitude += solar_anomaly
    ecliptic_longitude += 282.9372  # 180 + 102.9372
    ecliptic_longitude %= 360

    # solar_transit = (
    #     2451545.5 + solar_noon + 0.0053*math.sin(math.radians(solar_anomaly)) -
    #     0.0069*math.sin(math.radians(2*ecliptic_longitude))
    # )
    solar_transit = solar_noon
    del solar_noon
    solar_transit += 2451545.5
    numpy.radians(solar_anomaly, solar_anomaly)
    numpy.sin(solar_anomaly, solar_anomaly)
    solar_anomaly *= 0.0053
    solar_transit += solar_anomaly
    del solar_anomaly
    solar_transit_2 = ecliptic_longitude * 2
    numpy.radians(solar_transit_2, solar_transit_2)
    numpy.sin(solar_transit_2, solar_transit_2)
    solar_transit_2 *= 0.0069
    solar_transit -= solar_transit_2
    del solar_transit_2

    # declination = math.asin(math.sin(math.radians(ecliptic_longitude))*math.sin(math.radians(23.44)))
    declination = ecliptic_longitude
    del ecliptic_longitude
    numpy.radians(declination, declination)
    numpy.sin(declination, declination)
    declination *= math.sin(math.radians(23.44))
    numpy.arcsin(declination, declination)

    # hour_angle = math.acos(
    #     (math.sin(math.radians(-.83)) - math.sin(math.radians(lat)) * math.sin(declination)) /
    #     math.cos(math.radians(lat)) * math.cos(declination)
    # )
    hour_angle = numpy.radians(lat_arr)
    numpy.sin(hour_angle, hour_angle)
    hour_angle *= numpy.sin(declination)
    hour_angle -= math.sin(math.radians(-.83))
    hour_angle *= -1
    numpy.radians(lat_arr, lat_arr)
    numpy.cos(lat_arr, lat_arr)
    numpy.cos(declination, declination)
    lat_arr *= declination
    del declination
    hour_angle /= lat_arr
    del lat_arr
    numpy.arccos(hour_angle, hour_angle)

    # sunrise = solar_transit - math.degrees(hour_angle) / 360
    # sunset = solar_transit + math.degrees(hour_angle) / 360
    # return (sunset - sunrise) * 24
    numpy.degrees(hour_angle, hour_angle)
    hour_angle /= 360
    solar_transit = solar_transit.astype('float64')
    days = (solar_transit + hour_angle) - (solar_transit - hour_angle)
    days *= 24
    return days