File:Milankovic cycles 26-0 ka bp 1.svg

1. 3. Plot ice age orbital parameters
   5. with Python climlab
   7. many curves of parameters
   9. 19.10.2021 0000.0000

import numpy as np import matplotlib.pyplot as plt

from scipy.interpolate import make_interp_spline, BSpline

from climlab.solar.orbital import OrbitalTable from climlab.solar.insolation import daily_insolation from climlab import constants as const

yearsmax1=-26 numxlabels=14

yearsmax2=(-1*yearsmax1)+1

years = np.linspace(yearsmax1, 0,yearsmax2) yearlabels=np.linspace(yearsmax1,0)

xlabels=np.linspace(yearsmax1,0,numxlabels) xlabels2=xlabels*-1 xlabels2[numxlabels-1]=0

orb = OrbitalTable.interp(kyear=years)

1. print(orb)

kyear = orb['kyear'] ecc = orb['ecc'] long_peri = orb['long_peri']%360 obliquity = orb['obliquity'] precession = orb['precession'] insolation_65n=orb['65NJul'] sin_long_peri=np.sin(np.deg2rad(long_peri))

len1=len(kyear)

inso1=insolation_65n[len1-1].values insolation_65n_0=np.repeat(inso1,len1)

sinlo1=sin_long_peri[len1-1].values sin_long_peri_0=np.repeat(sinlo1,len1)

prece1=precession[len1-1].values precession_0=np.repeat(prece1,len1)

ecco1=ecc[len1-1].values ecc_0=np.repeat(ecco1,len1)

oblio1=obliquity[len1-1].values obliquity_0=np.repeat(oblio1,len1)

S65_15 = daily_insolation(lat=65, day=15, orb=orb) S65_90 = daily_insolation(lat=65, day=90, orb=orb) S65_180 = daily_insolation(lat=65, day=180, orb=orb)

years2 = np.linspace(years.min(), years.max(), 200) yearlabels2=np.linspace(yearsmax1,0,200)

1. define spline

spl = make_interp_spline(years, S65_15, k=3) S65_15_smoothed = spl(years2) spl = make_interp_spline(years, S65_90, k=3) S65_90_smoothed = spl(years2) spl = make_interp_spline(years, S65_180, k=3) S65_180_smoothed = spl(years2)

plt.figure(1)

plt.subplot(411)

1. plt.xticks(direction='in')

plt.xticks(ticks=xlabels, labels=[]) plt.yticks(fontsize=14) plt.title("Milanković cycles, last 26 000 years",color="#7f0000", size=18) plt.ylabel("Eccentricity", fontsize=15) plt.plot(years,ecc, color="#00007f" ) plt.plot(years,ecc_0, color="#00007f", linestyle=":" ) plt.subplot(412)

1. plt.xticks(fontsize=14)

plt.xticks(ticks=xlabels, labels=[]) plt.yticks(fontsize=14) plt.ylabel("Oliquity", fontsize=15) plt.plot(years,obliquity, color="#00007f") plt.plot(years,obliquity_0, color="#00007f", linestyle=":" ) plt.subplot(413)

1. plt.xticks(fontsize=14)

plt.xticks(ticks=xlabels, labels=[]) plt.yticks(fontsize=14) plt.ylabel("Precession", fontsize=15) plt.plot(years,precession, color="#00007f" ) plt.plot(years,precession_0, color="#00007f", linestyle=":" )

1. plt.subplot(514)
2. plt.xticks(fontsize=14)

1. plt.xticks(ticks=xlabels, labels=[])
2. plt.yticks(fontsize=14)
3. plt.ylabel("sin (Per.)", fontsize=15)
4. plt.plot(years,sin_long_peri, color="#00007f" )

plt.subplot(414)

1. plt.xticks(fontsize=14)

plt.xticks(ticks=xlabels,labels=xlabels2, fontsize=14) plt.yticks(fontsize=14) plt.ylabel("Insolation", fontsize=15) plt.plot(years,insolation_65n, color="#ff7f00" ) plt.plot(years,insolation_65n_0, color="#ff7f00", linestyle=":" )

plt.xlabel("Thousands of years ago ",size=16 )

plt.show()