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

1、太阳地球月亮运动轨迹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%画二维圆形

2、% (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;k

3、3 = 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)/r21.5; dy(3) = (-1)*miu_s*w(5)/r21.5; dy(4) = w(2);dy(5) = w(3);M 文件 4:RungeKutta_Mo

4、onEarth.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、5.972e24; miu_e = G * Me; r2 = w(4)”2+w(5)”2;dy(1) = 1;dy(2) = (-1)*miu_e*w(4)/r21.5; dy(3) = (-1)*miu_e*w(5)/r21.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;%地球的质量和半径

6、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_earth

7、2sun 存储地球相对 于太阳的轨迹数据ye = RungeKutta_EarthSun(w_e, h);% 4 阶 Runge-kutta 法解算地日微分方程w_e = ye;%用于Runge-kutta法解算微分方程的初值i = i+1;%解算一步则在步数上+1endi = 1;%用于记录解算的步数while w_m(1) = sim_time%解算月球相对于地球的轨迹trajectory_moon2earth (:, i) = w_m;% trajectory_moon2earth 存储月 球相 对于地球的轨迹数据ym = RungeKutta_MoonEarth(w_m, h);% 4

8、 阶 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;%放大系数，用于放大太阳、地球、月球的半

9、径，以便视觉 观测figure;% /01/画出地日三维动态轨迹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

10、,Re*size_enlge2);%初始位置的地球三维模型draw_ball(trajectory_earth2sun(4, j), trajectory_earth2sun(5, j), 0,Re*size_enlge2);%轨迹点上的地球三维模型pause(0.001);hold off;endfigure;% /02/画出地月三维动态轨迹for j=1:i-1draw_ball(0, 0, 0, Re*size_enlge/4); hold on;% 画地球的三维模型title(地月-动态轨迹)；xlabel(x/m);ylabel(y/m);zlabel(z/m);grid on;pl

11、ot(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);ho

12、ld off;endfigure;% /03/画出地月曰三维动态轨迹for j=1:i-1draw_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);

13、%画月球相对于太阳的轨迹线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/

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

15、 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,

16、:)+trajectory_moon2earth_up(2,:);figure;% /04/画出放大后的地月日三维动态轨迹for j=1:i-1draw_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, :)

17、, 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), traj

18、ectory_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;endfigure;% /05/画出地月日真实的二维动态轨迹draw_circle(0, 0, Rs*5); hold on;% 画太阳的二维模型title(地月日二维真实轨迹)；xlabel(x/m);ylabel(y/m);zlabel

19、(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_

20、moon2sun(2, j), k*);% 轨迹点上的月球点pause(0.001);hold off;endfigure;% /06/画出地月日放大的二维动态轨迹for j=1:i-1draw_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_moon2su

21、n_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;endM文件6:sun_earth_moon_2Dsingle.m %地月曰全运动动态仿真程序 clc;clear all;close

22、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_

23、e = ye;i = i+1;endwhile w_m(1) = sim_timetrajectory_moon2earth(:, i) = w_m;ym = RungeKutta_MoonEarth(w_m, h);w_m = ym;i = i+1;endtrajectory_moon2sun(1,:) = trajectory_moon2earth(4,:) +trajectory_earth2sun(4,:); trajectory_moon2sun(2,:)=trajectory_moon2earth(5,:) + trajectory_earth2sun(5,:); size_enlge = 40; size_up = 15;traje