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

太阳、地球、月亮运动轨迹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 111 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;% 解算一步则在