__author__ = "Jwely"
from datetime import datetime, timedelta
from numpy import radians, ndarray, sin, cos, degrees, arctan2, arcsin, tan, arccos
[docs]class solar:
"""
Object class for handling solar calculations. Many equations are taken from the
excel sheet at this url : [http://www.esrl.noaa.gov/gmd/grad/solcalc/calcdetails.html]
It requires a physical location on the earth and a datetime object
:param lat: decimal degrees latitude (float OR numpy array)
:param lon: decimal degrees longitude (float OR numpy array)
:param time_zone: float of time shift from GMT (such as "-5" for EST)
:param date_time_obj: either a timestamp string following fmt or a datetime obj
:param fmt: if date_time_obj is a string, fmt is required to interpret it
:param slope: slope of land at lat,lon for solar energy calculations
:param aspect: aspect of land at lat,lon for solar energy calculations
An instance of this class may have the following attributes:
=================== =========================================== ========
attribute description type
=================== =========================================== ========
lat latitude (array)
lon longitude (array)
tz time zone (scalar)
rdt reference datetime object (date_time_obj) (scalar)
slope slope, derivative of DEM (array)
aspect aspect (north is 0, south is 180) (array)
ajd absolute julian day (scalar)
ajc absolute julian century (scalar)
geomean_long geometric mean longitude of the sun (scalar)
geomean_anom geometric mean longitude anomaly of the sun (scalar)
earth_eccent eccentricity of earths orbit (scalar)
sun_eq_of_center the suns equation of center (scalar)
true_long true longitude of the sun (scalar)
true_anom true longitude anomaly of the sun (scalar)
app_long the suns apparent longitude (scalar)
oblique_mean_elip earth oblique mean ellipse (scalar)
oblique_corr correction to earths oblique elipse (scalar)
right_ascension suns right ascension angle (scalar)
declination solar declination angle (scalar)
equation_of_time equation of time (minutes) (scalar)
hour_angle_sunrise the hour angle at sunrise (array)
solar_noon LST of solar noon (array)
sunrise LST of sunrise time (array)
sunset LST of sunset time (array)
sunlight LST fractional days of sunlight (array)
true_solar LST for true solar time (array)
hour_angle total hour angle (array)
zenith zenith angle (array)
elevation elevation angle (array)
azimuth azimuthal angle (array)
rad_vector radiation vector (distance in AU) (scalar)
earth_distance earths distance to sun in meters (scalar)
norm_irradiance incident solar energy at earth distance (scalar)
=================== =========================================== ========
Units used by this class unless otherwise labeled
- angle = degrees
- distance = meters
- energy = watts or joules
- time = mostly in datetime objects. labeled in most cases.
Planned improvements
1. DONE. Inputs of numpy arrays for lat and lon needs to be allowed.
2. inputs of a numpy array DEM for slope/aspect effects on incident solar energy
Present performance
To process about one landsat tile (7300^2 matrix) requires 9GB of memory and
takes 45 seconds to process on a single 3.3GHz thread. It would be nice to get
the same output to run on ~5GB of memory so a 8GB system could handle it.
"""
def __init__(self, lat, lon, date_time_obj, time_zone = 0,
fmt = False, slope = None, aspect = None):
"""
Initializes critical spatial and temporal information for solar object.
"""
# empty list of class attributes
self.ajc = None # abs julian century (defined on __init__)
self.ajd = None # abs julian day (defined on __init__)
self.app_long = None
self.atmo_refraction = None
self.azimuth = None
self.declination = None
self.earth_distance = None
self.earth_eccent = None
self.elevation = None
self.elevation_noatmo = None
self.equation_of_time = None
self.frac_day = None
self.geomean_anom = None
self.geomean_long = None
self.hour_angle = None
self.hour_angle_sunrise = None
self.lat = lat # lattitude (E positive)- float
self.lat_r = radians(lat) # lattitude in radians
self.lon = lon # longitude (N positive)- float
self.lon_r = radians(lon) # longitude in radians
self.norm_irradiance = None
self.oblique_corr = None
self.oblique_mean_elip = None
self.rad_vector = None
self.rdt = None # reference datetime (defined on __init__)
self.right_ascension = None
self.solar_noon = None
self.solar_noon_time = None
self.sun_eq_of_center = None
self.sunlight = None
self.sunlight_time = None
self.sunrise = None
self.sunrise_time = None
self.sunset = None
self.sunset_time = None
self.true_anom = None
self.true_long = None
self.true_solar = None
self.true_solar_time = None
self.tz = None # time zone (defined on __init__)
self.zenith = None
# slope and aspect
self.slope = slope
self.aspect = aspect
# Constants as attributes
self.sun_surf_rad = 63156942.6 # radiation at suns surface (W/m^2)
self.sun_radius = 695800000. # radius of the sun in meters
self.orbital_period = 365.2563630 # num of days it takes earth to revolve
self.altitude = -0.01448623 # altitude of center of solar disk
# sets up the object with some subfunctions
self._set_datetime(date_time_obj, fmt, GMT_hour_offset = time_zone)
# specify if attributes are scalar floats or numpy array floats
if isinstance(lat, ndarray) and isinstance(lon, ndarray):
self.is_numpy = True
else:
self.is_numpy = False
return
def _set_datetime(self, date_time_obj, fmt = False, GMT_hour_offset = 0):
"""
sets the critical time information including absolute julian day/century.
Accepts datetime objects or a datetime string with format
:param date_time_obj: datetime object for time of solar calculations. Will also
accept string input with matching value for "fmt" param
:param fmt: if date_time_obj is input as a string, fmt allows it to be
interpreted
:param GMT_hour_offset: Number of hours from GMT for timezone of calculation area.
"""
# if input is datetime_obj set it
if isinstance(date_time_obj, datetime):
self.rdt = date_time_obj
self.rdt += timedelta(hours = -GMT_hour_offset)
elif isinstance(date_time_obj, str) and isinstance(fmt, str):
self.rdt = datetime.strptime(date_time_obj,fmt)
self.rdt += timedelta(hours = -GMT_hour_offset)
else:
raise Exception("bad datetime!")
self.tz = GMT_hour_offset
# uses the reference day of january 1st 2000
jan_1st_2000_jd = 2451545
jan_1st_2000 = datetime(2000,1,1,12,0,0)
time_del = self.rdt - jan_1st_2000
self.ajd = float(jan_1st_2000_jd) + float(time_del.total_seconds())/86400
self.ajc = (self.ajd - 2451545)/36525.0
return
[docs] def get_geomean_long(self):
""" :return geomean_long: geometric mean longitude of the sun"""
if not self.geomean_long is None:
return self.geomean_long
self.geomean_long = (280.46646 + self.ajc * (36000.76983 + self.ajc*0.0003032)) % 360
return self.geomean_long
[docs] def get_geomean_anom(self):
"""calculates the geometric mean anomoly of the sun"""
if not self.geomean_anom is None:
return self.geomean_anom
self.geomean_anom = (357.52911 + self.ajc * (35999.05029 - 0.0001537 * self.ajc))
return self.geomean_anom
[docs] def get_earth_eccent(self):
""" :return earth_eccent: precise eccentricity of earths orbit at referece datetime """
if not self.earth_eccent is None:
return self.earth_eccent
self.earth_eccent = 0.016708634 - self.ajc * (4.2037e-5 + 1.267e-7 * self.ajc)
return self.earth_eccent
[docs] def get_sun_eq_of_center(self):
""" :return sun_eq_of_center: the suns equation of center"""
if not self.sun_eq_of_center is None:
return self.sun_eq_of_center
if self.geomean_anom is None:
self.get_geomean_anom()
ajc = self.ajc
gma = radians(self.geomean_anom)
self.sun_eq_of_center = sin(gma) * (1.914602 - ajc*(0.004817 + 0.000014 * ajc)) + \
sin(2*gma) * (0.019993 - 0.000101 * ajc) + \
sin(3*gma) * 0.000289
return self.sun_eq_of_center
[docs] def get_true_long(self):
""" :return true_long: the tru longitude of the sun"""
if not self.true_long is None:
return self.true_long
if self.geomean_long is None:
self.get_geomean_long()
if self.sun_eq_of_center is None:
self.get_sun_eq_of_center()
self.true_long = self.geomean_long + self.sun_eq_of_center
return self.true_long
[docs] def get_true_anom(self):
""" :return true_anom: calculates the true anomaly of the sun"""
if not self.true_anom is None:
return self.true_anom
if self.geomean_long is None:
self.get_geomean_long()
if self.sun_eq_of_center is None:
self.get_sun_eq_of_center()
self.true_anom = self.geomean_anom + self.sun_eq_of_center
return self.true_anom
[docs] def get_rad_vector(self):
""" :return rad_vector: incident radiation vector to surface at ref_datetime (AUs)"""
if not self.rad_vector is None:
return self.rad_vector
if self.earth_eccent is None:
self.get_earth_eccent()
if self.true_anom is None:
self.get_true_anom()
ec = self.earth_eccent
ta = radians(self.true_anom)
self.rad_vector = (1.000001018*(1 - ec**2)) / (1 + ec *cos(ta))
return self.rad_vector
[docs] def get_app_long(self):
""" :return app_long: calculates apparent longitude of the sun"""
if not self.app_long is None:
return self.app_long
if self.true_long is None:
self.get_true_long()
stl = self.true_long
ajc = self.ajc
self.app_long = stl - 0.00569 - 0.00478 * sin(radians(125.04 - 1934.136 * ajc))
return self.app_long
[docs] def get_oblique_mean_elip(self):
""" :return oblique_mean_elip: oblique mean elliptic of earth orbit """
if not self.oblique_mean_elip is None:
return self.oblique_mean_elip
ajc = self.ajc
self.oblique_mean_elip = 23 + (26 + (21.448 - ajc * (46.815 + ajc * (0.00059 - ajc * 0.001813)))/60)/60
return self.oblique_mean_elip
[docs] def get_oblique_corr(self):
""" :return oblique_corr: the oblique correction """
if not self.oblique_corr is None:
return self.oblique_corr
if self.oblique_mean_elip is None:
self.get_oblique_mean_elip()
ome = self.oblique_mean_elip
ajc = self.ajc
self.oblique_corr = ome + 0.00256 * cos(radians(125.04 - 1934.136 * ajc))
return self.oblique_corr
[docs] def get_right_ascension(self):
""" :return right_ascension: the suns right ascension angle """
if not self.right_ascension is None:
return self.right_ascension
if self.app_long is None:
self.get_app_long()
if self.oblique_corr is None:
self.get_oblique_corr()
sal = radians(self.app_long)
oc = radians(self.oblique_corr)
self.right_ascension = degrees(arctan2(cos(oc) * sin(sal), cos(sal)))
return self.right_ascension
[docs] def get_declination(self):
""" :return declination: solar declination angle at ref_datetime"""
if not self.declination is None:
return self.declination
if self.app_long is None:
self.get_app_long()
if self.oblique_corr is None:
self.get_oblique_corr()
sal = radians(self.app_long)
oc = radians(self.oblique_corr)
self.declination = degrees(arcsin((sin(oc) * sin(sal))))
return self.declination
[docs] def get_equation_of_time(self):
""" :return equation_of_time: the equation of time in minutes """
if not self.equation_of_time is None:
return self.equation_of_time
if self.oblique_corr is None:
self.get_oblique_corr()
if self.geomean_long is None:
self.get_geomean_long()
if self.geomean_anom is None:
self.get_geomean_anom()
if self.earth_eccent is None:
self.get_earth_eccent()
oc = radians(self.oblique_corr)
gml = radians(self.geomean_long)
gma = radians(self.geomean_anom)
ec = self.earth_eccent
vary = tan(oc/2)**2
self.equation_of_time = 4 * degrees(vary * sin(2*gml) - 2 * ec * sin(gma) +
4 * ec * vary * sin(gma) * cos(2 * gml) -
0.5 * vary * vary * sin(4 * gml) -
1.25 * ec * ec * sin(2 * gma))
return self.equation_of_time
[docs] def get_hour_angle_sunrise(self):
""" :return hour_angle_sunrise: the hour angle of sunrise """
if not self.hour_angle_sunrise is None:
return self.hour_angle_sunrise
if self.declination is None:
self.get_declination()
d = radians(self.declination)
lat = self.lat_r
self.hour_angle_sunrise = degrees(arccos((cos(radians(90.833)) /
(cos(lat) * cos(d)) - tan(lat) * tan(d))))
return self.hour_angle_sunrise
[docs] def get_solar_noon(self):
""" :return solar_noon: solar noon in (local sidereal time LST)"""
if not self.solar_noon is None:
return self.solar_noon
if self.equation_of_time is None:
self.get_equation_of_time
lon = self.lon
eot = self.equation_of_time
tz = self.tz
self.solar_noon = (720 - 4 * lon - eot + tz * 60)/1440
# format this as a time for display purposes (Hours:Minutes:Seconds)
if self.is_numpy:
self.solar_noon_time = timedelta(days = self.solar_noon.mean())
else:
self.solar_noon_time = timedelta(days = self.solar_noon)
return self.solar_noon
[docs] def get_sunrise(self):
""" :return sunrise: returns the time of sunrise"""
if not self.sunrise is None:
return self.sunrise
if self.solar_noon is None:
self.get_solar_noon()
if self.hour_angle_sunrise is None:
self.get_hour_angle_sunrise()
sn = self.solar_noon
ha = self.hour_angle_sunrise
self.sunrise = (sn * 1440 - ha * 4)/1440
# format this as a time for display purposes (Hours:Minutes:Seconds)
if self.is_numpy:
self.sunrise_time = timedelta(days = self.sunrise.mean())
else:
self.sunrise_time = timedelta(days = self.sunrise)
return self.sunrise
[docs] def get_sunset(self):
""" :return sunset: returns the time of sunset"""
if not self.sunset is None:
return self.sunset
if self.solar_noon is None:
self.get_solar_noon()
if self.hour_angle_sunrise is None:
self.get_hour_angle_sunrise()
sn = self.solar_noon
ha = self.hour_angle_sunrise
self.sunset = (sn * 1440 + ha * 4)/1440
# format this as a time for display purposes (Hours:Minutes:Seconds)
if self.is_numpy:
self.sunset_time = timedelta(days = self.sunset.mean())
else:
self.sunset_time = timedelta(days = self.sunset)
return self.sunset
[docs] def get_sunlight(self):
""" :return sunlight: amount of daily sunlight in fractional days"""
if not self.sunlight is None:
return self.sunlight
if self.hour_angle_sunrise is None:
self.get_hour_angle_sunrise()
self.sunlight = 8 * self.hour_angle_sunrise / (60 * 24)
# format this as a time for display purposes (Hours:Minutes:Seconds)
if self.is_numpy:
self.sunlight_time = timedelta(days = self.sunlight.mean())
else:
self.sunlight_time = timedelta(days = self.sunlight)
return self.sunlight
[docs] def get_true_solar(self):
""" :return true_solar: true solar time at ref_datetime"""
if not self.true_solar is None:
return self.true_solar
if self.equation_of_time is None:
self.get_equation_of_time
lon = self.lon
eot = self.equation_of_time
# turn reference datetime into fractional days
frac_sec = (self.rdt - datetime(self.rdt.year, self.rdt.month, self.rdt.day)).total_seconds()
frac_hr = frac_sec / (60 * 60) + self.tz
frac_day = frac_hr / 24
self.frac_day = frac_day
# now get true solar time
self.true_solar = (frac_day * 1440 + eot + 4 * lon - 60 * self.tz) % 1440
# format this as a time for display purposes (Hours:Minutes:Seconds)
if self.is_numpy:
self.true_solar_time = timedelta(days = self.true_solar.mean() / (60*24))
else:
self.true_solar_time = timedelta(days = self.true_solar / (60*24))
return self.true_solar
[docs] def get_hour_angle(self):
""" :return hour_angle: returns hour angle at ref_datetime"""
if not self.hour_angle is None:
return self.hour_angle
if self.true_solar is None:
self.get_true_solar()
ts = self.true_solar
# matrix hour_angle calculations
if self.is_numpy:
ha = ts
ha[ha <= 0] = ha[ha <= 0]/4 + 180
ha[ha > 0] = ha[ha > 0]/4 - 180
self.hour_angle = ha
# scalar hour_angle calculations
else:
if ts <= 0:
self.hour_angle = ts/4 + 180
else:
self.hour_angle = ts/4 - 180
return self.hour_angle
[docs] def get_zenith(self):
""" :return zenith: returns solar zenith angle at ref_datetime"""
if not self.zenith is None:
return self.zenith
if self.declination is None:
self.get_declination()
if self.hour_angle is None:
self.get_hour_angle()
d = radians(self.declination)
ha = radians(self.hour_angle)
lat = self.lat_r
self.zenith = degrees(arccos(sin(lat) * sin(d) + cos(lat) * cos(d) * cos(ha)))
return self.zenith
[docs] def get_elevation(self):
""" :return elevation: returns solar elevation angle at ref_datetime"""
if not self.elevation is None:
return self.elevation
if self.zenith is None:
self.get_zenith()
# perform an approximate atmospheric refraction correction
# matrix hour_angle calculations
# these equations are hideous, but im not sure how to improve them without
# adding computational complexity
if self.is_numpy:
e = 90 - self.zenith
ar = e * 0
ar[e > 85] = 0
ar[(e > 5) & (e <=85)] = 58.1 / tan(radians(e[(e > 5) & (e <=85)])) - \
0.07 / tan(radians(e[(e > 5) & (e <=85)]))**3 + \
0.000086 / tan(radians(e[(e > 5) & (e <=85)]))**5
ar[(e > -0.575) & (e <= 5)] = 1735 + e[(e > -0.575) & (e <= 5)] * \
(103.4 + e[(e > -0.575) & (e <= 5)] * (-12.79 + e[(e > -0.575) & (e <= 5)] * 0.711))
ar[e <= -0.575] = -20.772 / tan(radians(e[e <= -0.575]))
# scalar hour_angle calculations
else:
e = 90 - self.zenith
er = radians(e)
if e > 85: ar = 0
elif e > 5: ar = 58.1 / tan(er) - 0.07 / tan(er)**3 + 0.000086 / tan(er)**5
elif e > -0.575: ar = 1735 + e * (103.4 + e * ( -12.79 + e * 0.711))
else: ar = -20.772 / tan(er)
self.elevation_noatmo = e
self.atmo_refraction = ar / 3600
self.elevation = self.elevation_noatmo + self.atmo_refraction
return self.elevation
[docs] def get_azimuth(self):
""" :return azimuth: returns solar azimuth angle at ref_datetime"""
if not self.azimuth is None:
return self.azimuth
if self.declination is None:
self.get_declination()
if self.hour_angle is None:
self.get_hour_angle()
if self.zenith is None:
self.get_zenith()
lat = self.lat_r
d = radians(self.declination)
ha = radians(self.hour_angle)
z = radians(self.zenith)
# matrix azimuth calculations
# these equations are hideous monsters, but im not sure how to improve them without
# adding computational complexity.
if self.is_numpy:
az = ha * 0
az[ha > 0] = (degrees(arccos(((sin(lat[ha > 0]) * cos(z[ha > 0])) - sin(d)) / (cos(lat[ha > 0]) * sin(z[ha > 0])))) + 180) % 360
az[ha <=0] = (540 - degrees(arccos(((sin(lat[ha <=0]) * cos(z[ha <=0])) -sin(d))/ (cos(lat[ha <=0]) * sin(z[ha <=0]))))) % 360
self.azimuth = az
else:
if ha > 0:
self.azimuth = (degrees(arccos(((sin(lat) * cos(z)) - sin(d)) / (cos(lat) * sin(z)))) + 180) % 360
else:
self.azimuth = (540 - degrees(arccos(((sin(lat) * cos(z)) -sin(d))/ (cos(lat) * sin(z))))) % 360
return self.azimuth
[docs] def get_earth_distance(self):
"""
:return earth_distance: distance between the earth and the sun at ref_datetime
"""
if self.rad_vector is None:
self.get_rad_vector()
# convert rad_vector length from AU to meters
self.earth_distance = self.rad_vector * 149597870700
return self.earth_distance
[docs] def get_norm_irradiance(self):
"""
Calculates incoming solar energy in W/m^2 to a surface normal to the sun.
inst_irradiance is calculated as = sun_surf_radiance\*(sun_radius / earth_distance)^2
and is then corrected as a function of solar incidence angle
:return norm_irradiance: the normal irradiance in W/m^2
"""
if not self.norm_irradiance is None:
return self.norm_irradiance
if self.earth_distance is None:
self.get_earth_distance()
ed = self.earth_distance
# calculate irradiance to normal surface at earth distance
self.norm_irradiance = self.sun_surf_rad * (self.sun_radius / ed)**2
return self.norm_irradiance
[docs] def get_inc_irradiance(self):
"""
calculates the actual incident solar irradiance at a given lat,lon coordinate
with adjustments for slope and aspect if they have been given. Not finished.
"""
print("this function is unfinished!")
return
[docs] def summarize(self):
""" prints attribute list and corresponding values"""
for key in sorted(self.__dict__.iterkeys()):
print("{0} {1}".format(key.ljust(20),sc.__dict__[key]))
return
[docs] def compute_all(self):
"""
Computes and prints all the attributes of this solar object. Spatial
averages are printed for numpy array type attributes.
"""
print("="*50)
print("Interogation of entire matrix of points.")
print("Some values displayed below are spatial averages")
print("="*50)
if self.is_numpy: # print means of lat/lon arrays
print("latitude, longitude \t{0}, {1}".format(self.lat.mean(), self.lon.mean()))
else:
print("latitude, longitude \t{0}, {1}".format(self.lat, self.lon))
print("datetime \t\t{0} (GMT)".format(self.rdt))
print("time zone \t\t{0} (GMT offset)".format(self.tz))
print("")
print("abs julian day \t\t{0}\t (day)".format(self.ajd))
print("abs julian century \t{0}\t (cen)".format(self.ajc))
print("suns geomean long \t{0}\t (deg)".format(self.get_geomean_long()))
print("suns geomean anom \t{0}\t (deg)".format(self.get_geomean_anom()))
print("earth eccentricity \t{0}".format(self.get_earth_eccent()))
print("suns eq of center \t{0}".format(self.get_sun_eq_of_center()))
print("suns true long \t\t{0}\t (deg)".format(self.get_true_long()))
print("suns true anom \t\t{0}\t (deg)".format(self.get_true_anom()))
print("suns apparent long \t{0}\t (deg)".format(self.get_app_long()))
print("earth obliq mean elip \t{0}\t (deg)".format(self.get_oblique_mean_elip()))
print("earth obliq correction\t{0}\t (deg)".format(self.get_oblique_corr()))
print("sun right ascension \t{0}\t (deg)".format(self.get_right_ascension()))
print("solar declination angle {0}\t (deg)".format(self.get_declination()))
print("equation of time \t{0}\t (min)".format(self.get_equation_of_time))
if self.is_numpy: # print means of hour angle array
print("hour angle sunrise\t{0}\t (deg)".format(self.get_hour_angle_sunrise().mean()))
else:
print("hour angle sunrise\t{0}\t (deg)".format(self.get_hour_angle_sunrise()))
print("")
self.get_solar_noon()
print("solar noon \t\t{0}\t (HMS - LST)".format(self.solar_noon_time))
self.get_sunrise()
print("sunrise \t\t{0}\t (HMS - LST)".format(self.sunrise_time))
self.get_sunset()
print("sunset \t\t{0}\t (HMS - LST)".format(self.sunset_time))
self.get_sunlight()
print("sunlight durration \t{0}\t (HMS)".format(self.sunlight_time))
self.get_true_solar()
print("true solar time \t{0}\t (HMS - LST)".format(self.true_solar_time))
print("")
if self.is_numpy: # print means of these array objects
print("hour angle \t\t{0}\t (deg)".format(self.get_hour_angle().mean()))
print("solar zenith angle \t{0}\t (deg)".format(self.get_zenith().mean()))
print("solar elevation angle \t{0}\t (deg)".format(self.get_elevation().mean()))
print("solar azimuth angle \t{0}\t (deg)".format(self.get_azimuth().mean()))
else:
print("hour angle \t\t{0}\t (deg)".format(self.get_hour_angle()))
print("solar zenith angle \t{0}\t (deg)".format(self.get_zenith()))
print("solar elevation angle \t{0}\t (deg)".format(self.get_elevation()))
print("solar azimuth angle \t{0}\t (deg)".format(self.get_azimuth()))
print("")
print("radiation vector \t{0}\t (AU)".format(self.get_rad_vector()))
print("earth sun distance \t{0}(m)".format(self.get_earth_distance()))
print("norm irradiance \t{0}\t (W/m*m)".format(self.get_norm_irradiance()))
print("="*50)
# testing
if __name__ == "__main__":
# use the current time and my time zone
my_datestamp = "20150515-120000" # date stamp
my_fmt = "%Y%m%d-%H%M%S" # datestamp format
my_tz = -4 # timezone (GMT/UTC) offset
my_lat = 37 # lat (N positive)
my_lon = -76.4 # lon (E positive)
sc = solar(my_lat, my_lon, my_datestamp, my_tz, my_fmt)
sc.compute_all()
sc.summarize()
## # numpy array test
## import numpy
## my_lat = numpy.array([[36, 36],[38,38]])
## my_lon = numpy.array([[-77.4,-75.4],[-77.4,-75.4]])
## sm = solar(my_lat, my_lon, my_datestamp, my_tz, my_fmt)
##
## sm.compute_all()