31:T404,

The orbit of a planet that would orbit alone around the Sun would be an immutable ellipse having the Sun located at one of its foci. This orbit would be completely determined by the following parameters: (i) the semi-major axis ${tex`a`} of the ellipse, (ii) the eccentricity ${tex`e`} of the ellipse, (iii) the longitude of the periapsis ${tex`\varpi`} (i.e., longitude of the perihelion in the case of the solar system), (iv) the inclination ${tex`I`} of its plane, and (v) the longitude of the ascending node ${tex`\Omega`}.

The presence of several planets around the Sun makes the precise determination of their motions much less simple. According to Newton's law of universal graviation, like the Sun, the planets act on each other proportionally to their masses and in inverse proportion to the squares of their respective distances.

32:T415,

The notebook is organised as follows:

Secular perturbation theory: N-body system

end of this section

Secular evolution of the solar system

section II

plot the evolution of the orbits

Other useful information is given in the Appendices:

Laplace coefficients

Matrix computations

Appendix B

^-7

i`, name: "Planet", relative_mass: html`m_i / m₀`, semi_major_axis: html`a_i [AU]`, eccentricity: html`e_i`, angular_velocity: html`n_i [°/y]`, longitude_perihelion: html`ϖ_i [°]`, inclination: html`I_i [°]`, longitude_ascending_node: html`Ω_i [°]`, }, align: { index: "left", name: "left", relative_mass: "left", semi_major_axis: "left", eccentricity: "left", angular_velocity: "left", longitude_perihelion: "left", inclination: "left", longitude_ascending_node: "left", }, width: { index: 50, name: 80, relative_mass: 110, semi_major_axis: 80, eccentricity: 85, angular_velocity: 70, longitude_perihelion: 80, inclination: 70, longitude_ascending_node: 70, }, maxWidth: 860, maxHeight: 240, height: 240 })35:T496,tex.block` \textbf{A} = \begin{bmatrix} A_{11} & A_{12} & \dotsb & A_{17} & A_{18} \\ A_{21} & A_{22} & \dotsb & A_{27} & A_{28} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ A_{71} & A_{72} & \dotsb & A_{77} &A_{78} \\ A_{81} & A_{82} & \dotsb & A_{87} &A_{88} \end{bmatrix} = \begin{bmatrix} ${(A[0][0]*1E3).toFixed(3)}\times 10^{-3} & ${(A[0][1]*1E4).toFixed(2)}\times 10^{-4} & \dotsb & ${(A[0][6]*1E8).toFixed(1)} \times 10^{-8} & ${(A[0][7]*1E9).toFixed(1)} \times 10^{-9} \\ ${(A[1][0]*1E5).toFixed(2)} \times 10^{-5} & ${(A[1][1]*1E3).toFixed(3)} \times 10^{-3} & \dotsb & ${(A[1][6]*1E8).toFixed(1)} \times 10^{-8} & ${(A[1][7]*1E9).toFixed(1)} \times 10^{-9} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ ${(A[6][0]*1E12).toFixed(1)} \times 10^{-12} & ${(A[6][1]*1E10).toFixed(1)} \times 10^{-10} & \dotsb & ${(A[6][6]*1E4).toFixed(3)} \times 10^{-4} & ${(A[6][7]*1E5).toFixed(2)} \times 10^{-5} \\ ${(A[7][0]*1E13).toFixed(1)} \times 10^{-13} & ${(A[7][1]*1E11).toFixed(1)} \times 10^{-11} & \dotsb & ${(A[7][6]*1E5).toFixed(2)} \times 10^{-5} & ${(A[7][7]*1E4).toFixed(3)} \times 10^{-4} \end{bmatrix} \,, `36:T496,tex.block` \textbf{B} = \begin{bmatrix} B_{11} & B_{12} & \dotsb & B_{17} & B_{18} \\ B_{21} & B_{22} & \dotsb & B_{27} & B_{28} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ B_{71} & B_{72} & \dotsb & B_{77} &B_{78} \\ B_{81} & B_{82} & \dotsb & B_{87} &B_{88} \end{bmatrix} = \begin{bmatrix} ${(B[0][0]*1E3).toFixed(3)}\times 10^{-3} & ${(B[0][1]*1E4).toFixed(2)}\times 10^{-4} & \dotsb & ${(B[0][6]*1E7).toFixed(1)} \times 10^{-7} & ${(B[0][7]*1E7).toFixed(1)} \times 10^{-7} \\ ${(B[1][0]*1E5).toFixed(2)} \times 10^{-5} & ${(B[1][1]*1E3).toFixed(3)} \times 10^{-3} & \dotsb & ${(B[1][6]*1E6).toFixed(1)} \times 10^{-6} & ${(B[1][7]*1E7).toFixed(1)} \times 10^{-7} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ ${(B[6][0]*1E10).toFixed(1)} \times 10^{-10} & ${(B[6][1]*1E8).toFixed(1)} \times 10^{-8} & \dotsb & ${(B[6][6]*1E4).toFixed(3)} \times 10^{-4} & ${(B[6][7]*1E5).toFixed(2)} \times 10^{-5} \\ ${(B[7][0]*1E11).toFixed(1)} \times 10^{-11} & ${(B[7][1]*1E9).toFixed(1)} \times 10^{-9} & \dotsb & ${(B[7][6]*1E5).toFixed(2)} \times 10^{-5} & ${(B[7][7]*1E4).toFixed(3)} \times 10^{-4} \end{bmatrix} \,. `37:T434,tex.block` \begin{array}{rlrcrlr} g_{1} &=& ${(g[0]).toFixed(3)} &\textrm{''/year} \:,& \qquad\qquad f_{1} &=& ${(f[0]).toFixed(3)} &\textrm{''/year} \:,\\ g_{2} &=& ${(g[1]).toFixed(3)} &\textrm{''/year} \:,& f_{2} &=& ${(f[1]).toFixed(3)} &\textrm{''/year} \:,\\ g_{3} &=& ${(g[2]).toFixed(2)} &\textrm{''/year} \:,& f_{3} &=& ${(f[2]).toFixed(2)} &\textrm{''/year} \:,\\ g_{4} &=& ${(g[3]).toFixed(2)} &\textrm{''/year} \:,& f_{4} &=& ${(f[3]).toFixed(2)} &\textrm{''/year} \:,\\ g_{5} &=& ${(g[4]).toFixed(3)} &\textrm{''/year} \:,& f_{5} &=& ${(f[4]).toFixed(3).replace(/^-0.000$/, '0.000')} &\textrm{''/year} \:,\\ g_{6} &=& ${(g[5]).toFixed(2)} &\textrm{''/year} \:,& f_{6} &=& ${(f[5]).toFixed(2)} &\textrm{''/year} \:, \\ g_{7} &=& ${(g[6]).toFixed(3)} &\textrm{''/year} \:,& f_{7} &=& ${(f[6]).toFixed(3)} &\textrm{''/year} \:, \\ g_{8} &=& ${(g[7]).toFixed(4)} &\textrm{''/year} \:,& f_{8} &=& ${(f[7]).toFixed(4)} &\textrm{''/year} \:. \end{array} `38:T58b,| ${tex`\bold{e}_{ij}`} | *i* = 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | | :----: | :---- | :---- | :---- | :---- | :---- | :---- | :---- | :---- | | ** *j* = 1 ** | ${E[0][0]} | ${E[0][1]} | ${E[0][2]} | ${E[0][3]} | ${E[0][4]} | ${E[0][5]} | ${E[0][6]} | ${E[0][7]} | | ** 2 ** | ${E[1][0]} | ${E[1][1]} | ${E[1][2]} | ${E[1][3]} | ${E[1][4]} | ${E[1][5]} | ${E[1][6]} | ${E[1][7]} | | ** 3 ** | ${E[2][0]} | ${E[2][1]} | ${E[2][2]} | ${E[2][3]} | ${E[2][4]} | ${E[2][5]} | ${E[2][6]} | ${E[2][7]} | | ** 4 ** | ${E[3][0]} | ${E[3][1]} | ${E[3][2]} | ${E[3][3]} | ${E[3][4]} | ${E[3][5]} | ${E[3][6]} | ${E[3][7]} | | ** 5 ** | ${E[4][0]} | ${E[4][1]} | ${E[4][2]} | ${E[4][3]} | ${E[4][4]} | ${E[4][5]} | ${E[4][6]} | ${E[4][7]} | | ** 6 ** | ${E[5][0]} | ${E[5][1]} | ${E[5][2]} | ${E[5][3]} | ${E[5][4]} | ${E[5][5]} | ${E[5][6]} | ${E[5][7]} | | ** 7 ** | ${E[6][0]} | ${E[6][1]} | ${E[6][2]} | ${E[6][3]} | ${E[6][4]} | ${E[6][5]} | ${E[6][6]} | ${E[6][7]} | | ** 8 ** | ${E[7][0]} | ${E[7][1]} | ${E[7][2]} | ${E[7][3]} | ${E[7][4]} | ${E[7][5]} | ${E[7][6]} | ${E[7][7]} |39:T58b,| ${tex`\bold{I}_{ij}`} | *i* = 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | | :----: | :---- | :---- | :---- | :---- | :---- | :---- | :---- | :---- | | ** *j* = 1 ** | ${I[0][0]} | ${I[0][1]} | ${I[0][2]} | ${I[0][3]} | ${I[0][4]} | ${I[0][5]} | ${I[0][6]} | ${I[0][7]} | | ** 2 ** | ${I[1][0]} | ${I[1][1]} | ${I[1][2]} | ${I[1][3]} | ${I[1][4]} | ${I[1][5]} | ${I[1][6]} | ${I[1][7]} | | ** 3 ** | ${I[2][0]} | ${I[2][1]} | ${I[2][2]} | ${I[2][3]} | ${I[2][4]} | ${I[2][5]} | ${I[2][6]} | ${I[2][7]} | | ** 4 ** | ${I[3][0]} | ${I[3][1]} | ${I[3][2]} | ${I[3][3]} | ${I[3][4]} | ${I[3][5]} | ${I[3][6]} | ${I[3][7]} | | ** 5 ** | ${I[4][0]} | ${I[4][1]} | ${I[4][2]} | ${I[4][3]} | ${I[4][4]} | ${I[4][5]} | ${I[4][6]} | ${I[4][7]} | | ** 6 ** | ${I[5][0]} | ${I[5][1]} | ${I[5][2]} | ${I[5][3]} | ${I[5][4]} | ${I[5][5]} | ${I[5][6]} | ${I[5][7]} | | ** 7 ** | ${I[6][0]} | ${I[6][1]} | ${I[6][2]} | ${I[6][3]} | ${I[6][4]} | ${I[6][5]} | ${I[6][6]} | ${I[6][7]} | | ** 8 ** | ${I[7][0]} | ${I[7][1]} | ${I[7][2]} | ${I[7][3]} | ${I[7][4]} | ${I[7][5]} | ${I[7][6]} | ${I[7][7]} |3a:T61e,html`

${figure_e_1 }

${figure_I_1}

${figure_e_2}

${figure_I_2}

${figure_e_3}

${figure_I_3}

${figure_e_4}

${figure_I_4}

${figure_e_5}

${figure_I_5}

${figure_e_6}

${figure_I_6}

${figure_e_7}

${figure_I_7}

${figure_e_8}

${figure_I_8}

`3b:T499,A = { // Create an array M [8x8] let M = new Array(8); for (let i = 0; i < 8; i++) { M[i] = new Array(8); } // Get the angular velocities for the 8 planets let planet_velocity = []; for(let i = 0; i<8; i++) { planet_velocity[i]=planet_data(i+1).angular_velocity; } // Get the relative masses for the 8 planets let planet_relative_mass = []; for(let i = 0; i<8; i++) { planet_relative_mass[i]=planet_data(i+1).relative_mass; } // Create an array V1 [8x1] let V1 = new Array(8); // Create an array M1 [8x8] let M1 = new Array(8); for (let i = 0; i < 8; i++) { M1[i] = new Array(8); } for (let i = 0; i < 8; i++) { for (let j = 0; j < 8; j++) { M1[i][j] = L_alpha[i][j]*L_alpha_bar[i][j]*L_1[i][j]; } } V1=math.multiply(M1,planet_relative_mass); // Save A_ij in M for (let i = 0; i < 8; i++) { for (let j = 0; j < 8; j++) { M[i][j] = (planet_velocity[i]/4)/((1+planet_relative_mass[j])); if(i!=j){ M[i][j] = - M[i][j]*planet_relative_mass[j]*L_alpha[i][j]*L_alpha_bar[i][j]*L_2[i][j]; } else { M[i][j] = M[i][j]*V1[i]; } } } return M; }3c:T49a,B = { // Create an array M [8x8] let M = new Array(8); for (let i = 0; i < 8; i++) { M[i] = new Array(8); } // Get the angular velocities for the 8 planets let planet_velocity = []; for(let i = 0; i<8; i++) { planet_velocity[i]=planet_data(i+1).angular_velocity; } // Get the relative masses for the 8 planets let planet_relative_mass = []; for(let i = 0; i<8; i++) { planet_relative_mass[i]=planet_data(i+1).relative_mass; } // Create an array V1 [8x1] let V1 = new Array(8); // Create an array M1 [8x8] let M1 = new Array(8); for (let i = 0; i < 8; i++) { M1[i] = new Array(8); } for (let i = 0; i < 8; i++) { for (let j = 0; j < 8; j++) { M1[i][j] = L_alpha[i][j]*L_alpha_bar[i][j]*L_1[i][j]; } } V1=math.multiply(M1,planet_relative_mass); // Save A_ij in M for (let i = 0; i < 8; i++) { for (let j = 0; j < 8; j++) { M[i][j] = (planet_velocity[i]/4)/((1+planet_relative_mass[j])); if(i!=j){ M[i][j] = M[i][j]*planet_relative_mass[j]*L_alpha[i][j]*L_alpha_bar[i][j]*L_1[i][j]; } else { M[i][j] = - M[i][j]*V1[i]; } } } return M; }3d:T499,

Build matrices ${tex`\textbf{A}`} and ${tex`\textbf{B}`} using data stored in arrays ${tex`\textbf{L}_{\alpha}`}, ${tex`\textbf{L}_{\bar{\alpha}}`}, ${tex`\textbf{L}_{1}`}, and ${tex`\textbf{L}_{2}`},

Compute the eigenvalues ${tex`\lambda_{j}`} and ${tex`\mu_{l}`} of matrices ${tex`\textbf{A}`} and ${tex`\textbf{B}`}, respectively,

Reorder the eigenvalues and define matrices ${tex`\textbf{M}_{\lambda}`} and ${tex`\textbf{M}_{\mu}`},

Store the ordered eigenvalues ${tex`g_{i}`} and ${tex`f_{i}`}, and compute the intermediate matrices ${tex`\tilde{\textbf{E}}`} and ${tex`\tilde{\textbf{I}}`},

Compute ${tex`\textbf{x}_{1}`} and ${tex`\textbf{x}_{2}`},

Find the intermediate angular phases ${tex`\tilde{\beta_{j}}`} and, ${tex`\tilde{\gamma_{j}}`},

Find the scaling factors ${tex`S_{j}`} and ${tex`T_{j}`},

Compute ${tex`\textbf{E}`} and ${tex`\textbf{I}`} (with positive terms on the diagonals) using ${tex`S_{j}`}, ${tex`T_{j}`}, ${tex`\tilde{\textbf{E}}`} and ${tex`\tilde{\textbf{I}}`},

Compute the angular phases ${tex`\beta_{j}`} and, ${tex`\gamma_{j}`}.

3e:T6d7,// Eccentricity e(t) for planet No. i = number function plot_e(number) { var data = [] var t, y const num_data = data_points for (let i = 0; i < num_data; i++) { t = -period/2 + i*period/num_data; y = math.pow( math.pow( E[0][number]*1e-5*math.sin((g[0]*t*math.PI/180/3600 + beta[0]*math.PI/180)) + E[1][number]*1e-5*math.sin((g[1]*t*math.PI/180/3600 + beta[1]*math.PI/180)) + E[2][number]*1e-5*math.sin((g[2]*t*math.PI/180/3600 + beta[2]*math.PI/180)) + E[3][number]*1e-5*math.sin((g[3]*t*math.PI/180/3600 + beta[3]*math.PI/180)) + E[4][number]*1e-5*math.sin((g[4]*t*math.PI/180/3600 + beta[4]*math.PI/180)) + E[5][number]*1e-5*math.sin((g[5]*t*math.PI/180/3600 + beta[5]*math.PI/180)) + E[6][number]*1e-5*math.sin((g[6]*t*math.PI/180/3600 + beta[6]*math.PI/180)) + E[7][number]*1e-5*math.sin((g[7]*t*math.PI/180/3600 + beta[7]*math.PI/180)) ,2) + math.pow( E[0][number]*1e-5*math.cos((g[0]*t*math.PI/180/3600 + beta[0]*math.PI/180)) + E[1][number]*1e-5*math.cos((g[1]*t*math.PI/180/3600 + beta[1]*math.PI/180)) + E[2][number]*1e-5*math.cos((g[2]*t*math.PI/180/3600 + beta[2]*math.PI/180)) + E[3][number]*1e-5*math.cos((g[3]*t*math.PI/180/3600 + beta[3]*math.PI/180)) + E[4][number]*1e-5*math.cos((g[4]*t*math.PI/180/3600 + beta[4]*math.PI/180)) + E[5][number]*1e-5*math.cos((g[5]*t*math.PI/180/3600 + beta[5]*math.PI/180)) + E[6][number]*1e-5*math.cos((g[6]*t*math.PI/180/3600 + beta[6]*math.PI/180)) + E[7][number]*1e-5*math.cos((g[7]*t*math.PI/180/3600 + beta[7]*math.PI/180)) ,2) ,0.5) data.push({time:t,eccentricity:y}) } return data } 3f:T6f1,// Inclination I(t) for planet No. i = number function plot_I(number) { var data = [] var t, y const num_data = data_points for (let i = 0; i < num_data; i++) { t = -period/2 + i*period/num_data; y = math.pow( math.pow( I[0][number]*1e-5*math.sin((f[0]*t*math.PI/180/3600 + gamma[0]*math.PI/180)) + I[1][number]*1e-5*math.sin((f[1]*t*math.PI/180/3600 + gamma[1]*math.PI/180)) + I[2][number]*1e-5*math.sin((f[2]*t*math.PI/180/3600 + gamma[2]*math.PI/180)) + I[3][number]*1e-5*math.sin((f[3]*t*math.PI/180/3600 + gamma[3]*math.PI/180)) + I[4][number]*1e-5*math.sin((f[4]*t*math.PI/180/3600 + gamma[4]*math.PI/180)) + I[5][number]*1e-5*math.sin((f[5]*t*math.PI/180/3600 + gamma[5]*math.PI/180)) + I[6][number]*1e-5*math.sin((f[6]*t*math.PI/180/3600 + gamma[6]*math.PI/180)) + I[7][number]*1e-5*math.sin((f[7]*t*math.PI/180/3600 + gamma[7]*math.PI/180)) ,2) + math.pow( I[0][number]*1e-5*math.cos((f[0]*t*math.PI/180/3600 + gamma[0]*math.PI/180)) + I[1][number]*1e-5*math.cos((f[1]*t*math.PI/180/3600 + gamma[1]*math.PI/180)) + I[2][number]*1e-5*math.cos((f[2]*t*math.PI/180/3600 + gamma[2]*math.PI/180)) + I[3][number]*1e-5*math.cos((f[3]*t*math.PI/180/3600 + gamma[3]*math.PI/180)) + I[4][number]*1e-5*math.cos((f[4]*t*math.PI/180/3600 + gamma[4]*math.PI/180)) + I[5][number]*1e-5*math.cos((f[5]*t*math.PI/180/3600 + gamma[5]*math.PI/180)) + I[6][number]*1e-5*math.cos((f[6]*t*math.PI/180/3600 + gamma[6]*math.PI/180)) + I[7][number]*1e-5*math.cos((f[7]*t*math.PI/180/3600 + gamma[7]*math.PI/180)) ,2) ,0.5)*180/math.PI data.push({time:t,inclination:y}) } return data } 40:Tae0,orbital_elements = [ { index: 1, name: "Mercury", relative_mass: 0.000000166004, // Relative to the Sun semi_major_axis: 0.38709893, // in AU angular_velocity: 1494.72939930044, // °/year eccentricity: 0.20563069, longitude_perihelion: 77.45645, // ° inclination: 7.00487, // ° longitude_ascending_node: 48.33167, }, { index: 2, name: "Venus", relative_mass: 0.000002447577, // Relative to the Sun semi_major_axis: 0.72333199, // in AU angular_velocity: 585.266193967071, // °/year eccentricity: 0.00677323, longitude_perihelion: 131.53298, // ° inclination: 3.39471, // ° longitude_ascending_node: 76.68069, }, { index: 3, name: "Earth-Moon", relative_mass: 0.000003003152, // Relative to the Sun semi_major_axis: 1.00000011, // in AU angular_velocity: 359.993728565271, // °/year eccentricity: 0.01671022, longitude_perihelion: 102.94719, // ° inclination: 0.00005, // ° longitude_ascending_node: 348.73936, }, { index: 4, name: "Mars", relative_mass: 0.000000322682, // Relative to the Sun semi_major_axis: 1.52366231, // in AU angular_velocity: 191.402952051006, // °/year eccentricity: 0.09341233, longitude_perihelion: 336.04084, // ° inclination: 1.85061, // ° longitude_ascending_node: 49.57854, }, { index: 5, name: "Jupiter", relative_mass: 0.000954497237, // Relative to the Sun semi_major_axis: 5.20336301, // in AU angular_velocity: 30.3490521835669, // °/year eccentricity: 0.04839266, longitude_perihelion: 14.75385, // ° inclination: 1.30530, // ° longitude_ascending_node: 100.55615, }, { index: 6, name: "Saturn", relative_mass: 0.000285786105, // Relative to the Sun semi_major_axis: 9.53707032, // in AU angular_velocity: 12.2211461425642, // °/year eccentricity: 0.05415060, longitude_perihelion: 92.43194, // ° inclination: 2.48446, // ° longitude_ascending_node: 113.71504, }, { index: 7, name: "Uranus", relative_mass: 0.000043653694, // Relative to the Sun semi_major_axis: 19.19126393, // in AU angular_velocity: 4.28466689476514, // °/year eccentricity: 0.04716771, longitude_perihelion: 170.96424, // ° inclination: 0.76986, // ° longitude_ascending_node: 74.22988, }, { index: 8, name: "Neptune", relative_mass: 0.000051495392, // Relative to the Sun semi_major_axis: 30.06896348, // in AU angular_velocity: 2.18440069773237, // °/year eccentricity: 0.00858587, longitude_perihelion: 44.97135, // ° inclination: 1.76917, // ° longitude_ascending_node: 131.72169, } ]41:[{"id":"f6a96ffd47b300eeb618aa449f0104147796cfc0d73da844ae204ce56edc3ae822bbd52a171d61448d510a43fef1b845b298aec14873fdd8d1f58754a93ade8e","href":"https://static.observableusercontent.com/files/f6a96ffd47b300eeb618aa449f0104147796cfc0d73da844ae204ce56edc3ae822bbd52a171d61448d510a43fef1b845b298aec14873fdd8d1f58754a93ade8e","name":"eccentricity_evolution.png","type":"image/png","lastModified":1663093893016.077,"size":142948}] 2f:{"title":"Secular perturbation theory","cells":[{"id":0,"value":"
Secular perturbation theory
\n
(Perturbation theory in celestial mechanics)
","pinned":false,"hidden":false,"mode":"html"},{"id":4276,"value":"
This is a continuation of a previous Observable notebook where we studied the secular evolution of the simplified 3-body system consisting of Jupiter, Saturn and the Sun.
","pinned":false,"hidden":false,"mode":"html"},{"id":1761,"value":"$31","pinned":false,"hidden":false,"mode":"html"},{"id":4221,"value":"introduction = md`
When considered over very long periods of time – over centuries or millions of years – interplanetary gravitational interactions have a profound influence on planetary orbits. In this notebook, a theory of celestial mechanics known as the orbital perturbation theory is used to examine the secular evolution of the eigth major planets of the solar system. This theory is also commonly referred to as Laplace-Lagrange secular perturbation theory.^{[[1](#Murray_Dermott),[2](#Fitzpatrick)]}
\n
\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":9515,"value":"eccentricity_secular_evolution1 = FileAttachment(\"eccentricity_evolution.png\").image()","pinned":false,"hidden":false,"mode":"ojs"},{"id":9462,"value":"
Secular evolution of the eccentricities of Mercury, Venus, Earth, Mars, Jupiter, and Saturn over a period of 20 million years. These plots were obtained thanks to the secular perturbation theory discussed in this notebook. \n
\n
[ llustration taken from Figure 2 ]\n
\n
","pinned":false,"hidden":false,"mode":"html"},{"id":4330,"value":"$32","pinned":false,"hidden":false,"mode":"html"},{"id":1973,"value":"three_body = section(\"I.\", \"Secular perturbation theory: N-body system\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":1320,"value":"
Our solar system consists of the Sun (our star, hereafter identified with index 0), and eight major planets: \n
\n
Mercury
\n
Venus
\n
Earth
\n
Mars
\n
Jupiter
\n
Saturn
\n
Uranus
\n
Neptune
\n
\n
\n
\nJupiter (planet No. 5) and Saturn (No. 6) are by far the heaviest planets of the solar system. If one limits the study to the 3-body problem, it is therefore appropriate to focus on these two planets since they are comparatively much less perturbed by the other planets. The simplified 3-body system is illustrated in Figure 1. \n
","pinned":false,"hidden":false,"mode":"html"},{"id":1041,"value":"$33","pinned":false,"hidden":false,"mode":"ojs"},{"id":2028,"value":"
\n
Figure 1: Model of the simplified 3-body system consisting of the Sun (⊙), Jupiter (♃), and Saturn (♄).
\n

","pinned":false,"hidden":false,"mode":"html"},{"id":4290,"value":"evolution_equations = md `
Evolution equations for the N-body system
\n
The perturbation theory developed for the simplified 3-body problem can be generalised to the N-body problem. Hereafter, we will pick up only some of the theoretical elements discussed in that notebook. The reader interested in the detailed theoretical developments is referred to the handbooks listed in the [References](#references) (e.g, Murray & Dermott, Chapter 6).^{[[1](#Murray_Dermott),[2](#Fitzpatrick)]}
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2158,"value":" md `
Disturbing function
\nThe N-body problem can be resolved by analyzing the gradient of a perturbing potential, also called the disturbing function, ${tex`\\mathcal{R}_{i}`}.\nThis disturbing function can be expanded into an infinite series whose terms can be classified either as secular, resonant, or short period, depending on the physical problem under consideration.^{[[1](#Murray_Dermott)]} Here, we will investigate the influence of interplanetary gravitational interactions on planetary orbits when integrated over long periods of time (i.e., secular perturbations).\n
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2185,"value":"
For this purpose, it is convenient to introduce alternative orbital\nelements for the eccentricity and inclination:
","pinned":false,"hidden":false,"mode":"html"},{"id":2187,"value":"tex.block`\n h_{i} = e_{i} \\sin(\\varpi_{i}) \\:, \\quad k_{i} = e_{i} \\cos(\\varpi_{i}) \\:, \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2195,"value":"tex.block`\n p_{i} = I_{i} \\sin(\\Omega_{i}) \\:, \\quad q_{i} = I_{i} \\cos(\\Omega_{i}) \\:. \\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2199,"value":"
These variables relate the eccentricity ${tex`e_{i}`} to the longitude of the perihelion ${tex`\\varpi_{i}`}, and the inclination ${tex`I_{i}`} to the longitude of the ascending node ${tex`\\Omega_{i}`} for each planet ${tex`i`} (where ${tex`i = 1...8`}).
","pinned":false,"hidden":false,"mode":"html"},{"id":2213,"value":"
It is also helpful to define the following parameters:
","pinned":false,"hidden":false,"mode":"html"},{"id":2215,"value":"tex.block`\n \\alpha_{ij} = \n\\begin{cases}\n \\frac{a_{i}}{a_j} \\quad &\\textrm{if} \\quad a_{i}a_{j} \\; (i>j) \\\\\n\\end{cases} \\:, \n\\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2235,"value":"
and
","pinned":false,"hidden":false,"mode":"html"},{"id":2237,"value":"tex.block`\n \\bar{\\alpha}_{ij} = \n\\begin{cases}\n \\frac{a_{i}}{a_j} \\quad &\\textrm{if} \\quad a_{i}a_{j} \\; (i>j) \\\\\n\\end{cases} \\:, \n\\\\\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2243,"value":"
where ${tex`a_{i}`} and ${tex`a_{j}`} are the semi-major axes of planets ${tex`i`} and ${tex`j`} (${tex`i,j = 1 ... 8`}, ${tex`i \\neq j`}).
","pinned":false,"hidden":false,"mode":"html"},{"id":2247,"value":"
The general secular part of the disturbing function can be written:
","pinned":false,"hidden":false,"mode":"html"},{"id":2250,"value":"tex.block`\n \\mathcal{R}_{i} = n_{i}a_{i}^2 \\Big[\\frac{1}{2}A_{ii}(h_{i}^2 + k_{i}^2)+ A_{ij}(h_{i}h_{j} + k_{i}k_{j}) \\\\\n \\qquad \\qquad \\qquad + \\frac{1}{2}B_{ii}(p_{i}^2 + q_{i}^2)+ B_{ij}(p_{i}p_{j} + q_{i}q_{j} \\Big] \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2269,"value":"
where ${tex`n_{i}`} is the orbital angular velocity of planet ${tex`i`}. The remaining quantities are constant elements from two matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`}:
","pinned":false,"hidden":false,"mode":"html"},{"id":2275,"value":"tex.block`\n A_{ii} = \\sum_{j, \\: j \\neq i} \\frac{n_i}{4}\\frac{m_{j}}{m_{0}+m_{i}}\\alpha_{ij}\\bar{\\alpha}_{ij} b_{3/2}^{(1)}(\\alpha_{ij}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2293,"value":"tex.block`\n A_{ij} = -\\frac{n_i}{4}\\frac{m_{j}}{m_{0}+m_{i}}\\alpha_{ij}\\bar{\\alpha}_{ij} b_{3/2}^{(2)}(\\alpha_{ij}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2298,"value":"tex.block`\n B_{ii} = -\\sum_{j, \\: j \\neq i} \\frac{n_i}{4}\\frac{m_{j}}{m_{0}+m_{i}}\\alpha_{ij}\\bar{\\alpha}_{ij} b_{3/2}^{(1)}(\\alpha_{ij}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2305,"value":"tex.block`\n B_{ij} = \\frac{n_i}{4}\\frac{m_{j}}{m_{0}+m_{i}}\\alpha_{ij}\\bar{\\alpha}_{ij} b_{3/2}^{(1)}(\\alpha_{ij}) \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2311,"value":"md `
In these equations, one can recognize the Laplace coefficients ${tex`b_s^{(m)}(\\alpha)`} given in [Appendix A](#Laplace_coefficients).
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2335,"value":"eigenfrequencies = md `
Eigenvalues and eigenvectors
The solutions of the equations of motions are
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2341,"value":"tex.block`\n h_{i} = \\sum_{j} e_{ij} \\sin(g_{j}t + \\beta_{j}) \\:, \\quad k_{i} = \\sum_{j} e_{ij} \\cos(g_{j}t + \\beta_{j}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2361,"value":"tex.block`\n p_{i} = \\sum_{j} I_{ij} \\sin(f_{j}t + \\gamma_{j}) \\:, \\quad q_{i} = \\sum_{j} I_{ij} \\cos(f_{j}t + \\gamma_{j}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2346,"value":"
where ${tex`t \\;`} is the time variable, \n
${tex`g_{j} \\;`} are the eigenvalues of matrix ${tex`\\textbf{A}`},\n
${tex`e_{ij} \\;`} are the eigenvectors of matrix ${tex`\\textbf{A}`},\n
${tex`f_{j} \\;`} are the eigenvalues of matrix ${tex`\\textbf{B}`},\n
${tex`I_{ij} \\;`} are the eigenvectors of matrix ${tex`\\textbf{B}`},\n
and ${tex`\\beta_{j}`} and ${tex`\\gamma_{j}`} are the phase angles.\n
","pinned":false,"hidden":false,"mode":"html"},{"id":2380,"value":"
The determination of the secular evolution of planetary orbits is then reduced to a numerical problem consisting in finding the eigenvalues and eigenvectors of matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`}. The final solution can be found by solving these equations for the initial condition (${tex`t=0`}) using the observation values for each planet ${tex`i`}: ${tex`e_{i}`}, ${tex`\\varpi_{i}`}, ${tex`I_{i}`}, and ${tex`\\Omega_{i}`}. All the details for matrix computations can be found in Appendix B.

","pinned":false,"hidden":false,"mode":"html"},{"id":1978,"value":"secular_perturbations = section(\"II.\", \"Secular evolution of the solar system\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":2396,"value":"md `
In a previous notebook, the study was limited to the 3-body system consisting of the Sun, Jupiter and Saturn. \nWe will now apply the secular perturbation theory to the 8 major planets orbiting the Sun. \nThe parameters used below are those of Julian year 2000, J2000 (see Murray & Dermott, Appendix \\\"Solar system data\\\"^{[[1](#Murray_Dermott)]}). In 2000 the orbital parameters were as follows:
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":2420,"value":"$34","pinned":false,"hidden":false,"mode":"ojs"},{"id":4242,"value":"
Note that we use the Earth-Moon system rather than planet Earth alone.
","pinned":false,"hidden":false,"mode":"html"},{"id":5454,"value":"
\n
II.1. Computation results
\n
All the computations details can be found in Appendix B.\nIn this section, we will simply summarize the main results.
","pinned":false,"hidden":false,"mode":"html"},{"id":6244,"value":"
The elements of matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`} are:
","pinned":false,"hidden":false,"mode":"html"},{"id":6246,"value":"$35","pinned":false,"hidden":false,"mode":"ojs"},{"id":6320,"value":"$36","pinned":false,"hidden":false,"mode":"ojs"},{"id":6439,"value":"
The eigenvalues of these matrices are:
","pinned":false,"hidden":false,"mode":"html"},{"id":6441,"value":"$37","pinned":false,"hidden":false,"mode":"ojs"},{"id":6635,"value":"
The phase angles are:
","pinned":false,"hidden":false,"mode":"html"},{"id":6638,"value":"tex.block`\n\\begin{array}{rlrcrlr}\n \\beta_{1} &=& ${(beta[0]).toFixed(1)} °, & & \\qquad\\qquad \\gamma_{1} &=& ${(gamma[0]).toFixed(1)} °,\\\\\n \\beta_{2} &=& ${(beta[1]).toFixed(1)} °, & & \\gamma_{2} &=& ${(gamma[1]).toFixed(1)} °,\\\\\n \\beta_{3} &=& ${(beta[2]).toFixed(1)} °, & & \\gamma_{3} &=& ${(gamma[2]).toFixed(1)} °,\\\\\n \\beta_{4} &=& ${(beta[3]).toFixed(1)} °, & & \\gamma_{4} &=& ${(gamma[3]).toFixed(1)} °,\\\\\n \\beta_{5} &=& ${(beta[4]).toFixed(1)} °, & & \\gamma_{5} &=& ${(gamma[4]).toFixed(1).replace(/^-0.000$/, '0.000')} °,\\\\\n \\beta_{6} &=& ${(beta[5]).toFixed(1)} °, & & \\gamma_{6} &=& ${(gamma[5]).toFixed(1)} °, \\\\\n \\beta_{7} &=& ${(beta[6]).toFixed(1)} °, & & \\gamma_{7} &=& ${(gamma[6]).toFixed(1)} °, \\\\\n \\beta_{8} &=& ${(beta[7]).toFixed(1)} °, & & \\gamma_{8} &=& ${(gamma[7]).toFixed(1)} °.\n\\end{array}\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6552,"value":"
The components of the eigenvectors, ${tex`e_{ij}`} and ${tex`I_{ij}`}, are listed in Tables 1 and 2.
","pinned":false,"hidden":false,"mode":"html"},{"id":6517,"value":"$38","pinned":false,"hidden":false,"mode":"md"},{"id":6580,"value":"
\n
Table 1: Components of the eccentricity eigenvectors, ${tex`e_{ij}`}. \n
All quantities have been multiplied by a factor ×10⁵: \n e.g., ${tex`e_{11}`} = ${(E[0][0]*1E-5).toFixed(5)}.
\n
\n
","pinned":false,"hidden":false,"mode":"html"},{"id":6617,"value":"$39","pinned":false,"hidden":false,"mode":"md"},{"id":6613,"value":"
\n
Table 2: Components of the inclination eigenvectors, ${tex`I_{ij}`}. \n
All quantities have been multiplied by a factor ×10⁵, values are in radian: \n e.g., ${tex`I_{11}`} = ${(I[0][0]*1E-5).toFixed(5)} rad.
\n
","pinned":false,"hidden":false,"mode":"html"},{"id":6674,"value":"md`\n
\n
For a comparison, we recommend that you consult the values provided by Murray & Dermott (Chapter 7).^{[[1](#Murray_Dermott)]}
\n
\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":823,"value":"graphs = md`
II.2. Secular evolution of the orbits of the 8 major planets
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":3549,"value":"
In section II.1, we have listed all the parameters needed to calculate ${tex`h_{i}(t)`}, ${tex`k_{i}(t)`}, ${tex`p_{i}(t)`}, and ${tex`q_{i}(t)`} for the 8 major planets. We can now derive the orbital elements of these planets at any time ${tex`t`}. For example, for the eccentricites we have
","pinned":false,"hidden":false,"mode":"html"},{"id":3560,"value":"tex.block`\n e_{i}(t) = \\sqrt{h_{i}^2(t)+k_{i}^2(t)} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":9443,"value":"
In Figure 2, we show the evolution of the variation of the orbital eccentricity and inclination\nof the 8 major planets of the solar system as a function of time.
\n
Although we must keep in mind that these are only approximations of the true behavior, these results provide many striking information about the long-term variations of the planetary orbits. One of these being that there is a noticeable coupling \nbetween Earth and Venus which show similar variations in their eccentricities and\ninclinations. Jupiter and Saturn also show a strong coupling, as we have already highlighted in the previous notebook covering the 3-body system (Jupiter, Saturn and the Sun).
","pinned":false,"hidden":false,"mode":"html"},{"id":9432,"value":"plot_secular_evolution = html`
\n
\n
Figure 2: Evolution of the orbital eccentricity ${tex`e_{i}(t)`} and inclination ${tex`I_{i}(t)`} of the 8 major planets of the solar system over a period of ${(math.round(period*1e-7,2)*10).toFixed(1)} million years centred on year 2000.
\n
\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":9550,"value":"viewof period = Inputs.range([500000, 50000000], {value: 20600000, step: 100000, label: \"Period (years):\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":9574,"value":"viewof data_points = Inputs.range([500, 5000], {value: 1000, step: 500, label: \"Points per graph: ⊛\"})","pinned":false,"hidden":false,"mode":"ojs"},{"id":9624,"value":"⊛ Be aware of the the aliasing effect when examining some of the graphs.
In particular, the periodic variations of Jupiter and Saturn are on a timescale of ∼ 50,000 years.\n
(This is a relatively \"high frequency\" for an observation period on the order of a few million years.)","pinned":false,"hidden":false,"mode":"html"},{"id":9458,"value":"$3a","pinned":false,"hidden":false,"mode":"ojs"},{"id":1874,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":1982,"value":"references = section(\"\", \"References\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":626,"value":"Murray_Dermott = reference(\"evolution_equations\", 1, \"C. D. Murray & S. F. Dermott, \\\"Solar System Dynamics,\\\" Cambridge University Press, 2000. DOI: 10.1017/CBO9781139174817\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":636,"value":"Fitzpatrick = reference(\"introduction\", 2, \"R. Fitzpatrick, \\\"An Introduction to Celestial Mechanics,\\\" Cambridge University Press, 2012. DOI: 10.1017/CBO9781139152310\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":152,"value":"Izsak1963 = reference(\"Laplace_series_expansion\", 3, \"I. G. Izsak, \\\"Laplace coefficients and their Newcomb derivatives,\\\" SAO Special Report 129, 1963. [PDF]\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":5331,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":1841,"value":"Laplace_coefficients = section(\"Appendix A:\", \"Laplace coefficients\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":764,"value":"md `
Laplace coefficients are key functions that need to be manipulated intensively when investigating the orbit of planets in the framework of the perturbation theory of planetary motion.^{[[1](#Murray_Dermott),[2](#Fitzpatrick)]} These coefficients look as follows:
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":12,"value":"tex.block`\nb_s^{(m)}(\\alpha) = \\frac{1}{\\pi} \\int_{0}^{2\\pi} \\frac{\\cos(m\\psi)}{ \\big( 1 - 2\\alpha\\cos(\\psi) + \\alpha^2 \\big)^s} \\,d\\psi \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":37,"value":"where ${tex`m`} is a natural number, \n
${tex`s`} is a half-integer,\n
and ${tex`\\alpha`} is a positive real number < 1.","pinned":false,"hidden":false,"mode":"html"},{"id":109,"value":"Laplace_series_expansion = md `
The series expansion of ${tex`b_s^{(m)}(\\alpha)`} is ^{[[3](#Izsak1963)]}:
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":116,"value":"tex.block`\nb_s^{(m)}(\\alpha) = 2\\alpha^m \\sum_{p=0}^{\\infin} \\; \\frac{\\Gamma(s+p)\\;\\Gamma(s+m+p)}{[\\Gamma(s)]^2 \\; p! \\; (m+p)!} \\; \\alpha^{2p} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5328,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":4955,"value":"matrix_computations = section(\"Appendix B:\", \"Matrix computations\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":4464,"value":"matrice = md`
The N-body problem finally reduces to the solving of a numerical problem consisting in finding the eigenvalues and the eigenvectors of the two matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`}. With 8 major planets in the solar system, we thus have 64 components in each matrix … For obvious reasons, it is useful to define a number of intermediate arrays, matrices and vectors. The dimensions of these objects are indicated in square brackets as [rows × columns].
`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4484,"value":"Let's define two vectors to group ${tex`h_{i},k_{i}`}, and ${tex`p_{i},q_{i}`}:","pinned":false,"hidden":false,"mode":"html"},{"id":4450,"value":"tex.block`\n\\textbf{v}_{0} = \\begin{bmatrix}\n h_{i}(0) \\\\\n k_{i}(0)\n\\end{bmatrix}_{[16 \\times 1]} \\:, \\qquad\n\\textbf{w}_{0} = \\begin{bmatrix}\n p_{i}(0) \\\\\n q_{i}(0)\n\\end{bmatrix}_{[16 \\times 1]} \\:, \\qquad i = 1...8 \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5472,"value":"v_0 = {\n let v_0 = [];\n for(let i = 0; i<8; i++)\n {\n // Store h_1...h_8 in rows v_0[0]...v_0[7]\n v_0[i]=planet_data(i+1).eccentricity*math.sin(planet_data(i+1).longitude_perihelion*math.PI/180);\n // Store k_1...k_8 in rows v_0[8]...v_0[15]\n v_0[i+8]=planet_data(i+1).eccentricity*math.cos(planet_data(i+1).longitude_perihelion*math.PI/180);\n }\n return v_0;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5497,"value":"w_0 = {\n let w_0 = [];\n for(let i = 0; i<8; i++)\n {\n // Store p_1...p_8 in rows w_0[0]...v_0[7]\n w_0[i]=planet_data(i+1).inclination*math.sin(planet_data(i+1).longitude_ascending_node*math.PI/180);\n // Store q_1...q_8 in rows w_0[8]...w_0[15]\n w_0[i+8]=planet_data(i+1).inclination*math.cos(planet_data(i+1).longitude_ascending_node*math.PI/180);\n }\n return w_0;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4499,"value":"
For parameters ${tex`\\alpha_{ij}`} and ${tex`\\bar{\\alpha}_{ij}`}, we define two arrays:
","pinned":false,"hidden":false,"mode":"html"},{"id":4506,"value":"tex.block`\n\\textbf{L}_{\\alpha} = \\begin{bmatrix}\n \\alpha_{ij} \\\\\n\\end{bmatrix}_{[8 \\times 8]} \\:, \\qquad\n\\textbf{L}_{\\bar{\\alpha}} = \\begin{bmatrix}\n \\bar{\\alpha_{ij}} \\\\\n\\end{bmatrix}_{[8 \\times 8]} \\:, \\qquad i,j = 1...8 \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5554,"value":"L_alpha = {\n \n // Get the semi-major axes for the 8 planets\n let planet_semi_major_axis = []; \n for(let i = 0; i<8; i++) { planet_semi_major_axis[i]=planet_data(i+1).semi_major_axis; }\n \n // Create an array L_alpha [8x8]\n let L_alpha = new Array(8);\n for (let i = 0; i < 8; i++) {\n L_alpha[i] = new Array(8);\n }\n \n // Save alpha_ij in L_alpha\n for (let i = 0; i < 8; i++) {\n for (let j = 0; j < 8; j++) {\n i!=j? L_alpha[i][j] = math.min(planet_semi_major_axis[i]/planet_semi_major_axis[j],planet_semi_major_axis[j]/planet_semi_major_axis[i]) : L_alpha[i][j] = 0; \n }\n }\n return L_alpha\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5589,"value":"L_alpha_bar = {\n \n // Get the semi-major axes for the 8 planets\n let planet_semi_major_axis = []; \n for(let i = 0; i<8; i++) { planet_semi_major_axis[i]=planet_data(i+1).semi_major_axis; }\n \n // Create an array L_alpha_bar [8x8]\n let L_alpha_bar = new Array(8);\n for (let i = 0; i < 8; i++) {\n L_alpha_bar[i] = new Array(8);\n }\n \n // Save alpha_ij in L_alpha\n for (let i = 0; i < 8; i++) {\n for (let j = 0; j < 8; j++) {\n i!=j? L_alpha_bar[i][j] = math.min(1,planet_semi_major_axis[i]/planet_semi_major_axis[j]) : L_alpha_bar[i][j] = 0; \n }\n }\n return L_alpha_bar\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4525,"value":"
For Laplace coefficients, we define two other arrays:
","pinned":false,"hidden":false,"mode":"html"},{"id":4529,"value":"tex.block`\n\\textbf{L}_{1} = \\begin{bmatrix}\n b_{3/2}^{(1)}(\\alpha_{ij}) \\\\\n\\end{bmatrix}_{[8 \\times 8]} \\:, \\qquad\n\\textbf{L}_{2} = \\begin{bmatrix}\n b_{3/2}^{(1)}(\\alpha_{ij}) \\\\\n\\end{bmatrix}_{[8 \\times 8]} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5604,"value":"L_1 = {\n \n // Create an array L [8x8]\n let L = new Array(8);\n for (let i = 0; i < 8; i++) {\n L[i] = new Array(8);\n }\n \n // Save alpha_ij in L_alpha\n for (let i = 0; i < 8; i++) {\n for (let j = 0; j < 8; j++) {\n i!=j? L[i][j] = b_m_s_alpha(1,1.5,L_alpha[i][j])[0] : L[i][j] = 0; \n }\n }\n return L\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5618,"value":"L_2 = {\n \n // Create an array L [8x8]\n let L = new Array(8);\n for (let i = 0; i < 8; i++) {\n L[i] = new Array(8);\n }\n \n // Save alpha_ij in L_alpha\n for (let i = 0; i < 8; i++) {\n for (let j = 0; j < 8; j++) {\n i!=j? L[i][j] = b_m_s_alpha(2,1.5,L_alpha[i][j])[0] : L[i][j] = 0; \n }\n }\n return L\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5622,"value":"
Matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`} can now be built easily using data stored in arrays\n ${tex`\\textbf{L}_{\\alpha}`}, ${tex`\\textbf{L}_{\\bar{\\alpha}}`}, ${tex`\\textbf{L}_{1}`}, and ${tex`\\textbf{L}_{2}`}.
","pinned":false,"hidden":false,"mode":"html"},{"id":5627,"value":"$3b","pinned":false,"hidden":false,"mode":"ojs"},{"id":5699,"value":"$3c","pinned":false,"hidden":false,"mode":"ojs"},{"id":5710,"value":"
The eigenvalues of matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`} are ${tex`\\lambda_{j}`} and ${tex`\\mu_{l}`}, respectively (${tex`j,l=1...8`}).
","pinned":false,"hidden":false,"mode":"html"},{"id":5704,"value":"lambda = {\n const prec = 1e-12 // Increase prec to ensure convergence of the eigenvalues \n const ans = math.eigs(A,prec) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}\n let lambda = math.multiply(ans.values,3600); // eigenvalues in ''/y\n return lambda;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5712,"value":"mu = {\n const prec = 1e-15 // Increase prec to ensure convergence of the eigenvalues\n const ans = math.eigs(B,prec) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}\n let mu = math.multiply(ans.values,3600); // eigenvalues in ''/y\n return mu\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4537,"value":"
When computing eigenvalues of a matrix, the ordering of the eigenvalues depends on the way the matrix is computed. It is therefore needed to manually sort the eigenvalues and match the corresponding eigenvectors. To perform these operations, we define two sorting matrices, ${tex`\\textbf{M}_{\\lambda}`} and ${tex`\\textbf{M}_{\\mu}`}:
","pinned":false,"hidden":false,"mode":"html"},{"id":4577,"value":"tex.block`\n\\begin{bmatrix}\n \\lambda_{j} \\\\\n\\end{bmatrix}_{[8 \\times 1]} \\xrightarrow{\\textbf{M}_{\\lambda}}\n\\begin{bmatrix}\n g_{i} \\\\\n\\end{bmatrix}_{[8 \\times 1]} \\;, \\qquad\n\\begin{bmatrix}\n g_{i} \\\\\n\\end{bmatrix} = \\textbf{M}_{\\lambda} \n\\begin{bmatrix}\n \\lambda_{j} \\\\\n\\end{bmatrix} \\:, \\qquad i,j = 1...8 \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4590,"value":"tex.block`\n\\begin{bmatrix}\n \\mu_{l} \\\\\n\\end{bmatrix}_{[8 \\times 1]} \\xrightarrow{\\textbf{M}_{\\mu}}\n\\begin{bmatrix}\n f_{i} \\\\\n\\end{bmatrix}_{[8 \\times 1]} \\;, \\qquad\n\\begin{bmatrix}\n f_{i} \\\\\n\\end{bmatrix} = \\textbf{M}_{\\mu} \n\\begin{bmatrix}\n \\mu_{l} \\\\\n\\end{bmatrix} \\:, \\qquad i,l = 1...8 \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4617,"value":"
Here, the indices ${tex`j`} and ${tex`l`} are chosen distinct to emphasize the fact that the eigenvalues must be sorted for the two matrices ${tex`\\textbf{A}`} and ${tex`\\textbf{B}`}.
","pinned":false,"hidden":false,"mode":"html"},{"id":5733,"value":"M_lambda = {\n\n // Create an array M filled with zeroes [8x8]\n let M = new Array(8).fill(0);\n for (let i = 0; i < 8; i++) {\n M[i] = new Array(8).fill(0);\n }\n \n M[0][7] = 1; // λ1 -> g8\n M[1][6] = 1; // λ2 -> g7\n M[2][4] = 1; // λ3 -> g5\n M[3][0] = 1; // λ4 -> g1\n M[4][1] = 1; // λ5 -> g2\n M[5][2] = 1; // λ6 -> g3\n M[6][3] = 1; // λ7 -> g4\n M[7][5] = 1; // λ8 -> g6\n\n return M;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5753,"value":"g = {\n let g = [];\n g = math.multiply(math.transpose(M_lambda),lambda);\n return g;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5759,"value":"M_mu = {\n\n // Create an array M filled with zeroes [8x8]\n let M = new Array(8).fill(0);\n for (let i = 0; i < 8; i++) {\n M[i] = new Array(8).fill(0);\n }\n \n M[0][4] = 1; // μ1 -> f5\n M[1][7] = 1; // μ2 -> f8\n M[2][6] = 1; // μ3 -> f7\n M[3][0] = 1; // μ4 -> f1\n M[4][1] = 1; // μ5 -> f2\n M[5][3] = 1; // μ6 -> f4\n M[6][2] = 1; // μ7 -> f3\n M[7][5] = 1; // μ8 -> f6\n\n return M;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5764,"value":"f = {\n let f = [];\n f = math.multiply(math.transpose(M_mu),mu);\n return f;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4597,"value":"
The corresponding eigenvectors are now sorted but they are still unscaled. \nIndeed, there is an infinity of eigenvectors that satisfy:","pinned":false,"hidden":false,"mode":"html"},{"id":4746,"value":"tex.block`\n\\textbf{A}\\begin{bmatrix}\n \\tilde{e}_{1j} \\\\\n \\vdots \\\\\n \\tilde{e}_{8j}\n\\end{bmatrix} = \\begin{bmatrix}\n \\tilde{e}_{1j} \\\\\n \\vdots \\\\\n \\tilde{e}_{8j}\n\\end{bmatrix} g_{j} \\,, \\qquad\n\\textbf{B}\\begin{bmatrix}\n \\tilde{I}_{1j} \\\\\n \\vdots \\\\\n \\tilde{I}_{8j}\n\\end{bmatrix} = \\begin{bmatrix}\n \\tilde{I}_{1j} \\\\\n \\vdots \\\\\n \\tilde{I}_{8j}\n\\end{bmatrix} f_{j} \\,, \\qquad j=1...8 \\,.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4765,"value":"
Therefore, we define the following intermediate matrices:
","pinned":false,"hidden":false,"mode":"html"},{"id":4599,"value":"tex.block`\n\\tilde{\\textbf{E}} = \\begin{bmatrix}\n \\tilde{e}_{ij} \\\\\n\\end{bmatrix}_{[8 \\times 8]} \\:, \\qquad\n\\tilde{\\textbf{I}} = \\begin{bmatrix}\n \\tilde{I}_{ij} \\\\\n\\end{bmatrix}_{[8 \\times 8]} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5785,"value":"E_tilde = {\n const prec = 1e-12 // Increase prec to ensure convergence of the eigenvalues \n const ans = math.eigs(A,prec) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}\n let E_tilde = math.transpose(math.multiply(math.transpose(M_lambda),math.transpose(ans.vectors))); // eigenvectors \n return E_tilde;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":5805,"value":"I_tilde = {\n const prec = 1e-15 // Increase prec to ensure convergence of the eigenvalues \n const ans = math.eigs(B,prec) // returns {values: [E1,E2...sorted], vectors: [v1,v2.... corresponding vectors as columns]}\n let I_tilde = math.transpose(math.multiply(math.transpose(M_mu),math.transpose(ans.vectors))); // eigenvectors\n return I_tilde;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4657,"value":"
To find the scaling factors, we must solve the equations of motions at time ${tex`t=0`}:
","pinned":false,"hidden":false,"mode":"html"},{"id":4727,"value":"tex.block`\n h_{i}(0) = \\sum_{j} e_{ij} \\sin(\\beta_{j}) \\:, \\quad k_{i}(0) = \\sum_{j} e_{ij} \\cos(\\beta_{j}) \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4733,"value":"tex.block`\n p_{i}(0) = \\sum_{j} I_{ij} \\sin(\\gamma_{j}) \\:, \\quad q_{i} = \\sum_{j} I_{ij} \\cos(\\gamma_{j}) \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4740,"value":"
Hence, we need to find the scaling factors ${tex`S_{j}`} and ${tex`T_{j}`} that satisfy
","pinned":false,"hidden":false,"mode":"html"},{"id":4781,"value":"tex.block`\n e_{ij} = S_{j}\\tilde{e}_{ij} \\:, I_{ij} = T_{j}\\tilde{I}_{ij} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4683,"value":"
One further request is that in the final matrices, the terms in the diagonal must be positive:
","pinned":false,"hidden":false,"mode":"html"},{"id":4684,"value":"tex.block`\n\\textbf{E} = \\begin{bmatrix}\n e_{ij} \\\\\n\\end{bmatrix}_{[8 \\times 8]} \\:, e_{ii} \\geq 0 \\:, \\qquad\n\\textbf{I} = \\begin{bmatrix}\n I_{ij} \\\\\n\\end{bmatrix}_{[8 \\times 8]} \\:, I_{ii} \\geq 0 \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4794,"value":"
To solve the problem, we introduce the following matrices and vectors:
","pinned":false,"hidden":false,"mode":"html"},{"id":4857,"value":"tex.block`\n\\textbf{v} = \\begin{bmatrix}\n h_{i}(0) \\\\\n k_{i}(0)\n\\end{bmatrix}_{[16 \\times 1]} =\n\\textbf{M}_{1} \\cdot \\textbf{x}_{1} \\;, \\qquad\n\n\\textbf{M}_{1} =\n\\begin{bmatrix}\n \\tilde{\\textbf{E}} & \\textbf{0} \\\\\n \\textbf{0} & \\tilde{\\textbf{E}}\n\\end{bmatrix}_{[16 \\times 16]} \\;, \\quad\n\n\\textbf{x}_{1} =\n\\begin{bmatrix}\n S_{1}\\sin(\\tilde{\\beta_{1}}) \\\\\n \\vdots \\\\\n S_{8}\\sin(\\tilde{\\beta_{8}}) \\\\\n S_{1}\\cos(\\tilde{\\beta_{1}}) \\\\\n \\vdots \\\\\n S_{8}\\cos(\\tilde{\\beta_{8}})\n\\end{bmatrix}_{[16 \\times 1]} \\:, \\quad\n\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4921,"value":"tex.block`\n\\textbf{w} = \\begin{bmatrix}\n p_{i}(0) \\\\\n q_{i}(0)\n\\end{bmatrix}_{[16 \\times 1]} =\n\\textbf{M}_{2} \\cdot \\textbf{x}_{2} \\;, \\qquad\n\n\\textbf{M}_{1} =\n\\begin{bmatrix}\n \\tilde{\\textbf{I}} & \\textbf{0} \\\\\n \\textbf{0} & \\tilde{\\textbf{I}}\n\\end{bmatrix}_{[16 \\times 16]} \\;, \\quad\n\n\\textbf{x}_{2} =\n\\begin{bmatrix}\n T_{1}\\sin(\\tilde{\\gamma_{1}}) \\\\\n \\vdots \\\\\n T_{8}\\sin(\\tilde{\\gamma_{8}}) \\\\\n T_{1}\\cos(\\tilde{\\gamma_{1}}) \\\\\n \\vdots \\\\\n T_{8}\\cos(\\tilde{\\gamma_{8}})\n\\end{bmatrix}_{[16 \\times 1]} \\:. \\quad\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":5941,"value":"M_1 = {\n // Create a matrix [16x16]\n // [[ E_tilde 0 ],\n // [ 0 E_tilde ]]\n let M = math.matrix()\n M.resize([16,16],0);\n \n M.subset(math.index([0,1,2,3,4,5,6,7],[0,1,2,3,4,5,6,7]), math.matrix(E_tilde));\n M.subset(math.index([8,9,10,11,12,13,14,15],[8,9,10,11,12,13,14,15]), math.matrix(E_tilde));\n\n return M._data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6097,"value":"M_2 = {\n // Create a matrix [16x16]\n // [[ I_tilde 0 ],\n // [ 0 I_tilde ]]\n let M = math.matrix()\n M.resize([16,16],0);\n \n M.subset(math.index([0,1,2,3,4,5,6,7],[0,1,2,3,4,5,6,7]), math.matrix(I_tilde));\n M.subset(math.index([8,9,10,11,12,13,14,15],[8,9,10,11,12,13,14,15]), math.matrix(I_tilde));\n\n return M._data\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":4929,"value":"
The solutions are found by solving
","pinned":false,"hidden":false,"mode":"html"},{"id":4900,"value":"tex.block`\n\\textbf{x}_{1} = \\big(\\textbf{M}_{1}^{\\textrm{T}}\\textbf{M}_{1}\\big)^{-1} \\textbf{M}_{1}^{\\textrm{T}} \\cdot \\textbf{v} \\:,\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":4912,"value":"tex.block`\n\\textbf{x}_{2} = \\big(\\textbf{M}_{2}^{\\textrm{T}}\\textbf{M}_{2}\\big)^{-1} \\textbf{M}_{2}^{\\textrm{T}} \\cdot \\textbf{w} \\:.\n`","pinned":false,"hidden":false,"mode":"ojs"},{"id":6101,"value":"x_1 = {\n let M = math.multiply(math.inv(math.multiply(math.transpose(M_1),M_1)),math.transpose(M_1));\n let x = math.multiply(M,v_0);\n return x;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6129,"value":"x_2 = {\n let M = math.multiply(math.inv(math.multiply(math.transpose(M_2),M_2)),math.transpose(M_2));\n let x = math.multiply(M,w_0)\n return x;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6137,"value":"beta_tilde = {\n let beta = [];\n for(let i = 0; i < 8; i++)\n {\n beta[i] = math.atan2(x_1[i],x_1[i+8]);\n }\n return beta;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6155,"value":"gamma_tilde = {\n let gamma = [];\n for(let i = 0; i < 8; i++)\n {\n gamma[i] = math.atan2(x_2[i],x_2[i+8]);\n }\n return gamma;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6158,"value":"S = {\n let S = [];\n for(let i = 0; i < 8; i++)\n {\n S[i] = math.pow(math.pow(x_1[i],2) + math.pow(x_1[i+8],2),0.5);\n }\n return S;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6162,"value":"T = {\n let T = [];\n for(let i = 0; i < 8; i++)\n {\n T[i] = math.pow(math.pow(x_2[i],2) + math.pow(x_2[i+8],2),0.5);\n }\n return T;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6166,"value":"E_beta = {\n // Create a matrix with scaled eigenvectors\n let M = math.transpose(math.multiply(E_tilde,math.diag(S)));\n\n // Get the signs of the terms on the diagonal (M_ii)\n let diag = math.sign(math.diag(M));\n let M_sign = math.diag(diag);\n\n // Transform M to have positive values on the diagonal\n M = math.multiply(M_sign,M)\n\n // Change beta accordingly\n let phase = [];\n for(let i = 0; i < 8; i++)\n {\n diag[i]<0 ? phase[i] = beta_tilde[i] + math.PI : phase[i] = beta_tilde[i];\n }\n \n return [M,phase]\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6190,"value":"E = math.round(math.multiply(E_beta[0],1E5),0)","pinned":false,"hidden":false,"mode":"ojs"},{"id":6192,"value":"beta = math.round(math.multiply(E_beta[1],180/math.PI),1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":6196,"value":"I_gamma = {\n // Create a matrix with scaled eigenvectors\n let M = math.transpose(math.multiply(I_tilde,math.diag(T)));\n\n // Get the signs of the terms on the diagonal (M_ii)\n let diag = math.sign(math.diag(M));\n let M_sign = math.diag(diag);\n\n // Transform M to have positive values on the diagonal\n M = math.multiply(M_sign,M)\n\n // Change gamma accordingly\n let phase = [];\n for(let i = 0; i < 8; i++)\n {\n diag[i]<0 ? phase[i] = gamma_tilde[i] + math.PI : phase[i] = gamma_tilde[i];\n }\n \n return [M,phase]\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":6200,"value":"I = math.round(math.multiply(I_gamma[0],1E5*math.PI/180),0)","pinned":false,"hidden":false,"mode":"ojs"},{"id":6203,"value":"gamma = math.round(math.multiply(I_gamma[1],180/math.PI),1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":4833,"value":"
\n
Summary of the computation steps
\n
\nTo perform the computation, we have followed the steps listed below:
","pinned":false,"hidden":false,"mode":"html"},{"id":4836,"value":"$3d","pinned":false,"hidden":false,"mode":"html"},{"id":10,"value":"
","pinned":false,"hidden":false,"mode":"html"},{"id":1914,"value":"other_appendix = section(\"Appendix C:\", \"Javascript / code + variables\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":873,"value":"math = require(\"mathjs@11.0.0\")","pinned":false,"hidden":false,"mode":"ojs"},{"id":197,"value":"function b_m_s_alpha (m,s,alpha, P_max = 100){ // P_max is an optional parameter\n let b_m_s_alpha = 0;\n let sum = 0;\n let sum_test = 0;\n let p = 0;\n for(p; p <= P_max && Number.isFinite(sum_test); p++){\n sum_test = math.gamma(s+p)/math.pow(math.gamma(s),2); \n sum_test /= math.factorial(p);\n sum_test *= math.gamma(s+m+p)*math.pow(alpha,2*p);\n sum_test /= math.factorial(m+p);\n Number.isFinite(sum_test)? sum += sum_test : null; // Check if the calculated value is a finite number\n }\n b_m_s_alpha = 2*math.pow(alpha,m)*sum;\n return [b_m_s_alpha,p-1];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":7256,"value":"$3e","pinned":false,"hidden":false,"mode":"ojs"},{"id":7298,"value":"$3f","pinned":false,"hidden":false,"mode":"ojs"},{"id":1838,"value":"section = (section, text) => {\n return html`
${section} ${text}
`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":147,"value":"// https://observablehq.com/@jashkenas/markdown-style-footnotes\nreference = (referencingCellName, number, text) => {\n return html`
↫ ${text}
`;\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2682,"value":"function planet_data (index) {\n return orbital_elements.filter(d => d.index==index)[0];\n}","pinned":false,"hidden":false,"mode":"ojs"},{"id":2413,"value":"$40","pinned":false,"hidden":false,"mode":"ojs"},{"id":6892,"value":"plot_height = 310","pinned":false,"hidden":false,"mode":"ojs"},{"id":6895,"value":"plot_width = 510","pinned":false,"hidden":false,"mode":"ojs"},{"id":6754,"value":"color_plot = \"#45526D\"","pinned":false,"hidden":false,"mode":"ojs"},{"id":5648,"value":"stroke_width = 1.2","pinned":false,"hidden":false,"mode":"ojs"},{"id":7146,"value":"label_years = \"→ Years \"","pinned":false,"hidden":false,"mode":"ojs"},{"id":7270,"value":"plot_e1 = plot_e(0)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7274,"value":"plot_e2 = plot_e(1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7280,"value":"plot_e3 = plot_e(2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7283,"value":"plot_e4 = plot_e(3)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7286,"value":"plot_e5 = plot_e(4)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7289,"value":"plot_e6 = plot_e(5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7292,"value":"plot_e7 = plot_e(6)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7295,"value":"plot_e8 = plot_e(7)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7312,"value":"plot_I1 = plot_I(0)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7315,"value":"plot_I2 = plot_I(1)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7317,"value":"plot_I3 = plot_I(2)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7319,"value":"plot_I4 = plot_I(3)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7321,"value":"plot_I5 = plot_I(4)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7323,"value":"plot_I6 = plot_I(5)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7325,"value":"plot_I7 = plot_I(6)","pinned":false,"hidden":false,"mode":"ojs"},{"id":7327,"value":"plot_I8 = plot_I(7)","pinned":false,"hidden":false,"mode":"ojs"},{"id":6677,"value":"figure_e_1 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0.115,0.255]), \n Plot.line(plot_e1, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_e1, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Mercury\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 0,\n dy: -24\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7112,"value":"figure_e_2 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0,0.07]), \n Plot.line(plot_e2, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_e2, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Venus\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 100,\n dy: 5\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7250,"value":"figure_e_3 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0,0.07]), \n Plot.line(plot_e3, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_e3, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Earth\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 80,\n dy: -10\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7403,"value":"figure_e_4 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"10.5px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0,0.16]), \n Plot.line(plot_e4, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_e4, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Mars\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 165,\n dy: -12\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7520,"value":"figure_e_5 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"10.5px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0.023,0.066]), \n Plot.line(plot_e5, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_e5, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Jupiter\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: -140,\n dy: -17\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7549,"value":"figure_e_6 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"10.5px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0.007,0.095]), \n Plot.line(plot_e6, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_e6, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Saturn\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: -120,\n dy: -18\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7610,"value":"figure_e_7 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"10.5px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0.00,0.09]), \n Plot.line(plot_e7, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_e7, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Uranus\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: -150,\n dy: -18\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7633,"value":"figure_e_8 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Eccentricity\"\n },\n style: { fontSize: \"10.5px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0.0035,0.0165]), \n Plot.line(plot_e8, Plot.windowY({k: 1, x: \"time\", y: \"eccentricity\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_e8, Plot.selectMaxY({\n x: \"time\",\n y: \"eccentricity\",\n text: () => \"Neptune\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 120,\n dy: -22\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":6732,"value":"figure_I_1 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([2.4,12.5]), \n Plot.line(plot_I1, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_I1, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Mercury\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 0,\n dy: -20\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7133,"value":"figure_I_2 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0,5.1]), \n Plot.line(plot_I2, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_I2, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Venus\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: -170,\n dy: -2\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7253,"value":"figure_I_3 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0,4.4]), \n Plot.line(plot_I3, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_I3, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Earth\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: -140,\n dy: 5\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7382,"value":"figure_I_4 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0,8]), \n Plot.line(plot_I4, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_I4, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Mars\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 80,\n dy: -8\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7537,"value":"figure_I_5 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([1.0,2.25]), \n Plot.line(plot_I5, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_I5, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Jupiter\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 70,\n dy: -18\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7563,"value":"figure_I_6 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0.42,2.86]), \n Plot.line(plot_I6, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_I6, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Saturn\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 35,\n dy: -15\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7616,"value":"figure_I_7 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0.25,3.15]), \n Plot.line(plot_I7, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_I7, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Uranus\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: -160,\n dy: -22\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"},{"id":7646,"value":"figure_I_8 = Plot.plot({\n height: plot_height,\n width: plot_width,\n x: {\n grid: true,\n label: label_years\n }, \n y: {\n grid: false,\n label: \"↑ Inclination (°)\"\n },\n style: { fontSize: \"11px\" },\n marks: [\n Plot.ruleX([-period/2,period/2]),\n Plot.ruleX([0], {strokeDasharray: \"4,3\"}),\n Plot.ruleY([0.64,2.68]), \n Plot.line(plot_I8, Plot.windowY({k: 1, x: \"time\", y: \"inclination\", strokeWidth:stroke_width, stroke: color_plot})),\n\n Plot.text(plot_I8, Plot.selectMaxY({\n x: \"time\",\n y: \"inclination\",\n text: () => \"Neptune\",\n fontSize: 13.5,\n fontWeight: 500,\n fill: color_plot,\n dx: 45,\n dy: -18\n })),\n ]\n})","pinned":false,"hidden":false,"mode":"ojs"}],"files":"$W41","stdlib":"1","titling":"auto"}