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

in #steemstem6 years ago (edited)

Pixnio.com link CC0 license

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

Sort:  

Worked perfectly for me! Btw the most accurate predictions are from JPL Horizons.

Thx for trying it out and thx for the tip. I will take a look at the JPL web site.

I couldn't get the exact times in Guide 9.1 or JPL, maybe refraction or could be the definition of sunset/sunrise. Is it the point the sun is completely below or completely above horizon, or is it when the sun's center is on the horizon. Either way great work on this!

I tried it out as well. There is a 15 minute difference between JPL and this cwlculation. I looked at JPL's parameter and they say it is for a vacuum so no diffraction effects are accounted for. I think this is the source of the differences.

I've started learning Object Oriented Programming in C++ and god Python seems so much clearer and logical!

Anyways, nice bit of code!

I'm wondering why did you start learning coding?

Fantastic use of python and it’s extreme versatility. I’m currently learning python as well as a beginner because it’s an easy introductory language.

My main goal is to conquer Swift eventually. My current endeavor is building a functional blockchain in python. It’s been quite the task. Object oriented programming is quite the abstract topic.

I’ve also always been a huge fan of complex mathematics and to put it simply, you make it look easy. Fantastic job. Keep up the outstanding work!!

Never heard of Swift until you mentioned it here. I see it is an Apple developed language. I'm mostly a PC guy myself.

Correct. Apple stepped away from objective-c in 2014 as their main language and have thus far concentrated solely on their in house language Swift. The syntax is beautiful to be honest. Very intuitive and if using x-code, storyboards and auto-layout are almost a work of magic, if you can wrap your head around it that is. Still object oriented though so a good foundation is python. Glad you took the time to look it up! I’m an Apple fan boy at heart, since 2000 when I started college.

I am really happy to see that you have received @utopian's support.
I upvoted your post because I like reading post with codes.
Congratulations!

Hi @procrastilearner!

Your post was upvoted by utopian.io in cooperation with steemstem - supporting knowledge, innovation and technological advancement on the Steem Blockchain.

Contribute to Open Source with utopian.io

Learn how to contribute on our website and join the new open source economy.

Want to chat? Join the Utopian Community on Discord https://discord.gg/h52nFrV

Coin Marketplace

STEEM 0.17
TRX 0.15
JST 0.028
BTC 62104.41
ETH 2404.22
USDT 1.00
SBD 2.49