MATLAB模拟地球月球公转自转代码

以下是MATLAB模拟地球月球公转自转的代码示例：

% 定义参数
G = 6.67428e-11;    % 万有引力常数
M_earth = 5.9722e24;    % 地球质量
M_moon = 7.342e22;    % 月球质量
R_earth = 6.371e6;    % 地球半径
R_moon = 1.737e6;    % 月球半径
D = 3.844e8;    % 地球和月球的距离
w_earth = 2*pi/(24*60*60);    % 地球自转角速度
w_moon = 2*pi/(27.3*24*60*60);    % 月球自转角速度

% 初始化状态
x_earth = [0; 0; 0];    % 地球位置
v_earth = [0; 0; 0];    % 地球速度
x_moon = [D; 0; 0];    % 月球位置
v_moon = [0; sqrt(G*M_earth/D); 0];    % 月球速度
theta_earth = 0;    % 地球自转角度
theta_moon = 0;    % 月球自转角度

% 模拟循环
t = 0;
dt = 60;
while t < 365*24*60*60
    % 计算地球和月球之间的引力
    r = x_moon - x_earth;
    F = G*M_earth*M_moon/norm(r)^2;
    a_earth = F*r/norm(r)/M_earth;
    a_moon = -F*r/norm(r)/M_moon;
    
    % 更新地球和月球的位置和速度
    v_earth = v_earth + a_earth*dt;
    x_earth = x_earth + v_earth*dt;
    v_moon = v_moon + a_moon*dt;
    x_moon = x_moon + v_moon*dt;
    
    % 更新地球和月球的自转角度
    theta_earth = theta_earth + w_earth*dt;
    theta_moon = theta_moon + w_moon*dt;
    
    % 绘制地球和月球的位置
    plot3(x_earth(1), x_earth(2), x_earth(3), 'bo');
    hold on;
    plot3(x_moon(1), x_moon(2), x_moon(3), 'ro');
    hold off;
    axis([-D D -D D -D D]);
    xlabel('x');
    ylabel('y');
    zlabel('z');
    title(sprintf('t=%.1f days', t/(24*60*60)));
    drawnow;
    
    t = t + dt;
end

注：以上代码仅为示例，仅考虑地球和月球之间的引力，未考虑其他因素，例如地球和月球的形状、地球的自转对大气运动的影响等。