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

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

《太阳地球月亮运动轨迹MATLAB仿真程序.docx》由会员分享，可在线阅读，更多相关《太阳地球月亮运动轨迹MATLAB仿真程序.docx（10页珍藏版）》请在冰豆网上搜索。

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

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

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

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

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

程序从第二页开始。

M文件1:

draw_ball.m

%画三维球体的函数

%（x0,y0,z0）为球心

%r为球半径

functiondraw_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为圆半径

functiondraw_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法解地日微分方程

functionyh=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;

%地日微分方程组

functiondy=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法解地月微分方程

functionyh=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;

%地月微分方程组

functiondy=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;clearall;closeall;

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=[0029535.61.52171522e110];%地球相对太阳的初始位置w_m=[w_e

（1）0990.32405500e30];%月球相对地球的初始位置

i=1;%用于记录解算的步数

whilew_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;%用于记录解算的步数

whilew_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/画出地日三维动态轨迹

forj=1:

i-1

draw_ball（0,0,0,Rs*size_enlge）;holdon;%画太阳的三维模型

title（'地日-动态轨迹'）;xlabel（'x/m'）;ylabel（'y/m'）;zlabel（'z/m'）;gridon;

plot（trajectory_earth2sun（4,:

）,trajectory_earth2sun（5,:

）,'-r'）;%画地球相对于太阳的轨迹线

axisequal;

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）;holdoff;

end

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

forj=1:

i-1

draw_ball（0,0,0,Re*size_enlge/4）;holdon;%画地球的三维模型

title（'地月-动态轨迹'）;xlabel（'x/m'）;ylabel（'y/m'）;zlabel（'z/m'）;gridon;

plot（trajectory_moon2earth（4,:

）,trajectory_moon2earth（5,:

）,'-r'）;%画月球相对于地球的轨迹线

axisequal;

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）;holdoff;

end

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

forj=1:

i-1

draw_ball（0,0,0,Rs*size_enlge）;holdon;%画太阳的三维模型

title（'地月日-动态轨迹'）;xlabel（'x/m'）;ylabel（'y/m'）;zlabel（'z/m'）;gridon;

plot（trajectory_earth2sun（4,:

）,trajectory_earth2sun（5,:

）,'-r'）;%画地球相对于太阳的轨

迹线

plot（trajectory_moon2sun（1,:

）,trajectory_moon2sun（2,:

）,'-b'）;%画月球相对于太阳的轨

迹线

axisequal;

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）;holdoff;

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/画出放大后的地月日三维动态轨迹

forj=1:

i-1

draw_ball（0,0,0,Rs*size_enlge）;holdon;%画太阳的三维模型

title（'地月日-动态轨迹（放大）'）;xlabel（'x/m'）;ylabel（'y/m'）;zlabel（'z/m'）;gridon;

plot（trajectory_earth2sun（4,:

）,trajectory_earth2sun（5,:

）,'-r'）;%画地球相对于太阳的轨

迹线

plot（trajectory_moon2sun_up（1,:

）,trajectory_moon2sun_up（2,:

）,'-b'）;%画月球相对于太

阳的轨迹线

axisequal;

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）;holdoff;

end

figure;%/05/画出地月日真实的二维动态轨迹

forj=1:

i-1

draw_circle（0,0,Rs*5）;holdon;%画太阳的二维模型

title（'地月日二维真实轨迹'）;xlabel（'x/m'）;ylabel（'y/m'）;zlabel（'z/m'）;gridon;

plot（trajectory_earth2sun（4,:

）,trajectory_earth2sun（5,:

）,'-r'）;%画地球相对于太阳的轨迹线

plot（trajectory_moon2sun（1,:

）,trajectory_moon2sun（2,:

）,'-b'）;%画月球相对于太阳的轨迹线

axisequal;

plot（trajectory_earth2sun（4,j）,trajectory_earth2sun（5,j）,'r*'）;%轨迹点上的地球点

plot（trajectory_moon2sun（1,j）,trajectory_moon2sun（2,j）,'k*'）;%轨迹点上的月球点

pause（0.001）;holdoff;

end

figure;%/06/画出地月日放大的二维动态轨迹

forj=1:

i-1

draw_circle（0,0,Rs*5）;holdon;%画太阳的二维模型

title（'地月日二维放大轨迹'）;xlabel（'x/m'）;ylabel（'y/m'）;zlabel（'z/m'）;gridon;

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）;holdoff;

end

M文件6:

sun_earth_moon_2Dsingle.m%地月日全运动动态仿真程序

clc;clearall;closeall;

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=[0029535.61.52171522e110];

w_m=[w_e

（1）0990.32405500e30];

whilew_e

（1）<=sim_time

trajectory_earth2sun（:

i）=w_e;

ye=RungeKutta_EarthSun（w_e,h）;

w_e=ye;

i=i+1;

end

i=1;

whilew_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;

forj=1:

i-1

draw_circle（0,0,Rs*5）;holdon;

title（'地月日二维真实轨迹'）;xlabel（'x/m'）;ylabel（'y/m'）;zlabel（'z/m'）;gridon;

plot（trajectory_earth2sun（4,:

）,trajectory_earth2sun（5,:

）,'-r'）;%plot（trajectory_moon2sun_up（1,:

）,trajectory_moon2sun_up（2,:

）,'-b'）;

axisequal;

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）;holdoff;

end