Merge pull request #10 from seanth/master

Calculate pressure and temperature from elevation

Pysolar/elevation.py

@@ -0,0 +1,62 @@
+#!/usr/bin/python
+
+#    Copyright Sean T. Hammond
+#
+#    This file is part of Pysolar.
+#
+#    Pysolar is free software; you can redistribute it and/or modify
+#    it under the terms of the GNU General Public License as published by
+#    the Free Software Foundation; either version 3 of the License, or
+#    (at your option) any later version.
+#
+#    Pysolar is distributed in the hope that it will be useful,
+#    but WITHOUT ANY WARRANTY; without even the implied warranty of
+#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+#    GNU General Public License for more details.
+#
+#    You should have received a copy of the GNU General Public License along
+#    with Pysolar. If not, see <http://www.gnu.org/licenses/>.
+
+"""Various elevation-related calculations
+
+"""
+import math
+
+def GetPressureWithElevation(h, Ps=101325.00, Ts=288.15, Tl=-0.0065, Hb=0.0, R=8.31432, g=9.80665, M=0.0289644):
+	#This function returns an estimate of the pressure in pascals as a function of elevation above sea level
+	#NOTE: This equation is only accurate up to 11,000 meters
+	#NOTE: results might be odd for elevations below 0 (sea level), like Dead Sea.
+	#h=elevation relative to sea level (m)
+	#Ps= static pressure (pascals) = 101325.00 P
+	#Ts= standard temperature (kelvin) = 288.15 K
+	#Tl= temperature lapse rate (kelvin/meter) = -0.0065 K/m
+	#Hb= height at the bottom of the layer = 0
+	#R= universal gas constant for air = 8.31432 N*m/s^2
+	#g= gravitational acceleration for earth = 9.80665 m/s^2
+	#M= Molar mass of Earth's atmosphere = 0.0289644 kg/mol
+	#P=Ps*(Ts/((Ts+Tl)*(h-Hb)))^((g*M)/(R*Tl))
+	#returns pressure in pascals
+	if h>11000.0: print"WARNING: Elevation used exceeds the recommended maximum elevation for this function (11,000m)"
+	theDenominator = Ts+(Tl*(h-Hb))
+	theExponent=(g*M)/(R*Tl)
+	return Ps*(Ts/theDenominator)**theExponent
+
+def GetTemperatureWithElevation(h, Ts=288.15, Tl=-0.0065):
+	#This function returns an estimate of temperature as a function above sea level
+	#NOTE: this is only accurate up to 11,000m
+	#NOTE: results might be odd for elevations below 0 (sea level), like Dead Sea.
+	#Ts= standard temperature (kelvin) = 288.15 K
+	#Tl= temperature lapse rate (kelvin/meter) = -0.0065 K/m
+	#returns temp in kelvin
+	return Ts+(h*Tl)
+
+def ElevationTest():
+	print "Elevation(m) Pressure(Pa) Temperature(K)"
+	h=0
+	for i in range(11):
+		P=GetPressureWithElevation(h)
+		T=GetTemperatureWithElevation(h)
+		print "%i %i %i" % (h, P, T)
+		h=h+1000
+
+

Pysolar/testsolar.py

@@ -22,6 +22,7 @@
 import solar
 import constants
 import julian
+import elevation
 import datetime
 import unittest
 
@@ -56,13 +57,15 @@ def setUp(self):
 		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.apparent_sun_longitude, self.true_ecliptic_obliquity, self.geocentric_latitude)
+		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)
 		self.topocentric_azimuth_angle = solar.GetTopocentricAzimuthAngle(self.topocentric_local_hour_angle, self.latitude, self.topocentric_sun_declination)
 		self.incidence_angle = solar.GetIncidenceAngle(self.topocentric_zenith_angle, self.slope, self.slope_orientation, self.topocentric_azimuth_angle)
+		self.pressure_with_elevation = elevation.GetPressureWithElevation(1567.7)
+		self.temperature_with_elevation = elevation.GetTemperatureWithElevation(1567.7)
 
 	def testGetJulianDay(self):
 		self.assertAlmostEqual(2452930.312847, self.jd, 6) # value from Reda and Andreas (2005)
@@ -134,6 +137,12 @@ def testGetTopocentricAzimuthAngle(self):
 	def testGetIncidenceAngle(self):
 		self.assertAlmostEqual(25.18700, self.incidence_angle, 3) # value from Reda and Andreas (2005)
 
+	def testPressureWithElevation(self):
+		self.assertAlmostEqual(83855.90228, self.pressure_with_elevation, 4)
+
+	def testTemperatureWithElevation(self):
+		self.assertAlmostEqual(277.9600, self.temperature_with_elevation, 4)
+
 suite = unittest.TestLoader().loadTestsFromTestCase(testSolar)
 unittest.TextTestRunner(verbosity=2).run(suite)