From 8b25b7e98d0997aed114c2c26e72c917aa60521a Mon Sep 17 00:00:00 2001
From: Bruce Flynn <brucef@ssec.wisc.edu>
Date: Thu, 11 Sep 2014 16:21:20 -0500
Subject: [PATCH] Fix some issues with solar position Bump version.
---
metobs/data/solar_position.py | 85 +++++++++++++++++------------------
setup.py | 2 +-
2 files changed, 43 insertions(+), 44 deletions(-)
diff --git a/metobs/data/solar_position.py b/metobs/data/solar_position.py
index ed01998..da81676 100644
--- a/metobs/data/solar_position.py
+++ b/metobs/data/solar_position.py
@@ -27,14 +27,14 @@ Simple usage is as follows::
solarDec: 4.81517855635
eqTime: -3.71773858624
solarConstant: 1369.49188965>
-
+
The 'SolarPosition' object is a list like structure that allows you to
access the data in a list like manner::
>>> for val in s:
print val
- ...:
+ ...:
45
-90
0.999456181732
@@ -45,18 +45,18 @@ access the data in a list like manner::
4.81517855635
-3.71773858624
1369.49188965
-
+
>>> s[1]
45
-
+
>>> s[sp.SP_COS_ZEN]
0.73358139287238222
Additionally, there are property attributes for each value::
-
+
>>> s.solarDec
4.8151785563547014
-
+
>>> s.earthSunDistance
0.99945618173240092
@@ -71,8 +71,7 @@ __docformat__ = 'restructuredtext en'
__version__ = '$Revision: 1.6 $'
# $Source: /cvsroot/TOOLS/dev/metobs/python/data/metobs/data/solar_position.py,v $
-from metobs import mytime
-from metobs.data import NaN, approx_eq
+from metobs.data import NaN, util
import math
_days = ["Sunday","Monday","Tuesday","Wednesday","Thursday","Friday","Saturday"]
@@ -97,9 +96,9 @@ def calcDayOfWeek(juld):
return _days[int(A)]
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.
-
+
>>> calcJD(1970,1,1)
2440587.5
"""
@@ -109,7 +108,7 @@ def calcJD(year, month, day):
A = math.floor(year/100)
B = 2 - A + math.floor(A/4)
return (math.floor(365.25 * (year + 4716))
- + math.floor(30.6001 * (month + 1))
+ + math.floor(30.6001 * (month + 1))
+ day + B - 1524.5)
def calcDateFromJD(jd):
@@ -125,16 +124,16 @@ def calcDateFromJD(jd):
else:
alpha = math.floor((z - 1867216.25)/36524.25)
A = z + 1 + alpha - math.floor(alpha / 4)
-
+
B = A + 1524
C = math.floor((B - 122.1) / 365.25)
D = math.floor(365.25 * C)
E = math.floor((B - D) / 30.6001)
-
+
day = int(B - D - math.floor(30.6001 * E) + f)
month = int((E < 14) and E - 1 or E - 13)
year = int((month > 2) and C - 4716 or C - 4715)
-
+
return "%d-%s-%d" % (day, _months[month-1], year)
def calcDayFromJD(jd):
@@ -193,9 +192,9 @@ def calcSunEqOfCenter(t):
sinm = math.sin(mrad)
sin2m = math.sin(mrad * 2)
sin3m = math.sin(mrad * 3)
-
- return (sinm * (1.914602 - t * (0.004817 + 0.000014 * t))
- + sin2m * (0.019993 - 0.000101 * t)
+
+ return (sinm * (1.914602 - t * (0.004817 + 0.000014 * t))
+ + sin2m * (0.019993 - 0.000101 * t)
+ sin3m * 0.000289)
def calcSunTrueLong(t):
@@ -286,7 +285,7 @@ def calcEquationOfTime(t):
Etime = (y * sin2l0 - 2.0 * e * sinm + 4.0 * e * y * sinm * cos2l0
- 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_LAT = 1
@@ -305,7 +304,7 @@ class SolarPosition(object):
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):
self._lat = lat
self._lon = lon
@@ -317,14 +316,14 @@ class SolarPosition(object):
self._solarDec = solarDec
self._eqTime = eqtime
self._date = date
-
+
self._cur_idx = 0
@property
- def latitude(self):
+ def latitude(self):
return self._lat
@property
- def longitude(self):
+ def longitude(self):
return self._lon
@property
def earthSunDistance(self):
@@ -349,7 +348,7 @@ class SolarPosition(object):
return self._eqTime
@property
def timenow(self):
- return mytime.datetime_to_epoch(self._date) / 3600
+ return util.to_unix_timestamp(self._date) / 3600
@property
def date(self):
return self._date
@@ -388,7 +387,7 @@ class SolarPosition(object):
val = self[self._cur_idx]
self._cur_idx += 1
return val
-
+
def __str__(self):
return """<SolarPosition
date: %s
@@ -411,10 +410,10 @@ def solar_position(date, latitude, longitude):
it assumed to be in localtime.
- `latitude`: Station latitude
- `longitude`: Station longitude
-
+
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
-
+
>>> from datetime import datetime
>>> date = datetime(1970,1,1,12,0,0)
>>> s = solar_position(date, 45, -90)
@@ -424,19 +423,19 @@ def solar_position(date, latitude, longitude):
>>> assert approx_eq(s.elevation, 22.038213553352705, _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]
- tzoffset = date.tzinfo.utcoffset(date)
- if tzoffset.days < 0:
- GMToffset = -(24 - tzoffset.seconds/3600)
- else:
- GMToffset = tzoffset.seconds/3600
+ GMToffset = 0
+ if date.tzinfo:
+ tzoffset = date.tzinfo.utcoffset(date)
+ if tzoffset.days < 0:
+ GMToffset = -(24 - tzoffset.seconds/3600)
+ else:
+ GMToffset = tzoffset.seconds/3600
timenow = hh + mm/60 + ss/3600 - GMToffset
JD = calcJD(year, month, day)
- T = calcTimeJulianCent(JD + timenow/24.0)
+ T = calcTimeJulianCent(JD + timenow/24.0)
solarDec = calcSunDeclination(T)
eqTime = calcEquationOfTime(T)
earthRadVec = calcSunRadVector(T)
@@ -453,9 +452,9 @@ def solar_position(date, latitude, longitude):
haRad = math.radians(hourAngle)
- csz = (math.sin(math.radians(latitude)) *
- math.sin(math.radians(solarDec)) +
- math.cos(math.radians(latitude)) *
+ csz = (math.sin(math.radians(latitude)) *
+ math.sin(math.radians(solarDec)) +
+ math.cos(math.radians(latitude)) *
math.cos(math.radians(solarDec)) * math.cos(haRad))
if csz > 1.0:
@@ -468,8 +467,8 @@ def solar_position(date, latitude, longitude):
math.sin(math.radians(zenith)) )
if abs(azDenom) > 0.001:
- azRad = ((( math.sin(math.radians(latitude)) *
- math.cos(math.radians(zenith)) ) -
+ azRad = ((( math.sin(math.radians(latitude)) *
+ math.cos(math.radians(zenith)) ) -
math.sin(math.radians(solarDec))) / azDenom)
if abs(azRad) > 1.0:
@@ -501,12 +500,12 @@ def solar_position(date, latitude, longitude):
else:
te = math.tan (math.radians(exoatmElevation))
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))
elif exoatmElevation > -0.575:
refractionCorrection = (1735.0 + exoatmElevation *
(-518.2 + exoatmElevation * (103.4 +
- exoatmElevation * (-12.79 +
+ exoatmElevation * (-12.79 +
exoatmElevation * 0.711) ) ))
else:
refractionCorrection = -20.774 / te
@@ -526,8 +525,8 @@ def solar_position(date, latitude, longitude):
azumuth = NaN
elevation = NaN
coszen = NaN
-
- return SolarPosition(latitude, longitude, earthRadVec, azimuth, solarZen, coszen,
+
+ return SolarPosition(latitude, longitude, earthRadVec, azimuth, solarZen, coszen,
elevation, solarDec, eqTime, date)
if __name__ == '__main__':
diff --git a/setup.py b/setup.py
index bec967f..76c2aa2 100644
--- a/setup.py
+++ b/setup.py
@@ -3,7 +3,7 @@ from setuptools import setup, find_packages
setup(
name='metobs.data',
- version='0.5dev',
+ version='0.5',
packages=find_packages(exclude=["*.tests", "*.tests.*", "tests.*", "tests"]),
namespace_packages=['metobs'],
include_package_data=True,
--
GitLab