Voici ce que j'ai fait:
- Sur la base de leurs masses, il est plus sûr de considérer initialement Jupiter et Saturne ainsi que Uranus. Il pourrait également être utile d'inclure la Terre dans l'analyse, d'obtenir des positions relatives, des angles d'observation, etc. Donc, je considérerai:
- Soleil
- Terre
- Jupiter
- Saturne
- Uranus
- Neptune
- Obtenez les paramètres gravitationnels standard (μ) pour chacun d'eux
- Obtenez les positions et vitesses initiales via JPL / HORIZONS pour toutes ces planètes. J'avais quelques données qui traînaient à partir de J2000.5, donc je viens d'utiliser les vecteurs d'état du 1er janvier 2000 à midi.
- Écrivez un intégrateur à N corps avec les outils MATLAB intégrés. Intégrez ce système solaire incomplet une fois sans Neptune et une fois avec Neptune inclus.
- Analysez et comparez!
Voici donc mes données et l'intégrateur à N corps:
function [t, yout_noNeptune, yout_withNeptune] = discover_Neptune()
% Time of integration (in years)
tspan = [0 97] * 365.25 * 86400;
% std. gravitational parameters [km/s²/kg]
mus_noNeptune = [1.32712439940e11; % Sun
398600.4415 % Earth
1.26686534e8 % Jupiter
3.7931187e7 % Saturn
5.793939e6]; % Uranus
mus_withNeptune = [mus_noNeptune
6.836529e6]; % Neptune
% Initial positions [km] and velocities [km/s] on 2000/Jan/1, 00:00
% These positions describe the barycenter of the associated system,
% e.g., sJupiter equals the statevector of the Jovian system barycenter.
% Coordinates are expressed in ICRF, Solar system barycenter
sSun = [0 0 0 0 0 0].';
sEarth = [-2.519628815461580E+07 1.449304809540383E+08 -6.175201582312584E+02,...
-2.984033716426881E+01 -5.204660244783900E+00 6.043671763866776E-05].';
sJupiter = [ 5.989286428194381E+08 4.390950273441353E+08 -1.523283183395675E+07,...
-7.900977458946710E+00 1.116263478937066E+01 1.306377465321731E-01].';
sSaturn = [ 9.587405702749230E+08 9.825345942920649E+08 -5.522129405702555E+07,...
-7.429660072417541E+00 6.738335806405299E+00 1.781138895399632E-01].';
sUranus = [ 2.158728913593440E+09 -2.054869688179662E+09 -3.562250313222718E+07,...
4.637622471852293E+00 4.627114800383241E+00 -4.290473194118749E-02].';
sNeptune = [ 2.514787652167830E+09 -3.738894534538290E+09 1.904284739289832E+07,...
4.466005624145428E+00 3.075618250100339E+00 -1.666451179600835E-01].';
y0_noNeptune = [sSun; sEarth; sJupiter; sSaturn; sUranus];
y0_withNeptune = [y0_noNeptune; sNeptune];
% Integrate the partial Solar system
% once with Neptune, and once without
options = odeset('AbsTol', 1e-8,...
'RelTol', 1e-10);
[t, yout_noNeptune] = ode113(@(t,y) odefcn(t,y,mus_noNeptune) , tspan, y0_noNeptune , options);
[~, yout_withNeptune] = ode113(@(t,y) odefcn(t,y,mus_withNeptune), t, y0_withNeptune, options);
end
% The differential equation
%
% dy/dt = d/dt [r₀ v₀ r₁ v₁ r₂ v₂ ... rₙ vₙ]
% = [v₀ a₀ v₁ a₁ v₂ a₂ ... vₙ aₙ]
%
% with
%
% aₓ = Σₘ -G·mₘ/|rₘ-rₓ|² · (rₘ-rₓ) / |rₘ-rₓ|
% = Σₘ -μₘ·(rₘ-rₓ)/|rₘ-rₓ|³
%
function dydt = odefcn(~, y, mus)
% Split up position and velocity
rs = y([1:6:end; 2:6:end; 3:6:end]);
vs = y([4:6:end; 5:6:end; 6:6:end]);
% Number of celestial bodies
N = size(rs,2);
% Compute interplanetary distances to the power -3/2
df = bsxfun(@minus, permute(rs, [1 3 2]), rs);
D32 = permute(sum(df.^2), [3 2 1]).^(-3/2);
D32(1:N+1:end) = 0; % (remove infs)
% Compute all accelerations
as = -bsxfun(@times, mus.', D32); % (magnitudes)
as = bsxfun(@times, df, permute(as, [3 2 1])); % (directions)
as = reshape(sum(as,2), [],1); % (total)
% Output derivatives of the state vectors
dydt = y;
dydt([1:6:end; 2:6:end; 3:6:end]) = vs;
dydt([4:6:end; 5:6:end; 6:6:end]) = as;
end
Voici le script du pilote que j'ai utilisé pour sortir de jolis tracés:
clc
close all
% Get coordinates from N-body simulation
[t, yout_noNeptune, yout_withNeptune] = discover_Neptune();
% For plot titles etc.
bodies = {'Sun'
'Earth'
'Jupiter'
'Saturn'
'Uranus'
'Neptune'};
% Extract positions
rs_noNeptune = yout_noNeptune (:, [1:6:end; 2:6:end; 3:6:end]);
rs_withNeptune = yout_withNeptune(:, [1:6:end; 2:6:end; 3:6:end]);
% Figure of the whole Solar sysetm, just to check
% whether everything went OK
figure, clf, hold on
for ii = 1:numel(bodies)
plot3(rs_withNeptune(:,3*(ii-1)+1),...
rs_withNeptune(:,3*(ii-1)+2),...
rs_withNeptune(:,3*(ii-1)+3),...
'color', rand(1,3));
end
axis equal
legend(bodies);
xlabel('X [km]');
ylabel('Y [km]');
title('Just the Solar system, nothing to see here');
% Compare positions of Uranus with and without Neptune
rs_Uranus_noNeptune = rs_noNeptune (:, 13:15);
rs_Uranus_withNeptune = rs_withNeptune(:, 13:15);
figure, clf, hold on
plot3(rs_Uranus_noNeptune(:,1),...
rs_Uranus_noNeptune(:,2),...
rs_Uranus_noNeptune(:,3),...
'b.');
plot3(rs_Uranus_withNeptune(:,1),...
rs_Uranus_withNeptune(:,2),...
rs_Uranus_withNeptune(:,3),...
'r.');
axis equal
xlabel('X [km]');
ylabel('Y [km]');
legend('Uranus, no Neptune',...
'Uranus, with Neptune');
% Norm of the difference over time
figure, clf, hold on
rescaled_t = t/365.25/86400;
dx = sqrt(sum((rs_Uranus_noNeptune - rs_Uranus_withNeptune).^2,2));
plot(rescaled_t,dx);
xlabel('Time [years]');
ylabel('Absolute offset [km]');
title({'Euclidian distance between'
'the two Uranuses'});
% Angles from Earth
figure, clf, hold on
rs_Earth_noNeptune = rs_noNeptune (:, 4:6);
rs_Earth_withNeptune = rs_withNeptune(:, 4:6);
v0 = rs_Uranus_noNeptune - rs_Earth_noNeptune;
v1 = rs_Uranus_withNeptune - rs_Earth_withNeptune;
nv0 = sqrt(sum(v0.^2,2));
nv1 = sqrt(sum(v1.^2,2));
dPhi = 180/pi * 3600 * acos(min(1,max(0, sum(v0.*v1,2) ./ (nv0.*nv1) )));
plot(rescaled_t, dPhi);
xlabel('Time [years]');
ylabel('Separation [arcsec]')
title({'Angular separation between the two'
'Uranuses when observed from Earth'});
que je décrirai ici étape par étape.
Tout d'abord, un tracé du système solaire pour vérifier que l'intégrateur à N corps fonctionne comme il se doit:
Agréable! Ensuite, je voulais voir la différence entre les positions d'Uranus avec et sans l'influence de Neptune. J'ai donc extrait uniquement les positions de ces deux Uranus et les ai tracées:
... ce n'est guère utile. Même lorsque vous effectuez un zoom avant important et que vous le faites pivoter, ce n'est tout simplement pas un tracé utile. J'ai donc regardé l'évolution de la distance euclidienne absolue entre les deux Uranus:
Cela commence à y ressembler davantage! Environ 80 ans après le début de notre analyse, les deux Uranus sont distants de près de 6 millions de kilomètres!
Aussi grand que cela puisse paraître, à plus grande échelle, cela pourrait se noyer dans le bruit lorsque nous prenons des mesures ici sur Terre. De plus, il ne raconte toujours pas toute l'histoire, comme nous le verrons dans un instant. Alors maintenant, regardons la différence angulaire entre les vecteurs d'observation, de la Terre vers les deux Uranus pour voir à quel point cet angle est grand, et s'il peut se démarquer au-dessus des seuils d'erreur d'observation:
... whoa! Différence de plus de 300 secondes d'arc, en plus de toutes sortes d'ondulations temporelles ondulantes. Cela semble bien dans les capacités d'observation de l'époque (bien que je ne puisse pas trouver une source fiable à ce sujet si rapidement; n'importe qui?)
Juste pour faire bonne mesure, j'ai également produit cette dernière intrigue en laissant Jupiter et Saturne hors de l'image. Bien que certains théorie de la perturbation a été développée dans le 17 e et 18 e siècles, il n'a pas été très bien développé et je doute même Le Verrier a Jupiter en considération (mais encore une fois, je peux me tromper, s'il vous plaît me corriger si vous en savez plus).
Voici donc la dernière intrigue sans Jupiter et Saturne:
Bien qu'il y ait des différences, elles sont minuscules et surtout non pertinentes pour découvrir Neptune.
Si je comprends bien, vous modélisez l'orbite d'Uranus comme une ellipse et souhaitez la comparer à l'orbite réelle d'Uranus perturbée par Neptune? Je n'ai pas de réponse, mais où puis-je trouver / visualiser les positions des planètes / étoiles / lunes / etc? explique comment utiliser SPICE, HORIZONS et d'autres outils pour trouver la vraie position d'Uranus à un moment donné + -15000 ans à partir de maintenant, y compris les paramètres elliptiques les plus adaptés (en utilisant la fonction "éléments orbitaux" d'HORIZONS).
Bien sûr, tout ce que vous ferez sera "circulaire" dans un certain sens, puisque la position calculée par HORIZONS d'Uranus dans le passé inclut déjà les perturbations de Neptune.
Si vous pouviez trouver des tableaux de prévisions de position d'Uranus ou quelque chose du passé, vous pourriez avoir quelque chose.
BTW, n'hésitez pas à me contacter (voir profil pour plus de détails) si ce projet s'étend au-delà d'une question d'échange de pile.
la source