@@ -188,4 +188,4 @@ if __name__ == '__main__': print 'Min error: ' + str(stats.tmin(alt_errors, None)) print 'Max error: ' + str(stats.tmax(alt_errors, None)) -# WriteComparisonsToCSV(comps, 'pysolar_v_usno.csv') + WriteComparisonsToCSV(comps, 'pysolar_v_usno.csv') query_usno.py @@ -73,7 +73,7 @@ def GetAltitude(latitude_deg, longitude_deg, utc_datetime, elevation = 0, temper geocentric_sun_right_ascension = GetGeocentricSunRightAscension(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude) geocentric_sun_declination = GetGeocentricSunDeclination(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude) local_hour_angle = GetLocalHourAngle(apparent_sidereal_time, longitude_deg, geocentric_sun_right_ascension) - parallax_sun_right_ascension = GetParallaxSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, geocentric_sun_declination, projected_axial_distance) + parallax_sun_right_ascension = GetParallaxSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, geocentric_sun_declination) topocentric_local_hour_angle = GetTopocentricLocalHourAngle(local_hour_angle, parallax_sun_right_ascension) topocentric_sun_declination = GetTopocentricSunDeclination(geocentric_sun_declination, projected_axial_distance, equatorial_horizontal_parallax, parallax_sun_right_ascension, local_hour_angle) topocentric_elevation_angle = GetTopocentricElevationAngle(latitude_deg, topocentric_sun_declination, topocentric_local_hour_angle) @@ -124,7 +124,7 @@ def GetAzimuth(latitude_deg, longitude_deg, utc_datetime, elevation = 0): geocentric_sun_right_ascension = GetGeocentricSunRightAscension(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude) geocentric_sun_declination = GetGeocentricSunDeclination(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude) local_hour_angle = GetLocalHourAngle(apparent_sidereal_time, longitude_deg, geocentric_sun_right_ascension) - parallax_sun_right_ascension = GetParallaxSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, geocentric_sun_declination, projected_axial_distance) + parallax_sun_right_ascension = GetParallaxSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, geocentric_sun_declination) topocentric_local_hour_angle = GetTopocentricLocalHourAngle(local_hour_angle, parallax_sun_right_ascension) topocentric_sun_declination = GetTopocentricSunDeclination(geocentric_sun_declination, projected_axial_distance, equatorial_horizontal_parallax, parallax_sun_right_ascension, local_hour_angle) return 180 - GetTopocentricAzimuthAngle(topocentric_local_hour_angle, latitude_deg, topocentric_sun_declination) @@ -145,10 +145,7 @@ def GetAzimuthFast(latitude_deg, longitude_deg, utc_datetime): return (180 - math.degrees(azimuth_rad)) def GetCoefficient(jme, constant_array): - total = 0 - for i in range(len(constant_array)): - total = total + (constant_array[i-1][0] * math.cos(constant_array[i-1][1] + (constant_array[i-1][2] * jme))) - return total + return sum([constant_array[i-1][0] * math.cos(constant_array[i-1][1] + (constant_array[i-1][2] * jme)) for i in range(len(constant_array))]) def GetDayOfYear(utc_datetime): year_start = datetime.datetime(utc_datetime.year, 1, 1,) @@ -261,14 +258,13 @@ def GetNutation(jde): return nutation -def GetParallaxSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, geocentric_sun_declination, projected_axial_distance): +def GetParallaxSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, geocentric_sun_declination): prd = projected_radial_distance ehp_rad = math.radians(equatorial_horizontal_parallax) lha_rad = math.radians(local_hour_angle) gsd_rad = math.radians(geocentric_sun_declination) - pad = projected_axial_distance a = -1 * prd * math.sin(ehp_rad) * math.sin(lha_rad) - b = math.cos(gsd_rad) - pad * math.sin(ehp_rad) * math.cos(lha_rad) + b = math.cos(gsd_rad) - prd * math.sin(ehp_rad) * math.cos(lha_rad) parallax = math.atan2(a, b) return math.degrees(parallax) @@ -326,14 +322,16 @@ def GetTopocentricSunDeclination(geocentric_sun_declination, projected_axial_dis gsd_rad = math.radians(geocentric_sun_declination) pad = projected_axial_distance ehp_rad = math.radians(equatorial_horizontal_parallax) - a = (math.sin(gsd_rad) - pad * math.sin(ehp_rad)) * math.cos(parallax_sun_right_ascension) - b = math.cos(gsd_rad) - (pad * math.sin(ehp_rad) * math.cos(local_hour_angle)) + psra_rad = math.radians(parallax_sun_right_ascension) + lha_rad = math.radians(local_hour_angle) + a = (math.sin(gsd_rad) - pad * math.sin(ehp_rad)) * math.cos(psra_rad) + b = math.cos(gsd_rad) - (pad * math.sin(ehp_rad) * math.cos(lha_rad)) return math.degrees(math.atan2(a, b)) -def GetTopocentricSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, projected_axial_distance, +def GetTopocentricSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude): gsd = GetGeocentricSunDeclination(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude) - psra = GetParallaxSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, gsd, projected_axial_distance) + psra = GetParallaxSunRightAscension(projected_radial_distance, equatorial_horizontal_parallax, local_hour_angle, gsd) gsra = GetGeocentricSunRightAscension(apparent_sun_longitude, true_ecliptic_obliquity, geocentric_latitude) return psra + gsra solar.py @@ -56,8 +56,8 @@ class testSolar(unittest.TestCase): self.projected_radial_distance = solar.GetProjectedRadialDistance(self.elevation, self.latitude) self.projected_axial_distance = solar.GetProjectedAxialDistance(self.elevation, self.latitude) self.topocentric_sun_right_ascension = solar.GetTopocentricSunRightAscension(self.projected_radial_distance, - self.equatorial_horizontal_parallax, self.local_hour_angle, self.projected_axial_distance, self.apparent_sun_longitude, self.true_ecliptic_obliquity, self.geocentric_latitude) - self.parallax_sun_right_ascension = solar.GetParallaxSunRightAscension(self.projected_radial_distance, self.equatorial_horizontal_parallax, self.local_hour_angle, self.geocentric_sun_declination, self.projected_axial_distance) + self.equatorial_horizontal_parallax, self.local_hour_angle, self.apparent_sun_longitude, self.true_ecliptic_obliquity, self.geocentric_latitude) + self.parallax_sun_right_ascension = solar.GetParallaxSunRightAscension(self.projected_radial_distance, self.equatorial_horizontal_parallax, self.local_hour_angle, self.geocentric_sun_declination) self.topocentric_sun_declination = solar.GetTopocentricSunDeclination(self.geocentric_sun_declination, self.projected_axial_distance, self.equatorial_horizontal_parallax, self.parallax_sun_right_ascension, self.local_hour_angle) self.topocentric_local_hour_angle = solar.GetTopocentricLocalHourAngle(self.local_hour_angle, self.parallax_sun_right_ascension) self.topocentric_zenith_angle = solar.GetTopocentricZenithAngle(self.latitude, self.topocentric_sun_declination, self.topocentric_local_hour_angle, self.pressure, self.temperature) @@ -117,7 +117,7 @@ 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.00036598821845849395, self.parallax_sun_right_ascension, 12) # value not validated + self.assertAlmostEqual(-0.00036599029186055283, 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) @@ -129,7 +129,7 @@ class testSolar(unittest.TestCase): 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, 3) # value from Reda and Andreas (2005) + self.assertAlmostEqual(194.34024, self.topocentric_azimuth_angle, 5) # value from Reda and Andreas (2005) def testGetIncidenceAngle(self): self.assertAlmostEqual(25.18700, self.incidence_angle, 3) # value from Reda and Andreas (2005) testsolar.py 36f7d48a45789e8e5d5e47525808db4332530f6b Brandon Stafford brandon@pingswept.org http://github.com/pingswept/pysolar/commit/f828ce01954571a4908c867148578079334752d9 f828ce01954571a4908c867148578079334752d9 2009-01-11T17:34:39-08:00 2009-01-11T17:34:39-08:00 Added corrections from Sean Taylor to GetParallaxSunRightAscension() and GetTopocentricSunDeclination() and updated test suite. Variation against USNO actually increased slightly, but Sean Taylor's changes are still right. The differences relative to the USNO are dominated by two data points, 2008-Mar-06 3:53:58 (largest variation in altitude relative to USNO) and 2008-Feb-10 21:09:53 (largest variation in azimuth relative to USNO). dcd51c4635c875a656799c68dd486f252f555be9 Brandon Stafford brandon@pingswept.org