%hide %latex This is an current work that we are doing at the Azimuth project. If you want to contribute contact us. 

FILE= 'INSOLN.LA2004.BTL.ASC' import numpy #from numpy cimport ndarray INSOLN = open(DATA+FILE, 'r') # This file is available for download from "http://www.imcce.fr/Equipes/ASD/insola/earth/La2004/index.html" # It has data on the orbit of the Earth from -51 million years to present. Era=[0,51000] # The last and first Kyear in the file. Era[a,b] corresponds to # the period from the year -1000*b to the year -1000*a. # We are counting backwards here, so b is positive and greater than a. #We initialize four arrays, for the moment initialized with meaningless numbers. Year = numpy.array(range(Era[1]+1)) #The year Ecc = numpy.array(range(Era[1]+1),dtype=float) #Eccentricity Obl = numpy.array(range(Era[1]+1),dtype=float) #Obliquity in radians Prec = numpy.array(range(Era[1]+1),dtype=float) #longitude of perihelion from moving equino x in radians #I assume that the equinox is the vernal equinox (Check this?). #We populate the four arrays with data we read from the file INSOLN. for i in range(Era[1]+1): x=INSOLN.readline() y=[float(a.replace("D","E")) for a in x.split()] Year[i]= 1000*y[0] Ecc[i]= y[1] Obl[i]= y[2] Prec[i]= y[3] #All data have been read into memory now, so to please the computer, we close the file INSOLN.close() 

def Insolation(Millenium,Latitude,OrbitLongitude): s=n(max(0,1-(sin(Latitude))^2-sin(OrbitLongitude)^2*sin(Obl[Millenium])^2)) p=n(sin(Latitude)*sin(Obl[Millenium])*sin(OrbitLongitude)) e=Ecc[Millenium] Omega=Prec[Millenium] x=s+p^2 if x<=0: return(0) else: return(n((1-e*cos(OrbitLongitude-Omega)/(1-e^2))^2*(p*arccos(-p/sqrt(x))+sqrt(s)))) # Latitude is given in radians, the North Pole is pi, # the South Pole is -pi. The input should lie between these two numbers, # but for now, we don't check if the input is valid. # # Orbitlongitude is given in radians. For instance, summer solstice is pi/2. # It is possible to convert this to a date of the year, but the conversion is not # entirely trivial. # # The formula is taken from Paillard "Climate and the orbital parameters of the Earth" # It's not valid close to the poles at the times when the sun does not rise or set, I'll look # into that later. 

# Now we want to draw curves over insolation, which should reproduce the curve # on to of figure 2 in Paillard and Parenin FirstYear=500 SummerDate = pi/2 NorthLat=65*pi/180 WinterDate = -pi/6 SouthLat=-60*pi/180 # The values pi/2 and -pi/6 are not exact, but not so far from June 21 and February 21. # For the purpose of checking that we are in the right ballpark, it should do fine. # We are outside the polar circle at all times, so the formula above # is correct. NorthInso = numpy.array(range(FirstYear+1),dtype=float) SouthInso = numpy.array(range(FirstYear+1),dtype=float) for i in range(FirstYear+1): NorthInso[i]=Insolation(i,NorthLat,SummerDate) SouthInso[i]=Insolation(i,SouthLat,WinterDate) 

# We reconstruct figure 2. The insolation is lower at -60 than at 65, since # the insolation on the North hemisphere is computed at summer solstice, # and the insolation at the South hemisphere is computed late in the Southern summer. # To get a nicer picture, we compensate fro this by a linear transformation. # graph1=plot_step_function([(i,0.7*NorthInso[i]) for i in range(FirstYear)],rgbcolor="red") graph2=plot_step_function([(i,SouthInso[i]) for i in range(FirstYear)]) 

graph1 + graph2 

# The following are parametervalues picked by Palliard and Parenin ppsV = 15 ppsC = 5 ppsA = 12 ppx = 1.3 ppy = 0.5 ppz = 0.8 ppalpha = 0.15 ppbeta = 0.5 ppgamma = 0.5 ppdelta = 0.4 ppa = 0.3 ppb = 0.7 ppc = 0.01 ppd = 0.27 

# The dynamical system def ppF(Millenium,ppV,ppA): return(ppa*ppV-ppb*ppA-ppc*SouthInso[Millenium]+ppd) def ppVr(Millenium,ppC): return(-ppx*ppC-ppy*NorthInso[Millenium]+ppz) def ppCr(Millenium,ppV,ppA): x=ppalpha*NorthInso[Millenium]-ppbeta*ppV+ppdelta if ppF(Millenium,ppV,ppA) >=0: return(x) else: return(x+ppgamma) # We start the dynamical system, and let it run # The initial values are chosen to be close to empirically noticed # long term averages. C = numpy.array(range(FirstYear+1),dtype=float) A = numpy.array(range(FirstYear+1),dtype=float) V = numpy.array(range(FirstYear+1),dtype=float) C[FirstYear]=1.295 A[FirstYear]=-1.455 V[FirstYear]=-1.46 for i in range(FirstYear): j=FirstYear-i-1 C[j]=C[j+1]+(ppCr(j,V[j+1],A[j+1])-C[j+1])/ppsC A[j]=A[j+1]+(V[j+1]-A[j+1])/ppsA V[j]=V[j+1]+(ppVr(j,C[j+1])-V[j+1])/ppsV 

# I don't know why I now graph -V instead of V, but it seems to work. # I'll be back on this. graph3=plot_step_function([(i,-V[i]) for i in range(FirstYear)]) graph3