System File 5: A python script to calculate solar position

 A python script to calculate solar position without use of ephemeris library.
The script is imperfect. When the sun goes down, the position coordinates become imaginary numbers. Imaginary numbers are part of high level mathematics and not necessary for this project but are a default numerical type when certain conditions are present.
I may return to using an ephemeris library if the errors become unacceptable.

CODE
#NOAA Solar Calculations spreadsheet, Day

import math
import datetime
pi=3.14159
r=360/2/pi
B3=36.3 # Lat
B4=-82.2 # Lon
B5=TZ=-5 # time zone
#B not used
B7=11/26/2017
Y=datetime.datetime.now().year
M=datetime.datetime.now().month
D=datetime.datetime.now().day
H2="Blank"
# UTC-5=EDT
UTclock=datetime.datetime.utcnow().strftime("%H:%M:%S %Z")
LTclock=datetime.datetime.now().strftime("%H:%M:%S")
(hh, mm, ss)=UTclock.split(':')
UT=(int(hh)*3600+int(mm)*60+int(ss))/3600
LT=UT+TZ
h=datetime.datetime.now().hour
m=datetime.datetime.now().minute
s=datetime.datetime.now().second
if h > 12:
h=h-12
time=(int(h)*3600+int(m)*60+int(s))/3600
#Converting string into datetime
from datetime import datetime
datetime_object = datetime.strptime('Jun 1 2005  1:33PM', '%b %d %Y %I:%M%p')
# forward and reverse transformation
JD=(367 * Y ) - (int(7 * (Y + int((M + 9) / 12)) / 4)) + (int(275 * M / 9) + D + 1721013.5) + (UT / 24)
UTx=24 * ( JD + (int(7 * (Y + int((M + 9) / 12)) / 4)) - (int(275 * M / 9) + D + 1721013.5) - (367 * Y ) )
D2=int(JD)
E2=LT/24
F2=JD
#G=Julian Century
G2=(F2-2451545)/36525
#H=blank
#I=Geom Mean Long Sun#I2(1)=math.mod(280.46646+G2*(36000.76983+G2*0.0003032),360)
I2a=(280.46646+G2*(36000.76983+G2*0.0003032))
I2b=360
I2=I2a-int(I2a/I2b)*I2b
#J=Geometric Mean Anomaly of the Sun
J2=357.52911+G2*(35999.05029-0.0001537*G2)
#K=Eccentricity of Earth orbit
K2=0.016708634-G2*(0.000042037+0.0000001267*G2)
L2=math.sin(J2/r)*(1.914602-G2*(0.004817+0.000014*G2))+math.sin(2*J2/r)*(0.019993-0.000101*G2)+math.sin(3*J2/r)*0.000289
#M=Sun True Long degrees
M2=I2+L2
#N=Sun True Anomaly degrees
N2=J2+L2
#O=Sun Radius Vector, A.U.
O2=(1.000001018*(1-K2*K2))/(1+K2*math.cos(N2))
#P=Sun Apparent Long
P2=M2-0.00569-0.00478*math.sin(125.04-1934.136*G2)
#Q=Mean Oblique Ecliptic degrees (23.44)
Q2=23+(26+((21.448-G2*(46.815+G2*(0.00059-G2*0.001813))))/60)/60
#R=Oblique Correction factor (23.44)
R2=Q2+0.00256*math.cos((125.04-1934.136*G2)/r) # new
#S=Right Ascension degrees (120.74)
S2=math.atan2(math.cos(R2/r)*math.sin(P2/r), math.cos(P2/r))*r
T2=math.asin(math.sin(R2/r)*math.sin(P2/r))*r
#U=var y (0.04)
U2=math.tan(R2/2/r)*math.tan(R2/2/r)
#V=EoT minutes (13.55)
V2=4*(U2*math.sin(2*I2/r)-2*K2*math.sin(J2/r)+4*K2*U2*math.sin(J2/r)*math.cos(2*I2/r)-0.5*U2*U2*math.sin(4*I2/r)-1.25*K2*math.sin(2*J2/r))*r
V2=4*(U2*math.sin(2*(I2/r))-2*K2*math.sin((J2/r))+4*K2*U2*math.sin((J2/r))*math.cos(2*(I2/r))-0.5*U2*U2*math.sin(4*(I2/r))-1.25*K2*K2*math.sin(2*(J2/r)))*r
#W=HA Sunrise degrees (75.24)
W2=(math.acos(math.cos(90.833/r)/(math.cos(B3/r)*math.cos(T2/r))-math.tan(B3/r)*math.tan(T2/r)))*r#X=Sun noon LST (0.51)
X2=(720-4*B4-V2+B5*60)/1440
#Y=Sunrise time (0.3)
Y2=X2-W2*4/1440
#Z=Sunset time (0.72)
Z2=X2+W2*4/1440
#AA2=sunlight duration (601.9)
AA2=8*W2
#AB2=TST (404.75)
AB2a=(E2*1440+V2+4*B4-60*B5)
AB2b=1440
AB2=AB2a-int(AB2a/AB2b)*AB2b
if AB2/4 < 0 : AC2=AB2/4+180
else: AC2=AB2/4-180
#AD2=Zenith (93.46)
AD2=math.acos(math.sin(B3/r)*math.sin(T2/r)+math.cos(B3/r)*math.cos(T2/r)*math.cos(AC2/r))*r
#AE2= Sun elevation in degrees (-3.46)
AE2=90-AD2
#AF2=Atmospheric refraction in degrees (0.1)
if AE2>85 : AF2=0
elif AE2>5: AF2=58.1/math.tan((AE2))-0.07/math.pow(math.tan((AE2)),3)+0.000086/math.pow(math.tan((AE2)),5)
elif AE2>-0.575: AF2=1735+AE2*(-518.2+AE2*(103.4+AE2*(-12.79+AE2*0.711)))
else : AF2=-20.772/math.tan(AE2)/3600
 #AG2=Sun Elevation corrected for refraction (-3.37)
AG2=AE2+AF2
#AH2=Sun Azimuth degrees clockwise from North (112.93)
AH2a1=math.acos(((math.sin((B3/r)) * math.cos((AD2/r))) - math.sin((T2/r))) / ((math.cos((B3/r)) * math.sin(AD2/r))))*r + 180
AH2b1=360
AH2a2=540-math.acos(((math.sin(B3/r)*math.cos(AD2/r))-math.sin(T2/r))/((math.cos(B3/r)*math.sin(AD2/r))))*r
AH2b2=AH2b1
if AC2 > 0 : AH2=AH2a1-int(AH2a1/AH2b1)*AH2b1
else : AH2=AH2a2-int(AH2a2/AH2b2)*AH2b2
# error compensation
#AH2=AH2+7
'''
print(Y2+",")
print(Z2,",")
print(AE2,",")
print(AH2,",")
'''
#print (str(X2)+","+str(Y2)+","+str(AE2)+","+str(AH2));
ae2=str(AE2)
ah2=str(AH2)
#=============================================================================
# Si, Sm calculations
r=360/2/3.14159
a=0.14
h=0.535
ah=a*h
alpha=AE2 #40.5603 # Sel
beta=45
Z=90-alpha
psi=180
theta=AH2 #175.9508 # Saz

#AM=1/math.cos(Z/r) #  theta is 90-Sel or 90-alpha - use Z
#for flat earth, no earth curvature
AM=1/(math.cos(Z/r)+0.50572*pow(96.07995-Z/r,-1.6364))

# flat earth, no terrestrial elevation
#Si=1.353*math.pow(0.7,math.pow(AM,0.678))
Si=1.353*((1-ah)*pow(0.7,pow(AM,0.678))+ah)
Sm=Si*(math.cos(alpha/r)*math.sin(beta/r)*math.cos((psi-theta)/r)+math.sin(alpha/r)*math.cos(beta/r))
####################################################################################

print(str(ae2)+','+str(ah2)+','+str(Si)+','+str(Sm));

END CODE