Killing Time With Recreational Math - Calculate Sunrise and Sunset Times Using Python

A little over 2 months ago I created a post called Kill Time With Recreational Math: Calculate Sunset and Sunrise Times.

In this post I converted the sunrise equation from Wikipedia into an Excel spreadsheet. That was quite the job but it was fun and well worth it as I now have a spreadsheet that does a nifty little calc for me. The job at the time was tricky because converting from degrees to radians and back again was not always obvious but I got it to work. Time offsets between Julian dates, our normal dates and the time differences between Greenwich UK and other locations was also tricky.

Learning Python has been on my to do list for a long time now and the best way to learn is to do. So, in today's post I have converted that spreadsheet into Python code.

The explanation of the variables used in the calculation can be found in the older post as well as on the Wikipedia page itself.

As well, a great free Python interpreter can be found at Python.org.

It is simple to install and the Integrated Development and Learning Environment (IDLE) is not that bad.

The Sunrise/Sunset Source Code

Here is the code, just copy and paste it into your favourite Python interpreter and try it out.

# **************************************************************************
# This code is released by Procrastilearner under the CC BY-SA 4.0 license.
# 
# Source for the sunrise calculation:
#     https://en.wikipedia.org/wiki/Sunrise_equation
# **************************************************************************
import time
import math
import datetime
import math
# *****************************************
# Some sample locations
# Toronto Ontario Canada
latitude_deg =43.65
longitude_deg = -79.38
timezone = -4.0 #Daylight Savings Time is in effect, this would be -5 for winter time

#Whitehorse Yukon Territories Canada
#latitude_deg =60.7
#longitude_deg = -135.1
#timezone = -7.0 #Daylight Savings Time is in effect, this would be -8 for winter time

#Paris France
#latitude_deg =48.85
#longitude_deg = 2.35
#timezone = 2.0 

#Hong Kong PRC
#latitude_deg =22.32
#longitude_deg =114.1
#timezone = 8.0 

#Perth Australia
#latitude_deg =-31.9
#longitude_deg =115.9
#timezone = 8.0 
# *****************************************

def date_to_jd(year,month,day):
    # Convert a date to Julian Day.
    # Algorithm from 'Practical Astronomy with your Calculator or Spreadsheet', 
    # 4th ed., Duffet-Smith and Zwart, 2011.
    # This function extracted from https://gist.github.com/jiffyclub/1294443
    if month == 1 or month == 2:
        yearp = year - 1
        monthp = month + 12
    else:
        yearp = year
        monthp = month
    # this checks where we are in relation to October 15, 1582, the beginning
    # of the Gregorian calendar.
    if ((year < 1582) or
        (year == 1582 and month < 10) or
        (year == 1582 and month == 10 and day < 15)):
        # before start of Gregorian calendar
        B = 0
    else:
        # after start of Gregorian calendar
        A = math.trunc(yearp / 100.)
        B = 2 - A + math.trunc(A / 4.)

    if yearp < 0:
        C = math.trunc((365.25 * yearp) - 0.75)
    else:
        C = math.trunc(365.25 * yearp)
    D = math.trunc(30.6001 * (monthp + 1))
    jd = B + C + D + day + 1720994.5
    return jd    
# end of date_to_jd    

pi=3.14159265359

latitude_radians = math.radians(latitude_deg)
longitude__radians = math.radians(longitude_deg)

jd2000 = 2451545 #the julian date for Jan 1 2000 at noon

currentDT = datetime.datetime.now()
current_year = currentDT.year
current_month = currentDT.month
current_day = currentDT.day
current_hour = currentDT.hour

jd_now = date_to_jd(current_year,current_month,current_day)

n = jd_now - jd2000 + 0.0008

jstar = n - longitude_deg/360

M_deg = (357.5291 + 0.98560028 * jstar)%360
M = M_deg * pi/180

C = 1.9148 * math.sin(M) + 0.0200 * math.sin(2*M) + 0.0003 * math.sin(3*M)

lamda_deg = math.fmod(M_deg + C + 180 + 102.9372,360)

lamda = lamda_deg * pi/180

Jtransit = 2451545.5 + jstar + 0.0053 * math.sin(M) - 0.0069 * math.sin(2*lamda)

earth_tilt_deg = 23.44
earth_tilt_rad = math.radians(earth_tilt_deg)

sin_delta = math.sin(lamda) * math.sin(earth_tilt_rad)
angle_delta = math.asin(sin_delta)

sun_disc_deg =  -0.83
sun_disc_rad = math.radians(sun_disc_deg)

cos_omega = (math.sin(sun_disc_rad) - math.sin(latitude_radians) * math.sin(angle_delta))/(math.cos(latitude_radians) * math.cos(angle_delta))

omega_radians = math.acos(cos_omega)
omega_degrees = math.degrees(omega_radians)

#Output section
print("------------------------------")
print("Today's date is " + currentDT.strftime("%Y-%m-%d"))
print("------------------------------")
#("%Y-%m-%d %H:%M")

print("Latitude =  " + str(latitude_deg))
print("Longitude = " + str(longitude_deg))
print("Timezone =  " + str(timezone))
print("------------------------------")

Jrise = Jtransit - omega_degrees/360
numdays = Jrise - jd2000
numdays =  numdays + 0.5 #offset because Julian dates start at noon
numdays =  numdays + timezone/24 #offset for time zone
sunrise = datetime.datetime(2000, 1, 1) + datetime.timedelta(numdays)
print("Sunrise is at " + sunrise.strftime("%H:%M"))

Jset = Jtransit + omega_degrees/360
numdays = Jset - jd2000
numdays =  numdays + 0.5 #offset because Julian dates start at noon
numdays =  numdays + timezone/24 #offset for time zone
sunset = datetime.datetime(2000, 1, 1) + datetime.timedelta(numdays)
print("Sunset is at  " + sunset.strftime("%H:%M"))
print("------------------------------")

Code Output And Code Validation

The code has been tested for a few different locations and compared against the calculated provided by Google. The answers from this code compare very well with the Google calculations.

How To Find Sunrise and Sunset Times:

Just go to Google and enter the following search term: "Toronto sunrise sunset" and the web site will give the current day's times (or whatever city you happen to live near).

How To Find Your Latitude and Longitude:

You most likely live somewhere else so to get your own sunrise and sunset times just enter your latitude in "latitude_deg", your longitude in "longitude_deg" and your time zone offset in "timezone". Go to Google Maps and left click on your location.

Your Time Zone Offset:

If you live west of Greenwich UK then the time zone offset will be negative. If you live east of London UK then the time zone offset will be positive.

Some Sample Locations Are Provided:

I have provided a few sample locations near the start of the code, these are the lines that are commented out. Just remove the comment hash character to calculate the times for those locations.

Caveat:

If you go far enough north the sun won't set in the summer and it won't rise in the winter. Vice versa if you go far enough south. In these locations the code may give you wacky answers or just choke and die. So, no guarantees if you live in Tuktoyaktuk Northwest Territories Canada.

Have fun, use and adapt as you see fit.

Thank you for reading my post.

Post Sources

[1] Sunrise equation
[2] Python.org