@@ -1,4 +1,4 @@ -#!/usr/bin/python + #!/usr/bin/python # Julian calendar calculations for calculating the position of the sun relative to the earth @@ -42,7 +42,7 @@ def GetJulianDay(utc_datetime): # based on NREL/TP-560-34302 by Andreas and Reda def GetJulianEphemerisCentury(julian_ephemeris_day): return (julian_ephemeris_day - 2451545.0) / 36525.0 -def GetJulianEphemerisDay(julian_day, delta_seconds): +def GetJulianEphemerisDay(julian_day, delta_seconds = 66.0): """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. A list of values for Delta-T can be found here: ftp://maia.usno.navy.mil/ser7/deltat.data""" julian.py @@ -19,6 +19,9 @@ # You should have received a copy of the GNU General Public License along # with Pysolar. If not, see <http://www.gnu.org/licenses/>. +import solar +import math + def GetAirMassRatio(altitude_deg): # from Masters, p. 412 # warning: pukes on input of zero @@ -34,8 +37,12 @@ def GetOpticalDepth(day): def GetRadiationDirect(utc_datetime, altitude_deg): # from Masters, p. 412 - day = GetDayOfYear(utc_datetime) - flux = GetApparentExtraterrestrialFlux(day) - optical_depth = GetOpticalDepth(day) - air_mass_ratio = GetAirMassRatio(altitude_deg) - return flux * math.exp(-1 * optical_depth * air_mass_ratio) + + if(altitude_deg > 0): + day = solar.GetDayOfYear(utc_datetime) + flux = GetApparentExtraterrestrialFlux(day) + optical_depth = GetOpticalDepth(day) + air_mass_ratio = GetAirMassRatio(altitude_deg) + return flux * math.exp(-1 * optical_depth * air_mass_ratio) + else: + return 0.0 radiation.py @@ -20,10 +20,11 @@ # with Pysolar. If not, see <http://www.gnu.org/licenses/>. import datetime +import radiation import solar def BuildTimeList(start_utc_datetime, end_utc_datetime, step_minutes): - + '''Create a list of sample points evenly spaced apart by step_minutes.''' step = step_minutes * 60 time_list = [] span = end_utc_datetime - start_utc_datetime @@ -31,9 +32,26 @@ def BuildTimeList(start_utc_datetime, end_utc_datetime, step_minutes): print span return map(lambda n: start_utc_datetime + dt * n, range((span.days * 86400 + span.seconds) / step)) -def SimulateSpan(latitude_deg, longitude_deg, start_utc_datetime, end_utc_datetime, step_minutes): +def SimulateSpan(latitude_deg, longitude_deg, start_utc_datetime, end_utc_datetime, step_minutes, elevation = 0, temperature_celsius = 25, pressure_millibars = 1013.25): + '''Simulate the motion of the sun over a time span and location of your choosing. - time_list = buildTimeList(start_utc_datetime, end_utc_datetime, step_minutes) + The start and end points are set by datetime objects, which can be created with + the standard Python datetime module like this: + import datetime + start = datetime.datetime(2008, 12, 23, 23, 14, 0) + ''' + time_list = BuildTimeList(start_utc_datetime, end_utc_datetime, step_minutes) - for time in time_list: - print 'Altitude: ' + str(solar.GetAltitude(latitude_deg, longitude_deg, time)) + ' Azimuth: ' + str(solar.GetAzimuth(latitude_deg, longitude_deg, time)) + angles_list = [( + time, + solar.GetAltitude(latitude_deg, longitude_deg, time, elevation, temperature_celsius, pressure_millibars), + solar.GetAzimuth(latitude_deg, longitude_deg, time, elevation) + ) for time in time_list] + power_list = [(time, alt, az, radiation.GetRadiationDirect(time, alt)) for (time, alt, az) in angles_list] + print power_list + +# xs = shade.GetXShade(width, 120, azimuth_deg) +# ys = shade.GetYShade(height, 120, altitude_deg) +# shaded_area = xs * ys +# shaded_percentage = shaded_area/area +# import simulate, datetime; s = datetime.datetime(2008,1,1); e = datetime.datetime(2008,1,5); simulate.SimulateSpan(42.0, -70.0, s, e, 30) simulate.py @@ -45,10 +45,42 @@ def EquationOfTime(day): b = (2 * math.pi / 364.0) * (day - 81) return (9.87 * math.sin(2 *b)) - (7.53 * math.cos(b)) - (1.5 * math.sin(b)) -def GetAberrationCorrection(r): # r is earth radius vector [astronomical units] - return -20.4898/(3600.0 * r) - -def GetAltitude(latitude_deg, longitude_deg, utc_datetime): +def GetAberrationCorrection(radius_vector): # r is earth radius vector [astronomical units] + return -20.4898/(3600.0 * radius_vector) + +def GetAltitude(latitude_deg, longitude_deg, utc_datetime, elevation = 0, temperature_celsius = 25, pressure_millibars = 1013.25): + '''See also the faster, but less accurate, GetAltitudeFast()''' + # location-dependent calculations + projected_radial_distance = GetProjectedRadialDistance(elevation, latitude_deg) + projected_axial_distance = GetProjectedAxialDistance(elevation, latitude_deg) + + # time-dependent calculations + jd = julian.GetJulianDay(utc_datetime) + jde = julian.GetJulianEphemerisDay(jd, 65) + jce = julian.GetJulianEphemerisCentury(jde) + jme = julian.GetJulianEphemerisMillenium(jce) + geocentric_latitude = GetGeocentricLatitude(jme) + geocentric_longitude = GetGeocentricLongitude(jme) + radius_vector = GetRadiusVector(jme) + aberration_correction = GetAberrationCorrection(radius_vector) + equatorial_horizontal_parallax = GetEquatorialHorizontalParallax(radius_vector) + nutation = GetNutation(jde) + apparent_sidereal_time = GetApparentSiderealTime(jd, jme, nutation) + true_ecliptic_obliquity = GetTrueEclipticObliquity(jme, nutation) + + # calculations dependent on location and time + apparent_sun_longitude = GetApparentSunLongitude(geocentric_longitude, nutation, aberration_correction) + 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) + 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) + refraction_correction = GetRefractionCorrection(pressure_millibars, temperature_celsius, topocentric_elevation_angle) + return topocentric_elevation_angle + refraction_correction + +def GetAltitudeFast(latitude_deg, longitude_deg, utc_datetime): # expect 19 degrees for solar.GetAltitude(42.364908,-71.112828,datetime.datetime(2007, 2, 18, 20, 13, 1, 130320)) @@ -67,7 +99,37 @@ def GetApparentSiderealTime(julian_day, jme, nutation): def GetApparentSunLongitude(geocentric_longitude, nutation, ab_correction): return geocentric_longitude + nutation['longitude'] + ab_correction -def GetAzimuth(latitude_deg, longitude_deg, utc_datetime): +def GetAzimuth(latitude_deg, longitude_deg, utc_datetime, elevation = 0): + + # location-dependent calculations + projected_radial_distance = GetProjectedRadialDistance(elevation, latitude_deg) + projected_axial_distance = GetProjectedAxialDistance(elevation, latitude_deg) + + # time-dependent calculations + jd = julian.GetJulianDay(utc_datetime) + jde = julian.GetJulianEphemerisDay(jd, 65) + jce = julian.GetJulianEphemerisCentury(jde) + jme = julian.GetJulianEphemerisMillenium(jce) + geocentric_latitude = GetGeocentricLatitude(jme) + geocentric_longitude = GetGeocentricLongitude(jme) + radius_vector = GetRadiusVector(jme) + aberration_correction = GetAberrationCorrection(radius_vector) + equatorial_horizontal_parallax = GetEquatorialHorizontalParallax(radius_vector) + nutation = GetNutation(jde) + apparent_sidereal_time = GetApparentSiderealTime(jd, jme, nutation) + true_ecliptic_obliquity = GetTrueEclipticObliquity(jme, nutation) + + # calculations dependent on location and time + apparent_sun_longitude = GetApparentSunLongitude(geocentric_longitude, nutation, aberration_correction) + 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) + 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) + +def GetAzimuthFast(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)) day = GetDayOfYear(utc_datetime) declination_rad = math.radians(GetDeclination(day)) @@ -94,11 +156,11 @@ def GetDayOfYear(utc_datetime): return delta.days def GetDeclination(day): - """The declination of the sun is the angle between + '''The declination of the sun is the angle between Earth's equatorial plane and a line between the Earth and the sun. The declination of the sun varies between 23.45 degrees and -23.45 degrees, hitting zero on the equinoxes and peaking on the solstices. - """ + ''' return 23.45 * math.sin((2 * math.pi / 365.0) * (day - 81)) def GetEquatorialHorizontalParallax(radius_vector): solar.py cb121fc2a107c5742eff64fe6580e3f579cbc75c Brandon Stafford brandon@pingswept.org http://github.com/pingswept/pysolar/commit/642720ed3dcfed8be39113809f4207d691fdf05d 642720ed3dcfed8be39113809f4207d691fdf05d 2008-03-16T13:14:10-07:00 2008-03-16T13:14:10-07:00 Fixed some broken stuff in radiation.py. Advanced simulation.py a little. 22e650da5100058688aea79fe4b047e5447c30e3 Brandon Stafford brandon@pingswept.org