#Makes plots of flux factor for a planet in an elliptical orbit
#This is a version of ellipseFlux.py, designed to aid interactive
#explorations of the effect of orbital parameter variations on
#the distribution of incident solar radiation.
#
#ToDo:  It's a little confusing that the latitude list is
#       expected in degrees, and everything else is radians.

Chapter = "7"
Figure = None

#Edit to give the location of the data on your system
#datapath = '/home/rtp1/GeoSci232/WorkbookDatasets/'
datapath ='/Users/rtp1/Havsornen/PlanetaryClimateBook/WorkbookDatasets/'
dataset = 'Chapter7Data/Milankovic/orbit91'

from ClimateUtilities import *
from solar import * #Solar flux computation library
import math 

#Get the orbital parameters data
orbitParams = readTable(datapath+dataset)
degrad = math.pi/180.

#Function to look up and return Earth's orbital parameters
#The argument is time before present, in thousands of years.
#The function returns a triple consisting of eccentricity,
#obliquity (degrees) and the angle between the position of
#the Northern Hemisphere Summer Solstice and the position of
#the perihelion.
#
def getOrbit(time):
    e = orbitParams['ECC'][time]
    obl = orbitParams['OBL'][time]
    prec = orbitParams['OMEGA'][time] + 90. #relative to NH solstice
    if prec > 360.:
        prec -= 360.
    return e,obl,prec

	


    

#Function to make a map of the flux vs latitude and time
#Input angles are in degrees
def fluxmap(e,obl,precess,lats):
    fluxes = []
    orb = orbit(e,precess*degrad)#Start at the solstice, to make the graphs nicer
    #The next statement produces a 2D array each row of which is the
    #flux vs. time for a given latitude. Recall that the orbitFlux()
    #function returns a 1D array.
    fluxes = [orbitFlux(lat*degrad,obl*degrad,precess*degrad,orb) for lat in lats]
    #Stuff results in a curve, and save
    #c = Curve()
    #c.addCurve(orb['t'],'Time')
    #for i in range(len(lats)):
    #    c.addCurve(fluxes[i],'%f'%lats[i])
    return numpy.array(fluxes)

#Function to return a Curve() giving the annual mean
#flux vs latitude. Input angles in degrees
def mean(e,obl,precess,lats):
    c = Curve()
    c.addCurve(lats,'lat','Latitude')
    orb = orbit(e)
    means = [AnnMeanFlux(lat*degrad,obl*degrad,precess*degrad,orb) for lat in lats]
    c.addCurve(means,'Ann','Annual Mean Flux')
    return(c)

#Function to create a plot of the seasonal cycle
#of insolation at a given latitude
def fluxSeas(e,obl,prec,lat):
    flux = Lsun*fluxmap(e,obl,prec,[lat])[0]
    t = numpy.array(range(len(flux)))/float(len(flux))
    c = Curve()
    c.addCurve(t,'t')
    c.addCurve(flux,'flux')
    c.Xlabel = 'Time of Year'
    return c
#Like fluxSeas, but does it for a list of precessions
#Added for Abisko lectures, but generally useful
def fluxSeas1(e,obl,precList,lat):
    c = Curve()
    start = True
    for prec in precList:
     flux = Lsun*fluxmap(e,obl,prec,[lat])[0]
     if start:
         t = numpy.array(range(len(flux)))/float(len(flux))
         c.addCurve(t,'t')
         c.Xlabel = 'Time of Year'
         start = False
     c.addCurve(flux,'flux%f'%prec)
    return c

#Examples of use:

#(a) Compute the incident solar radiation at 60N for present
#    Earth obliquity, at the NH summer solstice, the equinox,
#    and the time halfway between. Also compute the annual mean flux for a
#    CIRCULAR orbit
#**ToDo:  Use degrees and radians consistently in argument lists
#         These functions take radians, but some others take degrees
#          Convention: All low level functions like solar take radians
obl = 23.4*degrad
lat = 60.*degrad
Lsun = 1370.
S60_solstice = Lsun*solar(lat,0.,obl)
S60_mid = Lsun*solar(lat,45.*degrad,obl)
S60_equinox = Lsun*solar(lat,90.*degrad,obl)
S60_ann = Lsun*AnnMeanFluxCirc(lat,obl)
print S60_solstice,S60_mid,S60_equinox,S60_ann

#Now let's do a case for a non-circular orbit
#
#First get the orbital parameters for 5000 years ago
e,obl,prec = getOrbit(5)
#Define the set of latitudes for which we want to make the plot
lats = [-90. + 2.*i for i in range(91)]
#Compute the latitude-time map of the flux
f = Lsun*fluxmap(e,obl,prec,lats) #Input is in degrees
#Plot it
contour(f,y=lats)
#Now let's try a larger eccentricity
f = Lsun*fluxmap(.2,obl,prec,lats) #Input is in degrees
contour(f,y=lats)