Solar Geometry module: save into a file called solarGeom.py
import datetime import math #converts a date given as "mo/day" into a day of year def calc_dayOfYear(dateIn): daysInMonth= [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31] (month, day) = dateIn.split('/') dayOut = 0; for m in range(1,int(month)): dayOut = dayOut + daysInMonth[m-1] dayOut = dayOut + int(day) return dayOut #converts a time given in 24h time "hour:minutes" into a decimal hour def calc_hourDecimal(timeIn): (hour, minutes) = timeIn.split(':') return int(hour) + int(minutes)/60.0 #calculates alpha = ecliptic longitute, varying from alpha = 0 at the March #Equinox completing 360 in 1 year. Divides the year into intervals between s #solstice-equinox and treats alpha as uniformly varying throughout each interval #Output in degrees. def calc_alpha(dayIn): marchEquinox = calc_dayOfYear("3/20") juneSolstice = calc_dayOfYear("6/20") septemberEquinox = calc_dayOfYear("9/22") decemberSolstice = calc_dayOfYear("12/21") if (marchEquinox<=dayIn <= juneSolstice) : alphaOut = 90*(dayIn-marchEquinox)/(juneSolstice-marchEquinox) elif (juneSolstice < dayIn <=septemberEquinox): alphaOut = 90+90*(dayIn-juneSolstice)/(septemberEquinox-juneSolstice) elif (septemberEquinox< dayIn<=decemberSolstice): alphaOut = 180+90*(dayIn-septemberEquinox)/(decemberSolstice-septemberEquinox) else: alphaOut = 270+90*(dayIn+365-decemberSolstice)/(marchEquinox+365-decemberSolstice) return alphaOut #calculates the following solar position angles for given coordinates, #day of the year, local time. # declination # hour angle # altitude # azimuth # All output in degrees def calc_solarAngles(latitudeIn, longitudeIn, timezone, dayIn, hourIn): alpha = calc_alpha(dayIn) alpha_rad = math.radians(alpha) obliquity_rad = math.radians(23.44) #calc delta = declination angle delta_rad = math.asin(math.sin(alpha_rad)*math.sin(obliquity_rad)) delta = math.degrees(delta_rad) #calc omega = hour angle, angle between local longitude and solar noon ET = (0.0172 + 0.4281*math.cos(alpha_rad)-7.3515*math.sin(alpha_rad) \ -3.3495*math.cos(2*alpha_rad) - 9.3619*math.sin(2*alpha_rad))/60 omega = 15*(hourIn + longitudeIn/15 - timezone + ET - 12) omega_rad = math.radians(omega) #calc altitude phi_rad = math.radians(latitudeIn) altitude_rad = math.asin(math.cos(delta_rad)*math.cos(phi_rad)*math.cos(omega_rad) + \ math.sin(delta_rad)*math.sin(phi_rad)) altitude = math.degrees(altitude_rad) #calc azimuth azimuth_rad = math.acos((math.sin(delta_rad)*math.cos(phi_rad) - \ math.cos(delta_rad)*math.sin(phi_rad)*math.cos(omega_rad))/math.cos(altitude_rad)) if math.degrees(omega_rad) < 0: azimuth = math.degrees(azimuth_rad) else: azimuth = 360-math.degrees(azimuth_rad) return [delta, omega, altitude, azimuth] #determines the sun vector for given coordinates, day of the year, local time. def calc_sunVector(latitudeIn, longitudeIn, timezone, dayIn, hourIn): [d, o, altitude, azimuth] = calc_solarAngles(latitudeIn, longitudeIn, timezone, dayIn, hourIn) altitude_rad= math.radians(altitude) azimuth_rad = math.radians(azimuth) x = math.cos(altitude_rad)*math.sin(azimuth_rad) y = math.cos(altitude_rad)*math.cos(azimuth_rad) z = math.sin(altitude_rad) return [x, y, z] # Check to see if this file is being executed as the "Main" python # script instead of being used as a module by some other python script # This allows us to use the module which ever way we want. if __name__ == '__main__': print "Welcome to the solar geometry module!"
Example usage: To use the solarGeom module, you will need to first import it
import solarGeom as sg
This module can be used to compute the solar angles and the sun vector for any day and time (given in local time) for any location given in terms of coordinates and time zone relative to GMT.
import rhinoscriptsyntax as rs latitude = 40.75 longitude = -73.5 timezone = -5.0 day = sg.calc_dayOfYear("08/03") hour = sg.calc_hourDecimal("10:40") solarAngles = sg.calc_solarAngles(latitude, longitude, timezone, day, hour) print "Altitude = ", solarAngles[2] print "Azimuth = ", solarAngles[3] sun = sg.calc_sunVector(latitude, longitude, timezone, day, hour) R = 20 #scaling for the sun vector so that it is more visible on canvas sun[:] = [x*R for x in sun] rs.AddPoints(sun)
This can also be used to show the sun path throughout any range of times in the day.
startTime= "7:00" endTime = "19:30" hourStart = sg.calc_hourDecimal(startTime) hourEnd = sg.calc_hourDecimal(endTime) R = 20 N = 20 #number of sample points between starting time and end time for i in range(N+1): sun = sg.calc_sunVector(latitude, longitude, timezone, day, hourStart+i*(hourEnd-hourStart)/N) if sun[2] > 0: sun[:] = [x*R for x in sun] rs.AddPoint(sun)
This can also be used to show the sun path throughout a range of dates in the year, and for each date, a range of times in the day. It is interesting, and quite beautiful to see the sun's path traced out relative to a fixed time throughout the entire year.
dateStart = "1/1" dateEnd = "12/31" dayStart = sg.calc_dayOfYear(dateStart) dayEnd = sg.calc_dayOfYear(dateEnd) startTime= "7:00" endTime = "19:30" hourStart = sg.calc_hourDecimal(startTime) hourEnd = sg.calc_hourDecimal(endTime) R = 20 N = 20 #number of sample points between starting time and end time for dayNow in range(dayStart, dayEnd, 3): for i in range(N+1): sun = sg.calc_sunVector(latitude, longitude, timezone, dayNow, hourStart+i*(hourEnd-hourStart)/N) if sun[2] > 0: sun[:] = [x*R for x in sun] rs.AddPoint(sun)