2024年4月15日发(作者：)

太阳地球月亮运动轨迹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_earth2sun(5,:);

= trajectory_moon2earth(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,

Re*size_enlge^2);% 初始位置的地球三维模型

draw_ball(trajectory_earth2sun(4, j), trajectory_earth2sun(5,

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;% 画地球的三维模型

1), 0,

j), 0,

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,

Rm*size_enlge/4);%

初始位置的月球三维模型

draw_ball(trajectory_moon2earth(4,

Rm*size_enlge/4);%

轨迹点上的月球三维模型

pause(0.001);hold off;

end

figure;% /03/画出地月日三维动态轨迹

for j=1:i-1

1), trajectory_moon2earth(5,

j), trajectory_moon2earth(5,

1), 0,

j), 0,

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,

Re*size_enlge/2);% 初

1), trajectory_earth2sun(5, 1), 0,

始位置的地球三维模型

draw_ball(trajectory_moon2sun(1,

Rm*size_enlge/2);%

1), trajectory_moon2sun(2, 1), 0,

初始位置的月球三维模型

draw_ball(trajectory_earth2sun(4,

Re*size_enlge/2);% 轨迹

j), trajectory_earth2sun(5, j), 0,

点上的地球三维模型

draw_ball(trajectory_moon2sun(1,

Rm*size_enlge/2);% 轨

j), trajectory_moon2sun(2, j), 0,

迹点上的月球三维模型

pause(0.001);hold off;

end

% 放大地月之间的相对位置，以便三维视觉显示

size_up = 15;% 放大倍数，用于放大地月之间的相对位置

% trajectory_moon2earth_up存储放大后的月球到地球的位置数据

trajectory_moon2earth_up(1,:)

trajectory_earth2sun(4,:));

= size_up* (trajectory_moon2sun(1,:)-

trajectory_moon2earth_up(2,:)

trajectory_earth2sun(5,:));

= size_up* (trajectory_moon2sun(2,:)-

% 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),

Re*size_enlge*10);%

初始位置的地球三维模型

draw_ball(trajectory_moon2sun_up(1, 1), trajectory_moon2sun_up(2, 1),

Rm*size_enlge*10);% 初始位置的月球三维模型

draw_ball(trajectory_earth2sun(4, j), trajectory_earth2sun(5, j),

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;

0,

0,

0,

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,:)

trajectory_earth2sun(4,:));

= size_up* (trajectory_moon2sun(1,:)-

= size_up* trajectory_moon2earth_up(2,:)

(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