File:Longitude of perihelion of earth 140000 0 bp 1.svg

1. 2. Plot ice age orbital parameters
   3. with Python climlab
   5. 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

1. years = np.linspace(-5000, 0, 5001)

yearsmax1=-140 yearsmax2=(-1*yearsmax1)+1

fig, ax1 = plt.subplots()

1. fig, ax1=plt.gca()

1. ax2 = ax1.twinx()

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

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

print(yearsmax1)

1. quit(-1)

orb = OrbitalTable.interp(kyear=years)

print(orb)

1. quit(-1)

kyear = orb['kyear'] ecc = orb['ecc'] long_peri = orb['long_peri']%360 obliquity = orb['obliquity'] precession = orb['precession'] insolation_65n=orb['65NJul']

1. print (obliquity)

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)

1. print(S65)

1. ax1.set_xticks(years, yearlabels)

1. ax1.set_xticks(yearlabels2)

1. 1. plt.suptitle("Insolation at 65 N W/m2",size=18 )
2. plt.suptitle("Obliquity of Earth axis, degrees",size=18 )
3. plt.suptitle("Precession of Earth axis, degrees",size=18 )

plt.suptitle("Longitude of Earth perihelion", size=18)

ax1.set_xlabel("Thousands of years ago ",size=20 )

1. ax1.set_ylabel("Obliquity of Earth axis °",color="#7f0000", size=20)
2. ax1.set_ylabel("Precession of Earth axis",color="#7f0000", size=20)

ax1.set_ylabel("Longitude of Earth perihelion",color="#7f0000", size=20)

1. ax1.plot(years, obliquity)

ax1.plot(years, long_peri)

1. ax1.plot(years, ecc)
2. ax1.plot(years, precession)
3. ax1.plot(years, insolation_65n)

1. ax1.plot(years2, S65_15_smoothed, color="#0000ff", label="day 15 (winter)")
2. ax1.plot(years2, S65_180_smoothed, color="#ff0000", label="day 180 (summer)")
3. ax2.plot(years2, S65_90_smoothed, color="#00ff00", linestyle="--", label="day 90 (spring)")

1. ax1.set_ylabel("Insolation W/m2",color="#7f0000", size=20)
2. ax2.set_ylabel("Insolation W/m2",color="#007f00", size=20)

ax1.xaxis.set_tick_params(labelsize=16) ax1.yaxis.set_tick_params(labelsize=16)

1. ax2.yaxis.set_tick_params(labelsize=16)

ax1.set_xticks(xlabels) ax1.set_xticklabels(xlabels2)

1. ax1.set_xticklabels(["$0$"])
2. ax1.set_xticklabels(yearlabels)

1. ax1.axis["right"].major_ticklabels.set_visible(False)
2. ax1.axis["top"].major_ticklabels.set_visible(True)

ax1.legend(fontsize=16)

1. ax2.legend(fontsize=16)

plt.show()