simulate the solar system in matlab

main.m
% Define initial conditions
G = 6.67430e-11; % gravitational constant
Ms = 1.989e30; % mass of the Sun
Me = 5.97e24; % mass of the Earth
Mm = 7.34e22; % mass of the Moon
rse = 147.09e9; % distance between Sun and Earth
rem = 384.4e6; % distance between Earth and Moon
vse = 30.29e3; % velocity of Earth
vem = 1.02e3; % velocity of Moon
x = [0 rse+rse*cosd(23.5) rem+rse+rse*cosd(23.5)]'; % initial x positions
y = [0 0 rse*sind(23.5)+rem+sind(23.5)*rse]'; % initial y positions
vx = [0 vse*cosd(23.5) vem+vse*cosd(23.5)]'; % initial x velocities
vy = [0 vse*sind(23.5) vse*sind(23.5)+vem]'; % initial y velocities
m = [Ms Me Mm]'; % masses

% Set up simulation parameters
dt = 100; % time step in seconds
tmax = 365*24*60*60; % total simulation time in seconds
t = 0:dt:tmax;
n = length(m);

% Initialize position and velocity matrices
r = [x y];
v = [vx vy];

% Simulate motion using the Euler method
for i = 1:length(t)-1
    % Calculate gravitational forces
    F = zeros(n,2);
    for j = 1:n
        for k = 1:n
            if j ~= k
                rjk = r(k,:) - r(j,:);
                F(j,:) = F(j,:) + G*m(j)*m(k)*rjk/norm(rjk)^3;
            end
        end
    end
    % Update positions and velocities
    r = r + v*dt + 0.5*dt^2*F./m;
    v = v + dt*F./m;
end

% Plot results
figure
hold on
plot(r(:,1),r(:,2),'linewidth',2) % plot trajectory of Sun-Earth-Moon system
plot(x(1),y(1),'r.','markersize',40) % plot Sun
plot(x(2),y(2),'b.','markersize',20) % plot Earth
plot(x(3),y(3),'k.','markersize',10) % plot Moon
axis equal
xlabel('x (m)')
ylabel('y (m)')
title('Solar System Simulation')
legend('Sun-Earth-Moon trajectory','Sun','Earth','Moon')
1662 chars
55 lines