太阳地球月亮运动轨迹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