太阳地球月亮运动轨迹MATLAB仿真程序

太阳地球月亮运动轨迹MATLAB仿真程序太阳、地球、月亮运动轨迹MATLAB仿真程序
该代码包含了地球绕太阳运动，月亮绕地球运动的MATLAB轨迹仿真程序，实时显示地球、月亮的运动轨迹。

有兴趣的朋友可以购买，以供交流。

程序从第二页开始。

M文件1:draw_ball.m
% 画三维球体的函数
% (x0,y0,z0)为球心
% r为球半径
function draw_ball(x0, y0, z0, r) [x1, y1, z1]=sphere;
x = x1*r + x0;
y = y1*r + y0;
z = z1*r + z0;
surf(x,y,z);
M文件2:draw_circle.m
% 画二维圆形
% (x0,y0)为圆心
% r为圆半径
function draw_circle(x0, y0, r) theta = 0:pi/100:2*pi;
x = r*cos(theta)+x0;
y = r*sin(theta)+y0;
plot(x,y,'-r');
M文件3:RungeKutta_EarthSun.m % 四阶Runge-kutta法解地日微分方程
function yh = RungeKutta_EarthSun(w, h)
oldy = w;
k1 = dery(oldy);
midy = oldy + k1*h/2;
k2 = dery(midy);
midy = oldy + k2*h/2;
k3 = dery(midy);
midy = oldy + k3*h;
k4 = dery(midy);
yh = oldy + (k1 + 2*k2 + 2*k3 + k4)*h/6;
% 地日微分方程组
function dy = dery(w)
G = 6.674e-11; Ms = 1.989e30; miu_s = G * Ms; r2 = w(4)^2+w(5)^2;
dy(1) = 1;
dy(2) = (-1)*miu_s*w(4)/r2^1.5; dy(3) = (-1)*miu_s*w(5)/r2^1.5; dy(4) = w(2);
dy(5) = w(3);
M文件4:RungeKutta_MoonEarth.m % 四阶Runge-kutta法解地月微分方程
function yh = RungeKutta_MoonEarth(w, h)
oldy = w;
k1 = dery(oldy);
midy = oldy + k1*h/2;
k2 = dery(midy);
midy = oldy + k2*h/2;
k3 = dery(midy);
midy = oldy + k3*h;
k4 = dery(midy);
yh = oldy + (k1 + 2*k2 + 2*k3 + k4)*h/6;
% 地月微分方程组
function dy = dery(w)
G = 6.674e-11; Me = 5.972e24; miu_e = G * Me; r2 = w(4)^2+w(5)^2;
dy(1) = 1;
dy(2) = (-1)*miu_e*w(4)/r2^1.5; dy(3) = (-1)*miu_e*w(5)/r2^1.5; dy(4) = w(2);
dy(5) = w(3);
M文件5:sun_earth_moon.m
% 地月日全运动动态仿真程序
clc; clear all; close all;
G = 6.674e-11;% 引力常数
Ms = 1.989e30; Rs = 696300e3;% 太阳的质量和半径
Me = 5.972e24; Re = 6378e3;% 地球的质量和半径
Mm = 7.348e22; Rm = 3678e3;% 月球的质量和半径
sim_time = 3600*24*375;% 总仿真时间，375天，即约1年 h = 3600*24;%
仿真步长，24小时，即1天
w_e = [0 0 29535.6 1.52171522e11 0];% 地球相对太阳的初始位置 w_m =
[w_e(1) 0 990.32 405500e3 0];% 月球相对地球的初始位置
i = 1;% 用于记录解算的步数
while w_e(1) <= sim_time% 解算地球相对于太阳的轨迹
trajectory_earth2sun(:, i) = w_e;% trajectory_earth2sun存储地球相对于太阳的轨迹数据
ye = RungeKutta_EarthSun(w_e, h);% 4阶Runge-kutta法解算地日微分方程w_e = ye;% 用于Runge-kutta法解算微分方程的初值
i = i+1;% 解算一步则在步数上+1
end
i = 1;% 用于记录解算的步数
while w_m(1) <= sim_time% 解算月球相对于地球的轨迹
trajectory_moon2earth(:, i) = w_m;% trajectory_moon2earth存储月球相对于地球的轨迹数据
ym = RungeKutta_MoonEarth(w_m, h);% 4阶Runge-kutta法解算地月微分方程
w_m = ym;% 同上
i = i+1;% 同上
end
% 计算月球相对于太阳的轨迹，trajectory_moon2sun存储月球相对于太阳的轨迹数据 trajectory_moon2sun(1,:) = trajectory_moon2earth(4,:) + trajectory_earth2sun(4,:);
trajectory_moon2sun(2,:) = trajectory_moon2earth(5,:) +
trajectory_earth2sun(5,:);
size_enlge = 40;% 放大系数，用于放大太阳、地球、月球的半径，以便视觉观测
figure;% /01/画出地日三维动态轨迹
for j=1:i-1
draw_ball(0, 0, 0, Rs*size_enlge); hold on;% 画太阳的三维模型
title('地日-动态轨迹'); xlabel('x/m');ylabel('y/m');zlabel('z/m'); grid on;
plot(trajectory_earth2sun(4, :), trajectory_earth2sun(5, :), '-r');% 画地球相对于太阳的轨迹线
axis equal;
draw_ball(trajectory_earth2sun(4, 1), trajectory_earth2sun(5, 1), 0, Re*size_enlge^2);% 初始位置的地球三维模型
draw_ball(trajectory_earth2sun(4, j), trajectory_earth2sun(5, j), 0, Re*size_enlge^2);% 轨迹点上的地球三维模型
pause(0.001);hold off;
end
figure;% /02/画出地月三维动态轨迹
for j=1:i-1
draw_ball(0, 0, 0, Re*size_enlge/4); hold on;% 画地球的三维模型
title('地月-动态轨迹'); xlabel('x/m');ylabel('y/m');zlabel('z/m'); grid on;
plot(trajectory_moon2earth(4, :), trajectory_moon2earth(5, :), '-
r');% 画月球相对于地球的轨迹线
axis equal;
draw_ball(trajectory_moon2earth(4, 1), trajectory_moon2earth(5, 1), 0, Rm*size_enlge/4);%
初始位置的月球三维模型
draw_ball(trajectory_moon2earth(4, j), trajectory_moon2earth(5, j), 0, Rm*size_enlge/4);%
轨迹点上的月球三维模型
pause(0.001);hold off;
end
figure;% /03/画出地月日三维动态轨迹
for j=1:i-1
draw_ball(0, 0, 0, Rs*size_enlge); hold on;% 画太阳的三维模型
title('地月日-动态轨迹'); xlabel('x/m');ylabel('y/m');zlabel('z/m'); grid on;
plot(trajectory_earth2sun(4, :), trajectory_earth2sun(5, :), '-r');% 画地球相对于太阳的轨
迹线
plot(trajectory_moon2sun(1, :), trajectory_moon2sun(2, :), '-b');% 画月球相对于太阳的轨
迹线
axis equal;
draw_ball(trajectory_earth2sun(4, 1), trajectory_earth2sun(5, 1), 0, Re*size_enlge/2);% 初
始位置的地球三维模型
draw_ball(trajectory_moon2sun(1, 1), trajectory_moon2sun(2, 1), 0, Rm*size_enlge/2);%
初始位置的月球三维模型
draw_ball(trajectory_earth2sun(4, j), trajectory_earth2sun(5, j), 0, Re*size_enlge/2);% 轨迹
点上的地球三维模型
draw_ball(trajectory_moon2sun(1, j), trajectory_moon2sun(2, j), 0, Rm*size_enlge/2);% 轨
迹点上的月球三维模型
pause(0.001);hold off;
end
% 放大地月之间的相对位置，以便三维视觉显示
size_up = 15;% 放大倍数，用于放大地月之间的相对位置
% trajectory_moon2earth_up存储放大后的月球到地球的位置数据
trajectory_moon2earth_up(1,:) = size_up* (trajectory_moon2sun(1,:)-trajectory_earth2sun(4,:));
trajectory_moon2earth_up(2,:) = size_up* (trajectory_moon2sun(2,:)-trajectory_earth2sun(5,:));
% trajectory_moon2sun_up存储放大后的月球相对于太阳的轨迹数据
trajectory_moon2sun_up(1,:) =
(trajectory_earth2sun(4,:)+trajectory_moon2earth_up(1,:));
trajectory_moon2sun_up(2,:) =
(trajectory_earth2sun(5,:)+trajectory_moon2earth_up(2,:));
figure;% /04/画出放大后的地月日三维动态轨迹
for j=1:i-1
draw_ball(0, 0, 0, Rs*size_enlge); hold on;% 画太阳的三维模型
title('地月日-动态轨迹(放大)');
xlabel('x/m');ylabel('y/m');zlabel('z/m'); grid on;
plot(trajectory_earth2sun(4, :), trajectory_earth2sun(5, :), '-r');% 画地球相对于太阳的轨
迹线
plot(trajectory_moon2sun_up(1, :), trajectory_moon2sun_up(2, :), '-
b');% 画月球相对于太
阳的轨迹线
axis equal;
draw_ball(trajectory_earth2sun(4, 1), trajectory_earth2sun(5, 1), 0, Re*size_enlge*10);%
初始位置的地球三维模型
draw_ball(trajectory_moon2sun_up(1, 1), trajectory_moon2sun_up(2, 1), 0, Rm*size_enlge*10);% 初始位置的月球三维模型
draw_ball(trajectory_earth2sun(4, j), trajectory_earth2sun(5, j), 0, Re*size_enlge*10);% 轨
迹点上的地球三维模型
draw_ball(trajectory_moon2sun_up(1, j), trajectory_moon2sun_up(2, j), 0,
Rm*size_enlge*10);% 轨迹点上的月球三维模型
pause(0.001);hold off;
end
figure;% /05/画出地月日真实的二维动态轨迹
for j=1:i-1
draw_circle(0, 0, Rs*5); hold on;% 画太阳的二维模型
title('地月日二维真实轨迹');
xlabel('x/m');ylabel('y/m');zlabel('z/m'); grid on;
plot(trajectory_earth2sun(4, :), trajectory_earth2sun(5, :), '-r');% 画地球相对于太阳的轨迹线
plot(trajectory_moon2sun(1, :), trajectory_moon2sun(2, :), '-b');% 画月球相对于太阳的轨迹线
axis equal;
plot(trajectory_earth2sun(4, j), trajectory_earth2sun(5, j), 'r*');% 轨迹点上的地球点
plot(trajectory_moon2sun(1, j), trajectory_moon2sun(2, j), 'k*');% 轨迹点上的月球点
pause(0.001);hold off;
end
figure;% /06/画出地月日放大的二维动态轨迹
for j=1:i-1
draw_circle(0, 0, Rs*5); hold on;% 画太阳的二维模型
title('地月日二维放大轨迹');
xlabel('x/m');ylabel('y/m');zlabel('z/m'); grid on;
plot(trajectory_earth2sun(4, :), trajectory_earth2sun(5, :), '-r');% 画地球相对于太阳的轨迹线
plot(trajectory_moon2sun_up(1, :), trajectory_moon2sun_up(2, :), '-b');% 画月球相对于太阳的轨迹线
plot(trajectory_earth2sun(4, j), trajectory_earth2sun(5, j), 'r*');% 轨迹点上的地球点
plot(trajectory_moon2sun_up(1, j), trajectory_moon2sun_up(2, j),
'k*');% 轨迹点上的月球点
pause(0.001);hold off;
end
M文件6:sun_earth_moon_2Dsingle.m % 地月日全运动动态仿真程序
clc;clear all;close all;
G = 6.674e-11;
Ms = 1.989e30;Rs = 696300e3;
Me = 5.972e24;Re = 6378e3;
Mm = 7.348e22;Rm = 3678e3;
sim_time = 3600*24*375*3;
h = 3600*24;
i = 1;
w_e = [0 0 29535.6 1.52171522e11 0];
w_m = [w_e(1) 0 990.32 405500e3 0];
while w_e(1) <= sim_time
trajectory_earth2sun(:, i) = w_e;
ye = RungeKutta_EarthSun(w_e, h);
w_e = ye;
i = i+1;
end
i=1;
while w_m(1) <= sim_time
trajectory_moon2earth(:, i) = w_m;
ym = RungeKutta_MoonEarth(w_m, h);
w_m = ym;
i = i+1;
end
trajectory_moon2sun(1,:) = trajectory_moon2earth(4,:) +
trajectory_earth2sun(4,:); trajectory_moon2sun(2,:) =
trajectory_moon2earth(5,:) + trajectory_earth2sun(5,:); size_enlge = 40; size_up = 15;
trajectory_moon2earth_up(1,:) = size_up* (trajectory_moon2sun(1,:)-trajectory_earth2sun(4,:)); trajectory_moon2earth_up(2,:) = size_up* (trajectory_moon2sun(2,:)-trajectory_earth2sun(5,:));
trajectory_moon2sun_up(1,:) =
(trajectory_earth2sun(4,:)+trajectory_moon2earth_up(1,:));
trajectory_moon2sun_up(2,:) =
(trajectory_earth2sun(5,:)+trajectory_moon2earth_up(2,:)); figure;
for j=1:i-1
draw_circle(0, 0, Rs*5);hold on;
title('地月日二维真实轨迹');
xlabel('x/m');ylabel('y/m');zlabel('z/m');grid on;
plot(trajectory_earth2sun(4, :), trajectory_earth2sun(5, :), '-r');
% plot(trajectory_moon2sun_up(1, :), trajectory_moon2sun_up(2, :), '-b');
axis equal;
plot(trajectory_earth2sun(4, j), trajectory_earth2sun(5, j), 'r*');
plot(trajectory_moon2sun_up(1, j), trajectory_moon2sun_up(2, j),
'k*');
pause(0.001);hold off;
end。