Transliteration of "Arduino Uno and Solar Position Calculations", by David Brooks
In [1]:
import math
In [2]:
# Paper numbers vs. math.pi. Eventual difference is in the 4th decimal place.
paper_numbers = False
if paper_numbers:
PI = 3.141592654
TAU = 6.28318531
DEG_TO_RAD = 0.01745329
else:
# Using math.pi matches the Excel results from the paper
PI = math.pi
TAU = PI * 2
DEG_TO_RAD = TAU / 360
In [3]:
hour = 0
minute = 0
second = 0
day = 21
month = 6
year = 2015
zone = 5 # how many hours to subtract from UTC
In [4]:
# From the paper. Denotes a location in NJ near Philadelphia
Lat = 40 * DEG_TO_RAD
Lon = -75 * DEG_TO_RAD
In [5]:
def JulianDate(year, month, day):
if month <= 2:
year = year - 1
month = month + 12
A = int(year / 100) # implicit integer
B = 2 - A + int(A / 4)
JD_Whole = int((365.25 * (year + 4716)) + int(30.6001 * (month + 1)) + day + B - 1524)
# This look sketchy.
# is month zero-based? Then the table in the paper is misleading.
# is day zero-based?
JD_frac = (hour + minute / 60.0 + second / 3600.0) / 24.0 - 0.5
return int(JD_Whole), JD_frac
In [6]:
for hour in range(10, 25):
JD_whole, JD_frac = JulianDate(year, month, day)
T = ((JD_whole - 2451545) + JD_frac) / 36525.0
L0 = DEG_TO_RAD * (280.46645 + 36000.76983 * T) % 360.0
M = DEG_TO_RAD * (357.5291 + 35999.0503 * T) % 360.0
e = 0.016708617 - 0.000042037 * T
C = DEG_TO_RAD * ((1.9146 - 0.004847 * T) * math.sin(M) + (0.019993 - 0.000101 * T) * math.sin(2*M) + 0.00029 * math.sin(3 * M))
f = M + C
Ob1 = DEG_TO_RAD * (23.0 + 26.0/60.0 + 21.448/3600.0 - 46.815 / 3600 * T)
JDx = JD_whole - 2451545
GrHrAngle = 280.46061837 + (360 * JDx) % 360 + 0.98564736629 * JDx + 360.98564736629 * JD_frac
GrHrAngle = GrHrAngle % 360.0
L_true = (C + L0) % TAU
R = 1.000001018 * (1.0 - e*e) / (1.0 + e * math.cos(f))
RA = math.atan2(math.sin(L_true) * math.cos(Ob1), math.cos(L_true))
Decl = math.asin(math.sin(Ob1) * math.sin(L_true))
HrAngle = DEG_TO_RAD * GrHrAngle + Lon - RA
elev = math.asin(math.sin(Lat) * math.sin(Decl) + math.cos(Lat) * (math.cos(Decl) * math.cos(HrAngle)))
azimuth = PI + math.atan2(math.sin(HrAngle), math.cos(HrAngle) * math.sin(Lat) - math.tan(Decl) * math.cos(Lat))
print("{:>2}: elevation: {:>8.4f}, azimuth: {:>8.4f}".format(hour - zone, elev/DEG_TO_RAD, azimuth/DEG_TO_RAD))