Вот что я сделал:
- Исходя из их массы, наиболее безопасно изначально рассмотреть Юпитер и Сатурн, а также Уран. Также было бы полезно включить Землю в анализ, чтобы получить относительные положения, углы наблюдения и т. Д. Итак, я буду рассматривать:
- солнце
- Земля
- Юпитер
- Сатурн
- Уран
- Нептун
- Получите стандартные гравитационные параметры (μ) для всех них
- Получить начальные позиции и скорости через JPL / ГОРИЗОНТЫ для всех этих планет. У меня были некоторые данные из J2000.5, поэтому я просто использовал векторы состояний с 1 января 2000 года в полдень.
- Напишите N-body интегратор со встроенными инструментами MATLAB. Интегрируйте эту неполную Солнечную систему один раз без Нептуна и один раз с включенным Нептуном.
- Анализируйте и сравнивайте!
Итак, вот мои данные и 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
Вот скрипт драйвера, который я использовал для получения хороших графиков:
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'});
который я опишу здесь шаг за шагом.
Во-первых, график Солнечной системы, чтобы проверить, что интегратор N-тела работает как следует:
Ницца! Затем я хотел увидеть разницу между позициями Урана с влиянием Нептуна и без него. Итак, я выделил только позиции этих двух Уранов и составил их график:
... это вряд ли полезно. Даже если сильно увеличить масштаб и повернуть его, это просто бесполезный сюжет. Итак, я посмотрел на эволюцию абсолютного евклидова расстояния между двумя Уранами:
Это начинает больше походить на это! Примерно через 80 лет после начала нашего анализа два Урана находятся на расстоянии почти 6 миллионов км!
Как бы это ни звучало, в более широком масштабе это может заглушиться шумом, когда мы проводим измерения здесь, на Земле. Плюс, это все еще не рассказывает всю историю, как мы увидим через мгновение. Итак, давайте посмотрим на угловую разницу между векторами наблюдения от Земли к двум Уранам, чтобы увидеть, насколько велик этот угол, и может ли он выделяться выше пороговых значений ошибок наблюдения:
... Вау! Разница более чем в 300 угловых секунд, плюс всякие шаткие волнистые волнистые волны. Это кажется вполне подходящим для наблюдательных возможностей того времени (хотя я не могу найти надежный источник по этому вопросу так быстро, кто-нибудь?)
Просто для примера, я также создал этот последний сюжет, оставив Юпитер и Сатурн вне поля зрения. Хотя некоторая теория возмущений была разработана в 17- м и 18- м веках, она была не очень хорошо развита, и я сомневаюсь, что даже Ле Верье принял во внимание Юпитер (но опять же, я могу ошибаться; пожалуйста, поправьте меня, если вы знаете больше).
Итак, вот последний сюжет без Юпитера и Сатурна:
Хотя есть различия, они ничтожны и, самое главное, не имеют отношения к открытию Нептуна.
Если я правильно понимаю, вы моделируете орбиту Урана как эллипс и хотите сравнить ее с реальной орбитой Урана, возмущенной Нептуном? У меня нет ответа, но где я могу найти / визуализировать положения планет / звезд / лун / и т.д.? объясняет, как использовать SPICE, HORIZONS и другие инструменты, чтобы найти истинное положение Урана в данный момент времени + -15000 лет, включая наиболее подходящие эллиптические параметры (с помощью функции HORIZONS "orbital elements").
Конечно, все, что вы будете делать, будет в некотором смысле «круговым», так как горизонтальное положение Урана в прошлом уже включает возмущения Нептуна.
Если бы вы могли найти таблицы предсказаний положения Урана или что-то из прошлого, у вас могло бы быть что-то.
Кстати, не стесняйтесь обращаться ко мне (см. Профиль для деталей), если этот проект выходит за рамки вопроса обмена стека.
источник