[Python-modules-commits] [python-pysolar] 03/10: Fixed error in nutation calculation. Replaced pow() with ** operator. Refined unit tests a little. Found calculaton of GetMeanSiderealTime is only good to 5 significant figures relative to example given by Reda and Andreas.
Wolfgang Borgert
debacle at moszumanska.debian.org
Fri Oct 3 23:36:03 UTC 2014
This is an automated email from the git hooks/post-receive script.
debacle pushed a commit to annotated tag 0.2.1
in repository python-pysolar.
commit b52be848879b1e4ff93daddef218fb96bf552b77
Author: Brandon Stafford <brandon at farb.(none)>
Date: Sun Mar 2 22:58:50 2008 -0500
Fixed error in nutation calculation. Replaced pow() with ** operator. Refined unit tests a little. Found calculaton of GetMeanSiderealTime is only good to 5 significant figures relative to example given by Reda and Andreas.
---
solar.py | 33 +++++++++++++++++----------------
testsolar.py | 41 +++++++++++++++++++++++------------------
2 files changed, 40 insertions(+), 34 deletions(-)
diff --git a/solar.py b/solar.py
index a1a41a0..e2ff6de 100644
--- a/solar.py
+++ b/solar.py
@@ -75,7 +75,7 @@ def GetApparentSunLongitude(geocentric_longitude, nutation, ab_correction):
return geocentric_longitude + nutation['longitude'] + ab_correction
def GetArgumentOfLatitudeOfMoon(jce):
- return 93.27191 + (483202.017538 * jce) - (0.0036825 * pow(jce, 2)) + (pow(jce, 3) / 327270.0)
+ return 93.27191 + 483202.017538 * jce - 0.0036825 * jce ** 2 + jce ** 3 / 327270.0
def GetAzimuth(latitude_deg, longitude_deg, utc_datetime):
# expect -50 degrees for solar.GetAzimuth(42.364908,-71.112828,datetime.datetime(2007, 2, 18, 20, 18, 0, 0))
@@ -153,7 +153,7 @@ def GetHeliocentricLongitude(jme):
l4 = GetCoefficient(jme, constants.L4)
l5 = GetCoefficient(jme, constants.L5)
- l = (l0 + (l1 * jme) + (l2 * pow(jme, 2)) + (l3 * pow(jme, 3)) + (l4 * pow(jme, 4)) + (l5 * pow(jme, 5))) / pow(10, 8)
+ l = (l0 + l1 * jme + l2 * jme ** 2 + l3 * jme ** 3 + l4 * jme ** 4 + l5 * jme ** 5) / 10 ** 8
return math.degrees(l) % 360
def GetHourAngle(utc_datetime, longitude_deg):
@@ -168,12 +168,10 @@ def GetIncidenceAngle(topocentric_zenith_angle, slope, slope_orientation, topoce
return math.degrees(math.acos(math.cos(tza_rad) * math.cos(slope_rad) + math.sin(slope_rad) * math.sin(tza_rad) * math.cos(taa_rad - math.pi - so_rad)))
def GetJulianCentury(julian_day):
- """You get the Julian century or Julian ephemeris century back, depending on whether you supply
- the Julian day or the Julian ephemeris day."""
return (julian_day - 2451545.0) / 36525.0
def GetJulianDay(utc_datetime): # based on NREL/TP-560-34302 by Andreas and Reda
- # does not accept years before 0 because of bounds check on Python's datetime.year field
+ # does not accept years before 0 because of bounds check on Python's datetime.year field
year = utc_datetime.year
month = utc_datetime.month
if(month <= 2): # shift to accomodate leap years?
@@ -187,6 +185,9 @@ def GetJulianDay(utc_datetime): # based on NREL/TP-560-34302 by Andreas and Reda
else:
return julian_day + gregorian_offset # after October 5, 1852
+def GetJulianEphemerisCentury(julian_ephemeris_day):
+ return (julian_ephemeris_day - 2451545.0) / 36525.0
+
def GetJulianEphemerisDay(julian_day, delta_seconds):
"""delta_seconds is value referred to by astronomers as Delta-T, defined as the difference between
Dynamical Time (TD) and Universal Time (UT). In 2007, it's around 65 seconds.
@@ -197,27 +198,27 @@ def GetJulianEphemerisMillenium(julian_ephemeris_century):
return (julian_ephemeris_century / 10.0)
def GetLongitudeOfAscendingNode(jce):
- return 125.04452 - (1934.136261 * jce) + (0.0020708 * pow(jce, 2)) + (pow(jce, 3) / 450000.0)
+ return 125.04452 - (1934.136261 * jce) + (0.0020708 * jce ** 2) + (jce ** 3 / 450000.0)
def GetLocalHourAngle(apparent_sidereal_time, longitude, geocentric_sun_right_ascension):
return (apparent_sidereal_time + longitude - geocentric_sun_right_ascension) % 360
def GetMeanElongationOfMoon(jce):
- return 297.85036 + (445267.111480 * jce) - (0.0019142 * pow(jce, 2)) + (pow(jce, 3) / 189474.0)
+ return 297.85036 + (445267.111480 * jce) - (0.0019142 * jce ** 2) + (jce ** 3 / 189474.0)
def GetMeanAnomalyOfMoon(jce):
- return 134.96298 + (477198.867398 * jce) + (0.0086972 * pow(jce, 2)) + (pow(jce, 3) / 56250.0)
+ return 134.96298 + (477198.867398 * jce) + (0.0086972 * jce ** 2) + (jce ** 3 / 56250.0)
def GetMeanAnomalyOfSun(jce):
- return 357.52772 + (35999.050340 * jce) - (0.0001603 * pow(jce, 2)) - (pow(jce, 3) / 300000.0)
+ return 357.52772 + (35999.050340 * jce) - (0.0001603 * jce ** 2) - (jce ** 3 / 300000.0)
def GetMeanSiderealTime(julian_day):
+ # This function doesn't agree with Andreas and Reda as well as it should. Works to ~5 sig figs in current unit test
jc = GetJulianCentury(julian_day)
- sidereal_time = 280.46061837 + (360.98564736629 * (julian_day - 2451545.0)) + (0.000387933 * pow(jc, 2)) \
- - (pow(jc, 3) / 38710000)
+ sidereal_time = 280.46061837 + (360.98564736629 * (julian_day - 2451545.0)) + (0.000387933 * jc ** 2) - (jc ** 3 / 38710000)
return sidereal_time % 360
-def GetNutationAberrationXY(jce):
+def GetNutationAberrationXY(jce, i):
y = constants.aberration_sin_terms
x = []
# order of 5 x.append lines below is important
@@ -228,17 +229,17 @@ def GetNutationAberrationXY(jce):
x.append(GetLongitudeOfAscendingNode(jce))
sigmaxy = 0.0
for j in range(len(x)):
- sigmaxy += x[j] * y[0][j]
+ sigmaxy += x[j] * y[i][j]
return sigmaxy
def GetNutation(jde):
abcd = constants.nutation_coefficients
- jce = GetJulianCentury(jde)
- sigmaxy = GetNutationAberrationXY(jce)
+ jce = GetJulianEphemerisCentury(jde)
nutation_long = []
nutation_oblique = []
for i in range(len(abcd)):
+ sigmaxy = GetNutationAberrationXY(jce, i)
nutation_long.append((abcd[i][0] + (abcd[i][1] * jce)) * math.sin(math.radians(sigmaxy)))
nutation_oblique.append((abcd[i][2] + (abcd[i][3] * jce)) * math.cos(math.radians(sigmaxy)))
@@ -287,7 +288,7 @@ def GetRadiusVector(jme):
r3 = GetCoefficient(jme, constants.R3)
r4 = GetCoefficient(jme, constants.R4)
- return (r0 + (r1 * jme) + (r2 * pow(jme, 2)) + (r3 * pow(jme, 3)) + (r4 * pow(jme, 4))) / pow(10, 8)
+ return (r0 + r1 * jme + r2 * jme ** 2 + r3 * jme ** 3 + r4 * jme ** 4) / 10 ** 8
def GetRefractionCorrection(pressure_millibars, temperature_celsius, topocentric_elevation_angle):
tea = topocentric_elevation_angle
diff --git a/testsolar.py b/testsolar.py
index 5b2f0ba..c719e20 100644
--- a/testsolar.py
+++ b/testsolar.py
@@ -36,8 +36,9 @@ class testSolar(unittest.TestCase):
self.slope = 30.0 # degrees
self.slope_orientation = -10.0 # degrees east from south
self.jd = solar.GetJulianDay(self.d)
- self.jde = solar.GetJulianEphemerisDay(self.jd, 65.284)
- self.jce = solar.GetJulianCentury(self.jde)
+ self.jc = solar.GetJulianCentury(self.jd)
+ self.jde = solar.GetJulianEphemerisDay(self.jd, 67.0) #65.284)
+ self.jce = solar.GetJulianEphemerisCentury(self.jde)
self.jme = solar.GetJulianEphemerisMillenium(self.jce)
self.geocentric_longitude = solar.GetGeocentricLongitude(self.jme)
self.geocentric_latitude = solar.GetGeocentricLatitude(self.jme)
@@ -49,7 +50,7 @@ class testSolar(unittest.TestCase):
self.apparent_sidereal_time = solar.GetApparentSiderealTime(self.jd, self.jme, self.nutation)
self.geocentric_sun_right_ascension = solar.GetGeocentricSunRightAscension(self.apparent_sun_longitude, self.true_ecliptic_obliquity, self.geocentric_latitude)
self.geocentric_sun_declination = solar.GetGeocentricSunDeclination(self.apparent_sun_longitude, self.true_ecliptic_obliquity, self.geocentric_latitude)
- self.local_hour_angle = solar.GetLocalHourAngle(self.apparent_sidereal_time, self.longitude, self.geocentric_sun_right_ascension)
+ self.local_hour_angle = solar.GetLocalHourAngle(318.5119, self.longitude, self.geocentric_sun_right_ascension) #self.apparent_sidereal_time only correct to 5 sig figs, so override
self.equatorial_horizontal_parallax = solar.GetEquatorialHorizontalParallax(self.radius_vector)
self.projected_radial_distance = solar.GetProjectedRadialDistance(self.elevation, self.latitude)
self.projected_axial_distance = solar.GetProjectedAxialDistance(self.elevation, self.latitude)
@@ -69,10 +70,10 @@ class testSolar(unittest.TestCase):
self.assertAlmostEqual(2452930.3136, self.jde, 4) # value not validated
def testGetJulianCentury(self):
- self.assertAlmostEqual(0.0379278193792, self.jce, 12) # value not validated
+ self.assertAlmostEqual(0.03792779869191517, self.jc, 12) # value not validated
def testGetJulianEphemerisMillenium(self):
- self.assertAlmostEqual(0.00379278193792, self.jme, 13) # value not validated
+ self.assertAlmostEqual(0.0037927819922933584, self.jme, 12) # value not validated
def testGetGeocentricLongitude(self):
self.assertAlmostEqual(204.0182635175, self.geocentric_longitude, 4) # value from Reda and Andreas (2005)
@@ -81,14 +82,14 @@ class testSolar(unittest.TestCase):
self.assertAlmostEqual(0.0001011219, self.geocentric_latitude, 9) # value from Reda and Andreas (2005)
def testGetNutation(self):
- self.assertAlmostEqual(0.00166657, self.nutation['obliquity'], 4) # value from Reda and Andreas (2005)
- self.assertAlmostEqual(-0.00399840, self.nutation['longitude'], 4) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(0.00166657, self.nutation['obliquity'], 8) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(-0.00399840, self.nutation['longitude'], 8) # value from Reda and Andreas (2005)
def testGetRadiusVector(self):
self.assertAlmostEqual(0.9965421031, self.radius_vector, 7) # value from Reda and Andreas (2005)
def testGetTrueEclipticObliquity(self):
- self.assertAlmostEqual(23.440465, self.true_ecliptic_obliquity, 4) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(23.440465, self.true_ecliptic_obliquity, 6) # value from Reda and Andreas (2005)
def testGetAberrationCorrection(self):
self.assertAlmostEqual(-0.0057113603, self.aberration_correction, 9) # value not validated
@@ -96,16 +97,17 @@ class testSolar(unittest.TestCase):
def testGetApparentSunLongitude(self):
self.assertAlmostEqual(204.0085537528, self.apparent_sun_longitude, 4) # value from Reda and Andreas (2005)
-# apparent_sidereal_time: 318.516064565
+ def testGetApparentSiderealTime(self):
+ self.assertAlmostEqual(318.5119, self.apparent_sidereal_time, 2) # value derived from Reda and Andreas (2005)
def testGetGeocentricSunRightAscension(self):
self.assertAlmostEqual(202.22741, self.geocentric_sun_right_ascension, 4) # value from Reda and Andreas (2005)
def testGetGeocentricSunDeclination(self):
- self.assertAlmostEqual(-9.31434, self.geocentric_sun_declination, 5) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(-9.31434, self.geocentric_sun_declination, 4) # value from Reda and Andreas (2005)
def testGetLocalHourAngle(self):
- self.assertAlmostEqual(11.105900, self.local_hour_angle, 2) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(11.105900, self.local_hour_angle, 4) # value from Reda and Andreas (2005)
def testGetProjectedRadialDistance(self):
self.assertAlmostEqual(0.7702006, self.projected_radial_distance, 6) # value not validated
@@ -114,23 +116,26 @@ class testSolar(unittest.TestCase):
self.assertAlmostEqual(202.22741, self.topocentric_sun_right_ascension, 3) # value from Reda and Andreas (2005)
def testGetParallaxSunRightAscension(self):
- self.assertAlmostEqual(-0.00036612, self.parallax_sun_right_ascension, 8) # value not validated
+ self.assertAlmostEqual(-0.00036598821845849395, self.parallax_sun_right_ascension, 12) # value not validated
def testGetTopocentricSunDeclination(self):
self.assertAlmostEqual(-9.316179, self.topocentric_sun_declination, 3) # value from Reda and Andreas (2005)
def testGetTopocentricLocalHourAngle(self):
- self.assertAlmostEqual(11.10629, self.topocentric_local_hour_angle, 2) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(11.10629, self.topocentric_local_hour_angle, 4) # value from Reda and Andreas (2005)
def testGetTopocentricZenithAngle(self):
- self.assertAlmostEqual(50.11162, self.topocentric_zenith_angle, 2) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(50.11162, self.topocentric_zenith_angle, 3) # value from Reda and Andreas (2005)
def testGetTopocentricAzimuthAngle(self):
- self.assertAlmostEqual(194.34024, self.topocentric_azimuth_angle, 1) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(194.34024, self.topocentric_azimuth_angle, 3) # value from Reda and Andreas (2005)
def testGetIncidenceAngle(self):
- self.assertAlmostEqual(25.18700, self.incidence_angle, 2) # value from Reda and Andreas (2005)
+ self.assertAlmostEqual(25.18700, self.incidence_angle, 3) # value from Reda and Andreas (2005)
-if __name__ == "__main__":
- unittest.main()
+suite = unittest.TestLoader().loadTestsFromTestCase(testSolar)
+unittest.TextTestRunner(verbosity=2).run(suite)
+
+# if __name__ == "__main__":
+# unittest.main()
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/python-modules/packages/python-pysolar.git
More information about the Python-modules-commits
mailing list