Transliteration of "Arduino Uno and Solar Position Calculations", by David Brooks

http://www.instesre.org/ArduinoUnoSolarCalculations.pdf

``````

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
else:
# Using math.pi matches the Excel results from the paper
PI = math.pi
TAU = PI * 2

``````
``````

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

``````
``````

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))

``````
``````

5: elevation:   3.9347, azimuth:  62.4369
6: elevation:  14.4962, azimuth:  71.3769
7: elevation:  25.6184, azimuth:  80.0108
8: elevation:  37.0417, azimuth:  89.0132
9: elevation:  48.4861, azimuth:  99.4738
10: elevation:  59.4964, azimuth: 113.6657
11: elevation:  68.9263, azimuth: 137.1846
12: elevation:  73.4329, azimuth: 178.5533
13: elevation:  69.3859, azimuth: 220.8900
14: elevation:  60.1246, azimuth: 245.2711
15: elevation:  49.1647, azimuth: 259.8137
16: elevation:  37.7304, azimuth: 270.4153
17: elevation:  26.2972, azimuth: 279.4679
18: elevation:  15.1498, azimuth: 288.1025
19: elevation:   4.5468, azimuth: 297.0085

``````