In [1]:

    

%run basics
import math
import pytz
%matplotlib


    

    




Using matplotlib backend: Qt4Agg

In [9]:

    

def GetRadiationDirect(utc_datetime, altitude_deg):
    # from Masters, p. 412
    if(altitude_deg > 0):
            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)
            return flux * math.exp(-1 * optical_depth * air_mass_ratio)*math.sin(math.radians(altitude_deg))
    else:
            return 0.0

def GetDayOfYear(utc_datetime):
        year_start = datetime.datetime(utc_datetime.year, 1, 1, tzinfo=utc_datetime.tzinfo)
        delta = (utc_datetime - year_start)
        return delta.days

def GetApparentExtraterrestrialFlux(day):
    # from Masters, p. 412
    return 1160 + (75 * numpy.sin(numpy.radians((360./365) * (day - 275))))


    

In [3]:

    

# read the netCDF file
nc_name = qcio.get_filename_dialog()
ds = qcio.nc_read_series(nc_name)


    

In [4]:

    

# Ian's code
# set a value for "k", the extinction coefficient
k = 0.9
# get the latitude and longitude
lat_decdeg = float(ds.globalattributes["latitude"])
long_decdeg = float(ds.globalattributes["longitude"])
if "elevation" in ds.globalattributes:
    ALT_m = int(ds.globalattributes["elevation"])
else:
    ALT_m = 1.0
rec_length = int(ds.globalattributes["time_step"])
# Ian's original code queried AskGeo for the UTC offset in milliseconds
# and converted that to hours.  Here we use the time zone given in the
# netCDF file global attributes.
# get the time zone
time_zone = ds.globalattributes["time_zone"]
# get a pytz object
tz = pytz.timezone(time_zone)
# use the same "current" date as Ian
cdt = datetime.datetime(2013,6,1,0,0,0)
# get the offset to UTC time as a timedelta object
utcoffset = tz.utcoffset(cdt)
# get the offset as hours
GMT_zone = utcoffset.seconds/float(3600)
# get the datetime
ldt = ds.series["DateTime"]["Data"]
# get the day of the year from local time
#DOY = numpy.array([dt.timetuple().tm_yday for dt in ldt])
start_doy = ldt[0].timetuple().tm_yday
end_doy = ldt[-2].timetuple().tm_yday
DOY = numpy.arange(start_doy,end_doy+1,1)
# For each day calculate equation of time correction, solar noon, declination and TOA radiation
array_EqofTime=0.17*numpy.sin(4*numpy.pi*(DOY-80)/373)-0.129*numpy.sin(2*numpy.pi*(DOY-8)/355) # DiLaura (1984)
array_solar_noon=12+(GMT_zone*15.0-long_decdeg)/360*24-array_EqofTime # Me
array_decl=numpy.radians(23.4)*numpy.sin((DOY+284)/365.0*2*numpy.pi) # Oke (1987)
array_TOArad=(1+0.034*numpy.cos(DOY/365.25*2*numpy.pi))*1367.0 # Duffie and Beckman (1980)


    

In [7]:

    

flux = GetApparentExtraterrestrialFlux(DOY)


    

In [8]:

    

fig = plt.figure()
plt.plot(DOY,array_TOArad,'b.')
plt.plot(DOY,flux,'r+')
plt.show()


    

In [ ]: