%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))/pi)) # 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" # 

# Now we want to draw curves over insolation, which should reproduce the curve # on to of figure 2 in Paillard and Parenin FirstYear=3000000 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. # 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) # # These are preliminary values. For comparison purposes, we normalise both # time series to get average 0, standard deviation 1. # NorthMean=mean(NorthInso) NorthDev=std(NorthInso) SouthMean=mean(SouthInso) SouthDev=std(SouthInso) def NorthNorm(x): return((x-NorthMean)/NorthDev) def SouthNorm(x): return((x-SouthMean)/SouthDev) NorthInso=map(NorthNorm,NorthInso) SouthInso=map(SouthNorm,SouthInso) 

# We reconstruct the top of figure 2. # We concentrate on the last half million years. Window=[0,500] graph1=plot_step_function([(i,NorthInso[i]) for i in range(Window[1])],rgbcolor="red") graph2=plot_step_function([(i,SouthInso[i]) for i in range(Window[1])]) 

graph1 + graph2 

# The following are parameter values picked by Palliard and Parrenin ppsV = 15 ppsC = 5 ppsA = 12 ppx = 1.3 ppy = 0.5 ppz = 0.8 ppalpha = 0.15 ppbeta = 0.5 ppgamma = 0.7 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+8.5*Millenium/100000) 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 C = numpy.array(range(FirstYear+1),dtype=float) A = numpy.array(range(FirstYear+1),dtype=float) V = numpy.array(range(FirstYear+1),dtype=float) H = numpy.array(range(FirstYear+1),dtype=float) F = numpy.array(range(FirstYear+1),dtype=float) C[FirstYear]=0 A[FirstYear]=0 V[FirstYear]=0 for i in range(FirstYear): j=FirstYear-i-1 ppV=ppVr(j,C[j+1]) ppC=ppCr(j,V[j+1],A[j+1]) ppF1=ppF(j,V[j+1],A[j+1]) C[j]=C[j+1]+(ppC-C[j+1])/ppsC A[j]=A[j+1]+(V[j+1]-A[j+1])/ppsA V[j]=V[j+1]+(ppV-V[j+1])/ppsV F[j]=ppF1 if ppF1 >=0: H[j]=1 else: H[j]=0 

graph0=plot_step_function([(i,F[i]) for i in range(Window[1])]) 

# The following graph should agree with the bottom graph of figure 2. The last 300 years look fine, # the rest not so much. graph0 

graph1=plot_step_function([(i,-V[i]) for i in range(Window[1])]) graph2=plot_step_function([(i,-A[i]) for i in range(Window[1])],rgbcolor="red") 

# The following graph shopuld correspond to part of #figure 2" # The last 300,000 years look fine, before that there seems to be some disagreement. graph1+graph2 

graph2=plot_step_function([(i,C[i]) for i in range(Window[1])]) 

# Next graph should correspond to the black graph in the middle of figure 2. graph2 

a=animate([ circle((V[t],A[t]), 0.001, rgbcolor=(0,0,1),fill=true)for t in srange(100,0,-1)], xmin=-1.48,ymin=-1.48,xmax=-1.43,ymax=-1.43,figsize=[10,10]) 

# It seems that in this range, C[i] and V[i] are closely coupled. # graph4=plot_step_function([(C[i],V[i]) for i in range(500)],rgbcolor="green") graph4 

# V[i] and A[i] can vary in a more complex fashion, graph5=plot_step_function([(A[i],V[i]) for i in range(500)],rgbcolor="blue") graph5 

def fcol(F1): if F1>0: return((0,0,1)) else: return((1,0,0)) a1=animate([text(str(1000*t)+" BP ",(0.1,1.2))+ circle((A[t],V[t]), 0.03, rgbcolor=fcol(F[t]) ,fill=true)for t in srange(1100,900,-1)], xmin=0,ymin=0,xmax=1,ymax=1.3,figsize=[10,10]) 

# The following animation illustrates the dynamic around the mid-pleistocene transition. a1.show() 

a2=animate([text(str(1000*t)+" BP ",(0.1,1.2))+ circle((A[t],V[t]), 0.03, rgbcolor=fcol(F[t]) ,fill=true)for t in srange(200,0,-1)], xmin=0,ymin=0,xmax=1,ymax=1.3,figsize=[10,10]) 

# Here are the last 200,000 years, for comparison. a2.show()