Applied the 80/20 rule to a NOAA spreadsheet I found (that seems to have gone missing since the other day). It’s within a couple minutes for sunrise, noon and sunset and that’s close enough for what I’m doing.
Lots of room for improvement and not something I’d sign my name to at GitHub at this level. Thought I’d share here in case anybody else needed something basic.
- Cut and paste the code into a shared gateway script library called something like noaa_solar
- Edit for your lat, lon and tz and preferred root tag folder
- Call it from a gateway timer script with “shared.noaa_solar.begin()” say about once a minute
[code]’’’
2016-03-12
http://www.esrl.noaa.gov/gmd/grad/solcalc/sunrise.html
https://www.google.com/search?q=NOAA_Solar_Calculations_day&ie=utf-8&oe=utf-8#q=NOAA+Solar+Calculations+day
‘’’
import datetime, math
def begin():
def mkdir(par, dir):
''' make a folder if necessary '''
import system
if (system.tag.exists(par + "/" + dir) == 0):
system.tag.addTag(parentPath=par, name=dir, tagType="Folder")
return
def mktag(par, tag, typ, val):
''' make a tag if necessary then assign a value '''
import system
if (system.tag.exists(par + "/" + tag) == 0):
system.tag.addTag(parentPath=par, name=tag, tagType="MEMORY", dataType=typ)
system.tag.write(par + "/" + tag, val)
return
def mkloc(par, lat, lon, tz):
mktag(par, "Latitude", "Float8", lat)
mktag(par, "Longitude", "Float8", lon)
mktag(par, "TimeZone", "Float8", tz)
dt = datetime.datetime.now()
#dt = datetime.datetime(2016, 3, 11, 0, 6, 0) #for testing
mktag(par, "CurrentTime", "DateTime", dt)
#print "date %s " % dt
# E
TimePastLocalMidnight = (dt - dt.replace(hour=0, minute=0, second=0, microsecond=0)).seconds / 60.0 / 1440
Midnight = dt.replace(hour=0, minute=0, second=0, microsecond=0)
mktag(par, "TimePastLocalMidnight", "Float8", TimePastLocalMidnight)
#print "%s TimePastLocalMidnight" % TimePastLocalMidnight
# F
#ts = (dt - datetime.datetime(1970, 1, 1)).total_seconds()
#revised for compatibility with Ignition
dt = datetime.datetime(2016, 3, 11, 0, 6, 0) - datetime.datetime(1970, 1, 1)
ts = dt.days * 86400.0 + dt.seconds
JulianDay = ts / 86400.0 + 2440587.5 - tz / 24
mktag(par, "JulianDay", "Float8", JulianDay)
#print "%s JulianDay \n2457458.71" % JulianDay
# G
JulianCentury = (JulianDay - 2451545) / 36525
mktag(par, "JulianCentury", "Float8", JulianCentury)
#print "%s JulianCentury \n0.16190862" % JulianCentury
# H
GeomMeanLongSunDeg = 280.46646 + JulianCentury * (36000.76983 + JulianCentury * 0.0003032) % 360
mktag(par, "GeomMeanLongSunDeg", "Float8", GeomMeanLongSunDeg)
#print "%s GeomMeanLongSunDeg \n349.3015823" % GeomMeanLongSunDeg
# I
GeomMeanAnomSunDeg = 357.52911 + JulianCentury * (35999.05029 - 0.0001537 * JulianCentury)
mktag(par, "GeomMeanAnomSunDeg", "Float8", GeomMeanAnomSunDeg)
#print "%s GeomMeanAnomSunDeg \n6186.085812" % GeomMeanAnomSunDeg
# J
EccentEarthOrbit = 0.016708634 - JulianCentury * (0.000042037 + 0.0000001267 * JulianCentury)
mktag(par, "EccentEarthOrbit", "Float8", EccentEarthOrbit)
#print "%s EccentEarthOrbit \n0.016701825" % EccentEarthOrbit
# K
SunEqofCtr = math.sin(math.radians(GeomMeanAnomSunDeg)) * (1.914602 - JulianCentury * (0.004817 + 0.000014 * JulianCentury)) \
+ math.sin(math.radians(2 * GeomMeanAnomSunDeg)) * (0.019993 - 0.000101 * JulianCentury) \
+ math.sin(math.radians(3 * GeomMeanAnomSunDeg)) * 0.000289
mktag(par, "SunEqofCtr", "Float8", SunEqofCtr)
#print "%s SunEqofCtr \n1.764241916" % SunEqofCtr
# L
SunTrueLongDeg = GeomMeanLongSunDeg + SunEqofCtr
mktag(par, "SunTrueLongDeg", "Float8", SunTrueLongDeg)
#print "%s SunTrueLongDeg \n351.0658243" % SunTrueLongDeg
# M
SunTrueAnomDeg = GeomMeanAnomSunDeg + SunEqofCtr
mktag(par, "SunTrueAnomDeg", "Float8", SunTrueAnomDeg)
#print "%s SunTrueAnomDeg \n6187.850054" % SunTrueAnomDeg
# N
SunRadVectorAUs = (1.000001018 * (1 - EccentEarthOrbit * EccentEarthOrbit)) / (1 + EccentEarthOrbit * math.cos(math.radians(SunTrueAnomDeg)))
mktag(par, "SunRadVectorAUs", "Float8", SunRadVectorAUs)
#print "%s SunRadVectorAUs \n0.993466093" % SunRadVectorAUs
# O
SunAppLongDeg = SunTrueLongDeg - 0.00569 - 0.00478 * math.sin(math.radians(125.04 - 1934.136 * JulianCentury))
mktag(par, "SunAppLongDeg", "Float8", SunAppLongDeg)
#print "%s SunAppLongDeg \n351.0594597" % SunAppLongDeg
# P
MeanObliqEclipticDeg = 23 + (26 + ((21.448 - JulianCentury * (46.815 + JulianCentury * (0.00059 - JulianCentury * 0.001813))))/60)/60
mktag(par, "MeanObliqEclipticDeg", "Float8", MeanObliqEclipticDeg)
#print "%s MeanObliqEclipticDeg \n23.43718562" % MeanObliqEclipticDeg
# Q
ObliqCorrDeg = MeanObliqEclipticDeg + 0.00256 * math.cos(math.radians(125.04 - 1934.136 * JulianCentury))
mktag(par, "ObliqCorrDeg", "Float8", ObliqCorrDeg)
#print "%s ObliqCorrDeg \n23.43465125" % ObliqCorrDeg
# R
SunRtAscenDeg = math.degrees(math.atan2(math.cos(math.radians(SunAppLongDeg)),math.cos(math.radians(ObliqCorrDeg))*math.sin(math.radians(SunAppLongDeg))))
mktag(par, "SunRtAscenDeg", "Float8", SunRtAscenDeg)
#print "%s SunRtAscenDeg \n-8.213576963" % SunRtAscenDeg
# S
SunDeclinDeg = math.degrees(math.asin(math.sin(math.radians(ObliqCorrDeg))*math.sin(math.radians(SunAppLongDeg))))
mktag(par, "SunDeclinDeg", "Float8", SunDeclinDeg)
#print "%s SunDeclinDeg \n-3.543524949" % SunDeclinDeg
# T
vary = math.tan(math.radians(ObliqCorrDeg/2)) * math.tan(math.radians(ObliqCorrDeg/2))
mktag(par, "vary", "Float8", vary)
#print "%s vary \n0.043017009" % vary
# U
EqofTimeMin = 4 * math.degrees(vary * math.sin(2 * math.radians(GeomMeanLongSunDeg)) \
- 2 * EccentEarthOrbit * math.sin(math.radians(GeomMeanAnomSunDeg)) \
+ 4 * EccentEarthOrbit * vary * math.sin(math.radians(GeomMeanAnomSunDeg)) * math.cos(2 * math.radians(GeomMeanLongSunDeg)) \
- 0.5 * vary * vary * math.sin(4 * math.radians(GeomMeanLongSunDeg)) \
- 1.25 * EccentEarthOrbit * EccentEarthOrbit * math.sin(2 * math.radians(GeomMeanAnomSunDeg)))
mktag(par, "EqofTimeMin", "Float8", EqofTimeMin)
#print "%s EqofTimeMin \n-9.949653921" % EqofTimeMin
# V
HASunriseDeg = math.degrees(math.acos(math.cos(math.radians(90.833)) / (math.cos(math.radians(lat)) * math.cos(math.radians(SunDeclinDeg))) \
- math.tan(math.radians(lat)) * math.tan(math.radians(SunDeclinDeg))))
mktag(par, "HASunriseDeg", "Float8", HASunriseDeg)
#print "%s HASunriseDeg \n88.11050567" % HASunriseDeg
# W
SolarNoonLST = (720 - 4 * lon - EqofTimeMin + tz * 60) / 1440
#print Midnight + datetime.timedelta(seconds = 86400 * SolarNoonLST)
mktag(par, "SolarNoonLST", "DateTime", Midnight + datetime.timedelta(seconds = 86400 * SolarNoonLST))
#print "%s SolarNoonLST %s \n12:13:13" % (str(datetime.timedelta(seconds=SolarNoonLST * 86400)), SolarNoonLST)
# X
SunriseTimeLST = SolarNoonLST - HASunriseDeg * 4 / 1440
mktag(par, "SunriseTimeLST", "DateTime", Midnight + datetime.timedelta(seconds = 86400 * SunriseTimeLST))
#print "%s SunriseTimeLST %s \n6:20:47" % (str(datetime.timedelta(seconds=SunriseTimeLST * 86400)), SunriseTimeLST)
# Y
SunsetTimeLST = SolarNoonLST + HASunriseDeg * 4 / 1440
mktag(par, "SunsetTimeLST", "DateTime", Midnight + datetime.timedelta(seconds = 86400 * SunsetTimeLST))
#print "%s SunsetTimeLST %s \n18:05:40" % (str(datetime.timedelta(seconds=SunsetTimeLST * 86400)), SunsetTimeLST)
# Z
SunlightDurationMin = 8 * HASunriseDeg
mktag(par, "SunlightDurationMin", "Float8", SunlightDurationMin)
#print "%s SunlightDurationMin \n704.8840453" % SunlightDurationMin
# AA
TrueSolarTimeMin = (TimePastLocalMidnight * 1440 + EqofTimeMin + 4 * lon - 60 * tz) % 1440
mktag(par, "TrueSolarTimeMin", "Float8", TrueSolarTimeMin)
#print "%s TrueSolarTimeMin \n1432.778346" % TrueSolarTimeMin
# AB
#HourAngleDeg = IF(AA2/4<0,AA2/4+180,AA2/4-180)
HourAngleDeg = TrueSolarTimeMin/4 + 180 if TrueSolarTimeMin/4 < 0 else TrueSolarTimeMin/4 - 180
mktag(par, "HourAngleDeg", "Float8", HourAngleDeg)
#print "%s HourAngleDeg \n178.1945865" % HourAngleDeg
# AC
SolarZenithAngleDeg = math.degrees(math.acos(math.sin(math.radians(lat)) * math.sin(math.radians(SunDeclinDeg)) + \
math.cos(math.radians(lat)) * math.cos(math.radians(SunDeclinDeg)) * math.cos(math.radians(HourAngleDeg))))
mktag(par, "SolarZenithAngleDeg", "Float8", SolarZenithAngleDeg)
#print "%s SolarZenithAngleDeg \n143.4909661" % SolarZenithAngleDeg
# AD
SolarElevationAngleDeg = 90 - SolarZenithAngleDeg
mktag(par, "SolarElevationAngleDeg", "Float8", SolarElevationAngleDeg)
#print "%s SolarElevationAngleDeg \n-53.49096606" % SolarElevationAngleDeg
# AE
if SolarElevationAngleDeg > 85:
ApproxAtmosphericRefractionDeg = 0,
elif SolarElevationAngleDeg > 5:
ApproxAtmosphericRefractionDeg = 58.1 / math.tan(math.radians(SolarElevationAngleDeg)) - 0.07 / math.pow(math.tan(math.radians(SolarElevationAngleDeg)),3) + 0.000086 / math.pow(math.tan(math.radians(SolarElevationAngleDeg)),5)
elif SolarElevationAngleDeg > -0.575:
ApproxAtmosphericRefractionDeg = 1735 + SolarElevationAngleDeg * (-518.2 + SolarElevationAngleDeg * (103.4 + SolarElevationAngleDeg * (-12.79 + SolarElevationAngleDeg * 0.711)))
else:
ApproxAtmosphericRefractionDeg = -20.772 / math.tan(math.radians(SolarElevationAngleDeg))
ApproxAtmosphericRefractionDeg /= 3600
mktag(par, "ApproxAtmosphericRefractionDeg", "Float8", ApproxAtmosphericRefractionDeg)
#print "%s ApproxAtmosphericRefractionDeg \n0.004270983" % ApproxAtmosphericRefractionDeg
# AF
SolarElevationCorrectedForATMRefractionDeg = SolarElevationAngleDeg + ApproxAtmosphericRefractionDeg
mktag(par, "SolarElevationCorrectedForATMRefractionDeg", "Float8", SolarElevationCorrectedForATMRefractionDeg)
#print "%s SolarElevationCorrectedForATMRefractionDeg \n-53.48669508" % SolarElevationCorrectedForATMRefractionDeg
# AG
if HourAngleDeg > 0:
SolarAzimuthAngleDegCWfromN = math.degrees(math.acos(((math.sin(math.radians(lat)) * math.cos(math.radians(SolarZenithAngleDeg))) \
- math.sin(math.radians(SunDeclinDeg))) / (math.cos(math.radians(lat)) * math.sin(math.radians(SolarZenithAngleDeg)))))+180
else:
SolarAzimuthAngleDegCWfromN = 540 - math.degrees(math.acos(((math.sin(math.radians(lat)) * math.cos(math.radians(SolarZenithAngleDeg))) \
- math.sin(math.radians(SunDeclinDeg))) / (math.cos(math.radians(lat)) * math.sin(math.radians(SolarZenithAngleDeg)))))
SolarAzimuthAngleDegCWfromN = SolarAzimuthAngleDegCWfromN % 360
mktag(par, "SolarAzimuthAngleDegCWfromN", "Float8", SolarAzimuthAngleDegCWfromN)
#print "%s SolarAzimuthAngleDegCWfromN \n356.9703254" % SolarAzimuthAngleDegCWfromN
#print "ok"
return
root = "astro" # Ignition tag root folder
# check root exists
mkdir("", root)
# create or update locations
mkloc(root + "/phl", 39.950, -75.167, -5)
return
[/code]