# This script computes and plots seasonally varying
# and annual mean solar flux for the case of a circular
# orbit

from solar import *
from math import *
from ClimateUtilities import *

degrad = pi/180.
#Set the obliquity you want (in radians)
obliquity = degrad*90.
#
N = 90
phimax = pi/2.
dphi = 2.*phimax/(N-1)
phiList = [(dphi*i-phimax) for i in range(N)]
seasonList = [2*i*degrad for i in range(181)]

# The following statement shows how we can use list
# comprehension to make a list of a function evaluated
# at a list of arguments.
meanfluxList = [AnnMeanFluxCirc(phi,obliquity) for phi in phiList]
#
#Now plot the mean flux
c = Curve()
c.addCurve(phiList,'lat','Latitude')
c.addCurve(meanfluxList,'f','Annual Mean Flux')
c.switchXY = True
# Convert latitude to degrees for plotting.
c['lat'] = c['lat']/degrad
meanFluxPlot = plot(c)

# The following loop shows how a 2D Numpy array can be generated
# from lists of arguments, without needing to keep track of
# index counters.
fluxArray = [ ]
for phi in phiList:
    fluxArray.append([solar(phi,season,obliquity) for season in seasonList])

flux = numpy.array(fluxArray)

#Now make a contour plot of the flux
##seasonalPlot = contour(flux,x=[s/degrad for s in seasonList],
##                       y=[phi/degrad for phi in phiList])

#A little special expert code to make the tickmarks come out nicely
#Note that this explicitly uses PyNgl. If you are using a different
#graphic driver, this will be ignored.
r = Dummy()
r.tmXBMode = "Manual"
r.tmXBTickSpacingF = 60.
#
#Eliminate "colors" keyword argument to get the default rainbow colors
#Set colors to "grey" to get greyscale
#Set colors to "ncview_default" to get a nice looking rainbow scale
seasonalPlot = contour(flux,x=[s/degrad for s in seasonList],
                       y=[phi/degrad for phi in phiList],
                        resource = r,colors= 'ncview_default')