Poniższy notebook ma na celu zilustrowanie oscylacji parametrów orbity Ziemi, stanowiących główne źródło zmian nasłonecznienia, będącego siłą napędową zmian klimatu w ,,średniej’’ skali czasowej dziesiątek tysiecy  i milionów lat. Opieram się o obliczeniach Laskar 2010, demonstrujących parametry orbity Ziemi przez ostatnie 250 mln lat, przy czym tylko 20% z tego (50 mln lat wstecz) uważa się za wiarygodne. Przez ostatnie 10 000 lat dane początkowe są skalibrowane obserwacjami astromonicznymi.

(* http://vo.imcce.fr/insola/earth/online/earth/earth.html *)

Zaczynamy od wczytania danych. Nieco cierpliwości, trwa to ok. 30 sekund.

{czas, raw}=Import["http://vo.imcce.fr/insola/earth/online/earth/La2010/La2010d_alkhqp3L.dat"]//AbsoluteTiming;

czas

Typowa struktura PLIKU danych zawiera klasyczne parametry orbity eliptycznej, o których była mowa na wcześniejszym wykładzie. Spróbujemy tę wiedzę wykorzystać w praktyce. Dane są zapisane w krokach co 1000 lat.

raw[[1]] (* The file La2010X_alkhqp3.dat and La2010X_alkhqp3L.dat contain t,a,l,k,h,q,p
where t is the time from J2000 (in kyr)
a:semi-major axis
l:mean longitude (expressed in radians)
k:e cos (varpi)
h:e sin (varpi)
q:sin(i/2) cos (Omega)
p:sin(i/2) sin (Omega) *)

a=raw[[1;;,{1,2}]]; (* Wielka półoś orbity *)

ListPlot[a[[1;;112]],Joined->True,InterpolationOrder->0]

Dane są de facto zapisane w postaci liczb zespolonych, w postaci trygonometrycznej. Upraszczając, k + i h odpowiada wektorowi Rungego-Lenza, natomiast q + i p momentowi pędu.

(*Extract time (t),k,and h columns*)
t=raw[[All,1]]; (* Czas w tysiącach lat, kyr *)
k=raw[[All,4]];
h=raw[[All,5]];
q=raw[[All,6]];
p=raw[[All,7]];

Obliczamy na początek precesję orbity.

varpi=ArcTan[k,h];
(*Calculate varpi (longitude of perihelion in radians)*)
(*ArcTan[x,y] computes the arc tangent of y/x,giving a result in (-Pi,Pi]*)
(*We need ArcTan[k,h] to get atan2(h,k) which corresponds to varpi*)

varpi2=Arg[k+I h]; (* Inny sposób *)

varpi-varpi2//MinMax

Rozbieżność na poziomie 2 ULP (0 ULP na Mac ARM) oznacza, że z punktu widzenia arytmetyki liczb zmiennoprzecinkowych wyniki są prawie identyczne.

%/$MachineEpsilon

e=Sqrt[k^2+h^2]; (* Mimośród orbity *)

i=2 ArcSin[Sqrt[q^2+p^2]]; (* NAchylenie płaszczyzny orbity *)

Precesja podana jest w zakresie kątów -π do π; po wykonaniu pełnego obrotu elipsy wynik przekakuje o 2 π.

ListPlot[varpi[[1;;3012]],PlotMarkers->{Automatic, Tiny}]

Naiwnie oczekujemy, że elipsa (jej półoś) będzie obracała się ze stałą prędkością kątową. Nie jet to generalnie prawda. Co gorsza, tempo precesji bywa ujemne....

ListPlot[Differences[varpi[[1;;312]]],PlotMarkers->{Automatic, Tiny}]

Poniższa pomysłowa funkcja dodaje 2π tam gdzie dochodzi do przeskoku.

unwrap[list_]:=Module[{first=First[list]},
Accumulate@Prepend[Mod[Differences[list]+Pi,2 Pi]-Pi,first]
]; (* Brilliant, o3...,. absolutly brilliant. *)

ω=unwrap[varpi];

Przetworzone przez f. unwrap dane są już ciągłe.

ListPlot[ω[[1;;500]],Joined->True,PlotStyle->{Red,Medium}]

W bardzo długiej skali czasowej (50 mln lat na wykresie poniżej!) precesja faktycznie wydaje się jednostajna.

ListPlot[ω[[1;;50000]],Joined->True,PlotStyle->{Red,Medium}]

Spróbujemy do danych dopasować prostą.

dane=({t,ω}//Transpose)[[1;;50000]];

fit=FindFit[dane,precession*τ+ω0,{precession,ω0},τ]

Okres precesji orbity eliptycznej podręcznikowo to 112 tyś. lat.

2Pi/precession/.fit

Rzucimy okiem na tempo precesji policzone różnymi metodami dla różnych skal czasowych. Wynik zgadza się z ,,kanonem’’ tylko dla ostatnich 5000 lat.

ωF=Interpolation[Transpose[{t,ω}],InterpolationOrder->1];

eF=Interpolation[Transpose[{t,e}],InterpolationOrder->1];

Plot[{ωF[τ],0.04427862525854445`*τ,(π/56)τ},{τ,-10000,0}]

Plot[{ωF'[τ],eF[τ],ωF[τ]/τ,2Pi/112},{τ,-100,0}]

Plot[{ωF'[τ],eF[τ],ωF[τ]/τ,2Pi/112},{τ,-1000,0}]

Ponieważ proste metody analizy nie zadziałały, użyjemy analizy Fouriera. Dane dla e oraz  i wydają się quasi-periodyczne. Widać dudnienia, co sugeruje nałożone 2 bliskie częstości.

(* Beware of e->0 epochs! *)

ListPlot[e[[1;;5000]],Joined->True,InterpolationOrder->0,PlotStyle->{Blue,Thin}]

ListPlot[i[[1;;5000]],Joined->True,InterpolationOrder->0,PlotStyle->{Blue,Thin}]

Mimośród oscyluje od wartości niemal zerowych (3e-5) do około 0.07.

Norm/@({k,h}//Transpose)//MinMax

ScientificForm[%]

Typowy hodograf ,,wektora Rungego-Lenza’’.

ListPlot[({k,h}//Transpose)[[1;;1000]],Joined->True,PlotRange->0.1{{-1,1},{-1,1}}]

JAk wyżej, dla ,,montu pędu’’.

ListPlot[({
p,q}//Transpose)[[1;;1000]],Joined->True,PlotRange->0.05{{-1,1},{-1,1}}]

Przystępujemy do identyfikacja częstości w oscylacjach mimośrodu. Posiłkuję się przykładem z Manuala (Fourier->Applications->Frequency identification).

Length[e]

Rozważamy 20 milionów lat. W innych okresach mod dominujący może być inny.

n = 20000;
start=1;
pdata = e[[start;;n-1+start]];   (* e - ekscentryczność, i - inklinacja orbity *)
(*pdata=Sin[ω[[start;;n-1+start]]];*)
ListPlot[pdata,Joined->True,AspectRatio->1/3,ImageSize->900]

f = Abs[Fourier[pdata]];
peakMax= Last[TakeLargest[f,2]]; 
peaksize= Last[TakeLargest[f,4]]; (* Bierzemy 4 największe piki *)
peaks= Flatten[Position[f,x_ /; x≥ peaksize]]

peaksize

Show[ListPlot[f], Graphics[Table[{Red, Point[{peaks[[ii]], f[[peaks[[ii]]]]}]},{ii,2,4}]], PlotRange->{0,peakMax}]

f[[1;;12]]

Show[ListPlot[f], Graphics[Table[{Red,PointSize[Large], Point[{peaks[[ii]], f[[peaks[[ii]]]]}]},{ii,2,4}]], PlotRange->{{0,n/20},{0,peakMax}}]

Tu wybieramy mod do dokłądniejszej analizy, interesuje nas 2 lub trzeci pik.

pos=peaks[[3]]

fr = Abs[Fourier[pdata Exp[2 Pi I (pos - 2) N[Range[0, n - 1]]/n], FourierParameters->{0,2/n}]];
frpos = Position[fr, Max[fr]][[1,1]]

Show[ListPlot[fr], Graphics[{Red, Point[{frpos, fr[[frpos]]}]}], PlotRange->All]

N[n/(pos - 2 + 2 (frpos - 1)/n)]

Kod powyżej pozwala na identyfikację częstości:
1. Okresy oscylacji mimośrodu e:  405.186 tyś lat oraz 94.8685 tyś lat. W przybliżeniu 405 oraz 95 tyś. lat.
2. Okres oscylacji nachylenia orbity 109.948 tyś. lat oraz 97.9367 tyś. lat. W przybliżeniu 110 oraz 98 tyś. lat.

Poniżej znajduje się próba uzyskania podobnych jak wyżej danych orbitalnych wprost z efemeryd Ukłądu Słonecznego. Sięgają one ok. 10 000 lat wstecz i służą do kalibracji warunków początkowych dla numerycznej lub analitycznej symulacji oscylacji parametrów orbitalnych. Dane poniżej są równoważne ok. 10 punktom  symulacji Lascar (2010).

(*1. Choose a time grid.Negative=years BP*)
startYear=-8998;           (*yrs before 1950 CE*)
finishYear=2025;              (*“present”*)
step=25;             (*sampling interval*)

yrs=Range[startYear,finishYear,step];
dates=DateObject/@Map[{#,1,1}&,yrs];

ResourceFunction["HorizonsEphemerisData"]["OrbitalElements",{(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"]),<|"Dates"->(DateObject["Sat 1 Jan 1 BC", DateObject[{-1,1,1},"Day"]])|>},All]

ResourceFunction["HorizonsEphemerisData"]["OrbitalElements",{(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"]),<|"Dates"->(DateObject["Fri 1 Jan 8999 BC", DateObject[{-8999,1,1},"Day"]])|>},All]

ResourceFunction["HorizonsEphemerisData"]["OrbitalElements",{(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"]),<|"Dates"->(DateObject["Fri 1 Jan 8999 BC", DateObject[{-8999,1,1},"Day"]])|>},All,"Association"]

Length[dates]

e=Table[ResourceFunction["HorizonsEphemerisData"]["OrbitalElements",{(Entity["Earth", Entity["Planet","Earth"], "Entity[\"Planet\", \"Earth\"]", "planet"]),<|"Dates"->dates[[ii]]|>},All][[1;;,3]],{ii,1,Length[dates],1}]

ListPlot[Flatten[e],Joined->True,InterpolationOrder->0]