Comment découvrir Neptune depuis l'orbite d'Uranus (par simulation informatique)

11

Je voudrais démontrer l'existence d'une autre planète (Neptune) en étudiant l'écart entre l'observation de l'orbite d'Uranus et la prédiction mathématique, ce travail a été réalisé à partir du Verrier et j'aimerais comprendre sa méthode.

J'ai lu le chapitre 2, "La découverte de Neptune (1845-1846)", dans la biographie Le Verrier - Magnificent and Detestable Astronomer, mais il est trop approfondi et je n'ai pas très bien compris son travail.

J'étudie le problème des trois corps (Soleil, Uranus, Neptune) via Matlab et le problème des deux corps (Soleil, Uranus) en prenant la condition initiale d'ici:

http://nssdc.gsfc.nasa.gov/planetary/factsheet/uranusfact.html

J'ai essayé cette méthode: j'ai mis Uranus dans le périhélie avec le Max. la vitesse orbitale et moi calculons le demi-grand axe, et il est plus précis que celui obtenu en mettant Uranus et Neptune dans le périhélie avec leur Max respectif. vitesse orbitale.

Voici une photo sympa faite avec Matlab: Voici une photo sympa

Quelqu'un peut-il m'aider? ce que je dois faire et avec quelles données je dois comparer ma prédiction? Même un simple lien pourrait être utile.

Sergio Piccione
la source

Réponses:

11

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:

le système solaire

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:

Deux Uranus, avec et sans Neptune

... 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:

Evolution temporelle de la distance euclidienne 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:

Séparation angulaire entre les deux Uranus

... 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:

Séparation angulaire entre les deux Uranus, laissant Jupiter et Saturne hors de l'équation

Bien qu'il y ait des différences, elles sont minuscules et surtout non pertinentes pour découvrir Neptune.

Rody Oldenhuis
la source
Brillante réponse!
zephyr
4

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.

barrycarter
la source