%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() 

#Now we can draw the curves of Milankovitch cycles Firstyear=150 plot_step_function([(i,Obl[i]) for i in range(Firstyear)]) 

plot_step_function([(i,Ecc[i]) for i in range(Firstyear)]) 

plot_step_function([(i,Prec[i]) for i in range(Firstyear)]) 

#Now we compute the Imbrie and Imbrie forcing function. #First we set the constants according to their choice. Tw=10600 Tc=42500 alpha=-2 phi=2*pi*2000/22000 #Then we define the function denoted by "x" in Imbrie and Imbrie. #Since we only know the celestial mechanics functions as a table, we #can only define x for the values t=- 1000 * N where N is a natural number #less or equal to 51000. We define the function II(N) #(where II(N) = x(-1000*N)). def II(n): return(Obl[n]+alpha*Ecc[n]*sin(Prec[n]-phi)) # To see what we got, we plot the function. # This is supposed to agree with the equilibrium response plot_step_function([(i,II(i)) for i in range(Firstyear)]) 

#Now we compute the solution to the Imbrie and Imbrie differential equation. #Actually as a difference equation with steps=1000 years. #Because of the recursive nature of the solution, we will write it down as an array. #We initialise with the value 0.44, though this won't be important in the long run. InitialII=0.44 #Here is the step in the difference equation def IIStep(n): a=IISol[n]-II(n) if a > 0: return(IISol[n]-a*1000/Tc) else: return(IISol[n]-a*1000/Tw) #Next we define the array corresponding to IISol, and populate it. IISol = numpy.array(range(Era[1]+1),dtype=float) IISol[Era[1]]=InitialII for i in range(Era[1]): IISol[Era[1]-i-1]=IIStep(Era[1]-i) #Finally we take a look at the solution (in blue) #together with the equilibrium response (in red) graph1=plot_step_function([(i,IISol[i]) for i in range(Firstyear)],rgbcolor="red") graph2=plot_step_function([(i,II(i)) for i in range(Firstyear)]) graph1+graph2