#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:  Separate out the basic flux computation functions, and
#       make them into a separate module that can be imported
#       into various scripts. This version is a stand-alone version,
#       which incorporates the old solar_flux.py module.

Chapter = "8"
Figure = None

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

from ClimateUtilities import *
from math import * #Because some of the routines need sin, not math.sin, etc.
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

#----------------------------------------------------------------
#Basic functions to get latitude of the Sun, zenith angle, and
#so forth. Input angles are all in radians

#Latitude of the Sun
def sunlat( season, gamma):
	return asin( cos(season)*sin(gamma))
#Hour angle giving the terminator position
def h( phi, delta):
	temp = -sin(phi)*sin(delta)/(cos(phi)*cos(delta));
	if temp < -1. :
		return pi
	elif temp > 1.:
		 return 0.
	else:
		 return acos(temp)

#Flux factor as a function of latitude and lat of the Sun
def flux( phi, delta):
	return (2.*sin(h(phi,delta))*cos(phi)*cos(delta) +
	     2.*h(phi,delta)*sin(phi)*sin(delta))/(2.*3.14159)

#  A convenience function, with a calling sequence simpler to use
#  but otherwise, it does the same as the function flux(phi,delta)
def solar(phi,season,obliquity):
	return flux(phi,sunlat(season,obliquity))

	
# The following function computes the annual
# mean flux, but it is only valid for a circular orbit
#
def AnnMeanFluxCirc( phi, obliquity):
        avg = 0.
        n = 180 # n should be even for symmetry
        ds = 2.*pi/n
	for season in [ds/2. + i*ds for i in range(n)]:
                avg += flux(phi,sunlat(season,obliquity))
	return avg/180.

#----------------------------------------------------------------

    
#----------------------------------------------------------------
#
#Functions for computing fluxes for an elliptical orbit
#
#Angular velocity function. e is the eccentricity,
#J is the orbital angular momentum and a is the
#semi-major axis
def kappadot(t,kappa,params):
    e = params.e
    J = params.J
    a = params.a
    return (J/a**2)*(1+e*math.cos(kappa))**2/(1-e**2)**2

#Distance as function of semi-major axis and angle
def r(a,e,kappa):
    return a*(1.-e*e)/(1.+e*math.cos(kappa))

#Function to do the flux computation
def fluxmap(e,obl,precess):
    #Set up the parameter object
    params = Dummy()
    params.e = e
    params.J = 1.
    params.a = 1.
    #Now do the integration

    kappaStart = precess #Start at the solstice, to make the graphs nicer
    #Set up the integrator and install parameter object
    #dt = 2.*math.pi/6000
    dt = 2.*math.pi/600
    m = integrator(kappadot,0.,kappaStart,dt)
    m.setParams(params)

    tlist = []
    fluxes = [[] for lat in lats] #Sets up a list of empty lists of len(lats)
    means = [0. for lat in lats]  #List of mean values
    meanr = 0.
    navg = 0
    kappa = kappaStart
    count = 0
    while kappa < 2.*pi +kappaStart:
        t,kappa = m.next()
        L = 1./r(params.a,params.e,kappa)**2
        tlist.append(t)
        navg += 1
        meanr += r(params.a,params.e,kappa)
        for j in range(len(lats)):
            if count%10 == 0:
              fluxes[j].append(L*flux(degrad*lats[j],sunlat(kappa-precess,obl)))
            means[j] += fluxes[j][-1] # Add the preceding value to the mean
        count += 1

    #Save Nondimensional length of year
    year = tlist[-1]
    # Normalize time so length of year= 1
    tlist = [tt/year for tt in tlist]
    # Compute the means
    means =  [sum/navg for sum in means]
    meanr = meanr/navg
    #Stuff results in a curve, and save
    c = Curve()
    c.addCurve(tlist,'Time')
    for i in range(len(lats)):
        c.addCurve(fluxes[i],'%f'%lats[i])
    return Numeric.array(fluxes)



#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
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)] #A Global variable
#Compute the latitude-time map of the flux
f = Lsun*fluxmap(e,obl*degrad,prec*degrad) #Convert to radians!
#Plot it
contour(f,y=lats)
#Now let's try a larger eccentricity
f = Lsun*fluxmap(.2,obl*degrad,prec*degrad)
contour(f,y=lats)