Cómo descubrir Neptuno desde la órbita de Urano (por simulación por computadora)

11

Me gustaría demostrar la existencia de otro planeta (Neptuno) estudiando la discrepancia entre la observación de la órbita de Urano y la predicción matemática, este trabajo fue realizado por Le Verrier y me gustaría entender su método.

He leído el Capítulo 2, "El descubrimiento de Neptuno (1845-1846)", en la biografía Le Verrier - Astrónomo magnífico y detestable, pero es demasiado profundo y no entendí muy bien su trabajo.

Estoy estudiando el problema de tres cuerpos (Sol, Urano, Neptuno) a través de Matlab y el problema de dos cuerpos (Sol, Urano) tomando la condición inicial a partir de aquí:

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

He intentado este método: puse a Urano en el Perihelio con el Max. velocidad orbital y calculo el eje semi-mayor, y es más preciso que el obtenido al poner a Urano y Neptuno en el Perihelio con sus respectivos Max. velocidad orbital.

Aquí una foto genial hecha con Matlab: Aquí una foto genial

Alguien puede ayudarme? ¿Qué tengo que hacer y con qué datos tengo que comparar mi predicción? Incluso un simple enlace podría ser útil.

Sergio Piccione
fuente

Respuestas:

11

Aquí esta lo que hice:

  • Según sus masas, es más seguro considerar inicialmente a Júpiter y Saturno, así como a Urano. También podría ser fructífero incluir la Tierra en el análisis, para obtener posiciones relativas, ángulos de observación, etc. Por lo tanto, consideraré:
    • Dom
    • Tierra
    • Júpiter
    • Saturno
    • Urano
    • Neptuno
  • Obtenga los parámetros gravitacionales estándar (μ) para todos ellos
  • Obtenga posiciones y velocidades iniciales a través de JPL / HORIZONS para todos estos planetas. Tenía algunos datos de J2000.5, así que solo usé los vectores de estado del 1 de enero de 2000 al mediodía.
  • Escriba un integrador de N-body con herramientas MATLAB incorporadas. Integre este sistema solar incompleto una vez sin Neptuno, y una vez con Neptuno incluido.
  • Analizar y comparar!

Entonces, aquí están mis datos y el integrador de N-body:

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

Aquí está la secuencia de comandos del controlador que usé para obtener algunas buenas tramas:

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 describiré aquí paso a paso.

Primero, un diagrama del sistema solar para verificar que el integrador de N-body funcione como debería:

el sistema solar

¡Agradable! Luego, quería ver la diferencia entre las posiciones de Urano con y sin la influencia de Neptuno. Entonces, extraje solo las posiciones de esos dos Urano, y los tracé:

Dos Uranos, con y sin Neptuno.

... eso no es útil. Incluso al hacer un gran acercamiento y rotar al máximo, esto no es una trama útil. Así que miré la evolución de la distancia euclidiana absoluta entre los dos Urano:

Evolución temporal de la distancia euclidiana entre los dos Urano

¡Eso empieza a parecerse más! ¡Aproximadamente 80 años después del comienzo de nuestro análisis, los dos Urano están separados por casi 6 millones de kilómetros!

Por grande que pueda parecer, en la escala más amplia de las cosas, esto podría ahogarse en el ruido cuando tomamos medidas aquí en la Tierra. Además, todavía no cuenta toda la historia, como veremos en un momento. Entonces, veamos la diferencia angular entre los vectores de observación, desde la Tierra hacia los dos Urano, para ver qué tan grande es ese ángulo, y si puede sobresalir por encima de los umbrales de error de observación:

Separación angular entre los dos uranos.

... whoa! Mucho más de 300 segundos de arco de diferencia, además de todo tipo de ondulaciones tambaleantes y ondulantes. Eso parece estar dentro de las capacidades de observación de la época (aunque no puedo encontrar una fuente confiable sobre esto tan rápido; ¿alguien?)

Solo por si acaso, también produje esa última trama dejando a Júpiter y Saturno fuera de escena. Aunque algunos teoría de la perturbación se había desarrollado en el 17 º y 18 º siglos, no fue muy bien desarrollada y dudo incluso de Le Verrier tomó Júpiter en cuenta (pero de nuevo, podría estar equivocado, por favor corríjanme si sabe más).

Entonces, aquí está la última trama sin Júpiter y Saturno:

Separación angular entre los dos Urano, dejando a Júpiter y Saturno fuera de la ecuación.

Aunque existen diferencias, son mínimas, y lo más importante, irrelevantes para descubrir Neptuno.

Rody Oldenhuis
fuente
Brillante respuesta!
zephyr
4

Si entiendo correctamente, ¿estás modelando la órbita de Urano como una elipse y quieres compararla con la órbita real de Urano perturbada por Neptuno? No tengo una respuesta, pero ¿Dónde puedo encontrar / visualizar las posiciones de los planetas / estrellas / lunas / etc.? explica cómo usar SPICE, HORIZONS y otras herramientas para encontrar la verdadera posición de Uranus en un momento dado + -15000 años a partir de ahora, incluyendo los parámetros elípticos que mejor se ajustan (usando la función "elementos orbitales" de HORIZONS).

Por supuesto, todo lo que hagas será "circular" en algún sentido, ya que la posición calculada por HORIZONTES de Urano en el pasado ya incluye las perturbaciones de Neptuno.

Si pudiera encontrar tablas de predicciones de posición de Urano o algo del pasado, podría tener algo.

Por cierto, no dude en ponerse en contacto conmigo (vea el perfil para más detalles) si este proyecto se extiende más allá de una pregunta de intercambio de pila.

barrycarter
fuente