Skip to content
Snippets Groups Projects
Commit 8b25b7e9 authored by Bruce Flynn's avatar Bruce Flynn
Browse files

Fix some issues with solar position

Bump version.
parent b0f8042b
No related branches found
No related tags found
No related merge requests found
...@@ -27,14 +27,14 @@ Simple usage is as follows:: ...@@ -27,14 +27,14 @@ Simple usage is as follows::
solarDec: 4.81517855635 solarDec: 4.81517855635
eqTime: -3.71773858624 eqTime: -3.71773858624
solarConstant: 1369.49188965> solarConstant: 1369.49188965>
The 'SolarPosition' object is a list like structure that allows you to The 'SolarPosition' object is a list like structure that allows you to
access the data in a list like manner:: access the data in a list like manner::
>>> for val in s: >>> for val in s:
print val print val
...: ...:
45 45
-90 -90
0.999456181732 0.999456181732
...@@ -45,18 +45,18 @@ access the data in a list like manner:: ...@@ -45,18 +45,18 @@ access the data in a list like manner::
4.81517855635 4.81517855635
-3.71773858624 -3.71773858624
1369.49188965 1369.49188965
>>> s[1] >>> s[1]
45 45
>>> s[sp.SP_COS_ZEN] >>> s[sp.SP_COS_ZEN]
0.73358139287238222 0.73358139287238222
Additionally, there are property attributes for each value:: Additionally, there are property attributes for each value::
>>> s.solarDec >>> s.solarDec
4.8151785563547014 4.8151785563547014
>>> s.earthSunDistance >>> s.earthSunDistance
0.99945618173240092 0.99945618173240092
...@@ -71,8 +71,7 @@ __docformat__ = 'restructuredtext en' ...@@ -71,8 +71,7 @@ __docformat__ = 'restructuredtext en'
__version__ = '$Revision: 1.6 $' __version__ = '$Revision: 1.6 $'
# $Source: /cvsroot/TOOLS/dev/metobs/python/data/metobs/data/solar_position.py,v $ # $Source: /cvsroot/TOOLS/dev/metobs/python/data/metobs/data/solar_position.py,v $
from metobs import mytime from metobs.data import NaN, util
from metobs.data import NaN, approx_eq
import math import math
_days = ["Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday"] _days = ["Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday"]
...@@ -97,9 +96,9 @@ def calcDayOfWeek(juld): ...@@ -97,9 +96,9 @@ def calcDayOfWeek(juld):
return _days[int(A)] return _days[int(A)]
def calcJD(year, month, day): def calcJD(year, month, day):
"""Get the interval of time in days and fractions of a day, since January 1, """Get the interval of time in days and fractions of a day, since January 1,
4713 BC Greenwich noon, Julian proleptic calendar. 4713 BC Greenwich noon, Julian proleptic calendar.
>>> calcJD(1970,1,1) >>> calcJD(1970,1,1)
2440587.5 2440587.5
""" """
...@@ -109,7 +108,7 @@ def calcJD(year, month, day): ...@@ -109,7 +108,7 @@ def calcJD(year, month, day):
A = math.floor(year/100) A = math.floor(year/100)
B = 2 - A + math.floor(A/4) B = 2 - A + math.floor(A/4)
return (math.floor(365.25 * (year + 4716)) return (math.floor(365.25 * (year + 4716))
+ math.floor(30.6001 * (month + 1)) + math.floor(30.6001 * (month + 1))
+ day + B - 1524.5) + day + B - 1524.5)
def calcDateFromJD(jd): def calcDateFromJD(jd):
...@@ -125,16 +124,16 @@ def calcDateFromJD(jd): ...@@ -125,16 +124,16 @@ def calcDateFromJD(jd):
else: else:
alpha = math.floor((z - 1867216.25)/36524.25) alpha = math.floor((z - 1867216.25)/36524.25)
A = z + 1 + alpha - math.floor(alpha / 4) A = z + 1 + alpha - math.floor(alpha / 4)
B = A + 1524 B = A + 1524
C = math.floor((B - 122.1) / 365.25) C = math.floor((B - 122.1) / 365.25)
D = math.floor(365.25 * C) D = math.floor(365.25 * C)
E = math.floor((B - D) / 30.6001) E = math.floor((B - D) / 30.6001)
day = int(B - D - math.floor(30.6001 * E) + f) day = int(B - D - math.floor(30.6001 * E) + f)
month = int((E < 14) and E - 1 or E - 13) month = int((E < 14) and E - 1 or E - 13)
year = int((month > 2) and C - 4716 or C - 4715) year = int((month > 2) and C - 4716 or C - 4715)
return "%d-%s-%d" % (day, _months[month-1], year) return "%d-%s-%d" % (day, _months[month-1], year)
def calcDayFromJD(jd): def calcDayFromJD(jd):
...@@ -193,9 +192,9 @@ def calcSunEqOfCenter(t): ...@@ -193,9 +192,9 @@ def calcSunEqOfCenter(t):
sinm = math.sin(mrad) sinm = math.sin(mrad)
sin2m = math.sin(mrad * 2) sin2m = math.sin(mrad * 2)
sin3m = math.sin(mrad * 3) sin3m = math.sin(mrad * 3)
return (sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) return (sinm * (1.914602 - t * (0.004817 + 0.000014 * t))
+ sin2m * (0.019993 - 0.000101 * t) + sin2m * (0.019993 - 0.000101 * t)
+ sin3m * 0.000289) + sin3m * 0.000289)
def calcSunTrueLong(t): def calcSunTrueLong(t):
...@@ -286,7 +285,7 @@ def calcEquationOfTime(t): ...@@ -286,7 +285,7 @@ def calcEquationOfTime(t):
Etime = (y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0 Etime = (y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0
- 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m) - 0.5 * y * y * sin4l0 - 1.25 * e * e * sin2m)
return math.degrees(Etime) * 4.0 return math.degrees(Etime) * 4.0
SP_DATE = 0 SP_DATE = 0
SP_LAT = 1 SP_LAT = 1
...@@ -305,7 +304,7 @@ class SolarPosition(object): ...@@ -305,7 +304,7 @@ class SolarPosition(object):
VAL_CNT = 12 VAL_CNT = 12
def __init__(self, lat, lon, earthRadVec, azimuth, solarZen, coszen, def __init__(self, lat, lon, earthRadVec, azimuth, solarZen, coszen,
elevation, solarDec, eqtime, date): elevation, solarDec, eqtime, date):
self._lat = lat self._lat = lat
self._lon = lon self._lon = lon
...@@ -317,14 +316,14 @@ class SolarPosition(object): ...@@ -317,14 +316,14 @@ class SolarPosition(object):
self._solarDec = solarDec self._solarDec = solarDec
self._eqTime = eqtime self._eqTime = eqtime
self._date = date self._date = date
self._cur_idx = 0 self._cur_idx = 0
@property @property
def latitude(self): def latitude(self):
return self._lat return self._lat
@property @property
def longitude(self): def longitude(self):
return self._lon return self._lon
@property @property
def earthSunDistance(self): def earthSunDistance(self):
...@@ -349,7 +348,7 @@ class SolarPosition(object): ...@@ -349,7 +348,7 @@ class SolarPosition(object):
return self._eqTime return self._eqTime
@property @property
def timenow(self): def timenow(self):
return mytime.datetime_to_epoch(self._date) / 3600 return util.to_unix_timestamp(self._date) / 3600
@property @property
def date(self): def date(self):
return self._date return self._date
...@@ -388,7 +387,7 @@ class SolarPosition(object): ...@@ -388,7 +387,7 @@ class SolarPosition(object):
val = self[self._cur_idx] val = self[self._cur_idx]
self._cur_idx += 1 self._cur_idx += 1
return val return val
def __str__(self): def __str__(self):
return """<SolarPosition return """<SolarPosition
date: %s date: %s
...@@ -411,10 +410,10 @@ def solar_position(date, latitude, longitude): ...@@ -411,10 +410,10 @@ def solar_position(date, latitude, longitude):
it assumed to be in localtime. it assumed to be in localtime.
- `latitude`: Station latitude - `latitude`: Station latitude
- `longitude`: Station longitude - `longitude`: Station longitude
Test values were results of running the SolarPosition Java code in cvs at Test values were results of running the SolarPosition Java code in cvs at
https://cvs.ssec.wisc.edu/cgi-bin/cvsweb.cgi/java/tower/SolarPosition.java https://cvs.ssec.wisc.edu/cgi-bin/cvsweb.cgi/java/tower/SolarPosition.java
>>> from datetime import datetime >>> from datetime import datetime
>>> date = datetime(1970,1,1,12,0,0) >>> date = datetime(1970,1,1,12,0,0)
>>> s = solar_position(date, 45, -90) >>> s = solar_position(date, 45, -90)
...@@ -424,19 +423,19 @@ def solar_position(date, latitude, longitude): ...@@ -424,19 +423,19 @@ def solar_position(date, latitude, longitude):
>>> assert approx_eq(s.elevation, 22.038213553352705, _ACCEPTABLE_DELTA) >>> assert approx_eq(s.elevation, 22.038213553352705, _ACCEPTABLE_DELTA)
>>> assert approx_eq(s.cosZen, 0.3752248974802099, _ACCEPTABLE_DELTA) >>> assert approx_eq(s.cosZen, 0.3752248974802099, _ACCEPTABLE_DELTA)
""" """
if not date.tzinfo:
date = mytime.set_tz(date, tz=mytime.LocalTimezone)
year, month, day, hh, mm, ss = date.timetuple()[:6] year, month, day, hh, mm, ss = date.timetuple()[:6]
tzoffset = date.tzinfo.utcoffset(date) GMToffset = 0
if tzoffset.days < 0: if date.tzinfo:
GMToffset = -(24 - tzoffset.seconds/3600) tzoffset = date.tzinfo.utcoffset(date)
else: if tzoffset.days < 0:
GMToffset = tzoffset.seconds/3600 GMToffset = -(24 - tzoffset.seconds/3600)
else:
GMToffset = tzoffset.seconds/3600
timenow = hh + mm/60 + ss/3600 - GMToffset timenow = hh + mm/60 + ss/3600 - GMToffset
JD = calcJD(year, month, day) JD = calcJD(year, month, day)
T = calcTimeJulianCent(JD + timenow/24.0) T = calcTimeJulianCent(JD + timenow/24.0)
solarDec = calcSunDeclination(T) solarDec = calcSunDeclination(T)
eqTime = calcEquationOfTime(T) eqTime = calcEquationOfTime(T)
earthRadVec = calcSunRadVector(T) earthRadVec = calcSunRadVector(T)
...@@ -453,9 +452,9 @@ def solar_position(date, latitude, longitude): ...@@ -453,9 +452,9 @@ def solar_position(date, latitude, longitude):
haRad = math.radians(hourAngle) haRad = math.radians(hourAngle)
csz = (math.sin(math.radians(latitude)) * csz = (math.sin(math.radians(latitude)) *
math.sin(math.radians(solarDec)) + math.sin(math.radians(solarDec)) +
math.cos(math.radians(latitude)) * math.cos(math.radians(latitude)) *
math.cos(math.radians(solarDec)) * math.cos(haRad)) math.cos(math.radians(solarDec)) * math.cos(haRad))
if csz > 1.0: if csz > 1.0:
...@@ -468,8 +467,8 @@ def solar_position(date, latitude, longitude): ...@@ -468,8 +467,8 @@ def solar_position(date, latitude, longitude):
math.sin(math.radians(zenith)) ) math.sin(math.radians(zenith)) )
if abs(azDenom) > 0.001: if abs(azDenom) > 0.001:
azRad = ((( math.sin(math.radians(latitude)) * azRad = ((( math.sin(math.radians(latitude)) *
math.cos(math.radians(zenith)) ) - math.cos(math.radians(zenith)) ) -
math.sin(math.radians(solarDec))) / azDenom) math.sin(math.radians(solarDec))) / azDenom)
if abs(azRad) > 1.0: if abs(azRad) > 1.0:
...@@ -501,12 +500,12 @@ def solar_position(date, latitude, longitude): ...@@ -501,12 +500,12 @@ def solar_position(date, latitude, longitude):
else: else:
te = math.tan (math.radians(exoatmElevation)) te = math.tan (math.radians(exoatmElevation))
if (exoatmElevation > 5.0): if (exoatmElevation > 5.0):
refractionCorrection = (58.1 / te - 0.07 / (te*te*te) refractionCorrection = (58.1 / te - 0.07 / (te*te*te)
+ 0.000086 / (te*te*te*te*te)) + 0.000086 / (te*te*te*te*te))
elif exoatmElevation > -0.575: elif exoatmElevation > -0.575:
refractionCorrection = (1735.0 + exoatmElevation * refractionCorrection = (1735.0 + exoatmElevation *
(-518.2 + exoatmElevation * (103.4 + (-518.2 + exoatmElevation * (103.4 +
exoatmElevation * (-12.79 + exoatmElevation * (-12.79 +
exoatmElevation * 0.711) ) )) exoatmElevation * 0.711) ) ))
else: else:
refractionCorrection = -20.774 / te refractionCorrection = -20.774 / te
...@@ -526,8 +525,8 @@ def solar_position(date, latitude, longitude): ...@@ -526,8 +525,8 @@ def solar_position(date, latitude, longitude):
azumuth = NaN azumuth = NaN
elevation = NaN elevation = NaN
coszen = NaN coszen = NaN
return SolarPosition(latitude, longitude, earthRadVec, azimuth, solarZen, coszen, return SolarPosition(latitude, longitude, earthRadVec, azimuth, solarZen, coszen,
elevation, solarDec, eqTime, date) elevation, solarDec, eqTime, date)
if __name__ == '__main__': if __name__ == '__main__':
......
...@@ -3,7 +3,7 @@ from setuptools import setup, find_packages ...@@ -3,7 +3,7 @@ from setuptools import setup, find_packages
setup( setup(
name='metobs.data', name='metobs.data',
version='0.5dev', version='0.5',
packages=find_packages(exclude=["*.tests", "*.tests.*", "tests.*", "tests"]), packages=find_packages(exclude=["*.tests", "*.tests.*", "tests.*", "tests"]),
namespace_packages=['metobs'], namespace_packages=['metobs'],
include_package_data=True, include_package_data=True,
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment