function three_body_planar () global G m1 m2 m3 G = 6.67384e-11; % Masse Sonne, Erde, Jupiter x 500 [kg] m1 = 2e30; m2 = 6e24; m3 = 500*2e27; % AWB Ortsposition [m] x01 = [0; 0]; x02 = [0; 1.5e11]; x03 = [5*1.5e11; 0]; % AWB Geschwindigkeiten [m/s] v01 = [0; 0]; tmp = 2*pi*x02(2)/(365*24*3600); v02 = [-tmp; 0]; tmp = 2*pi*x03(1)/(12*365*24*3600); v03 = [0; tmp]; % AWB in Vektoren zusammenfassen x0 = [x01(1); x02(1); x03(1); x01(2); x02(2); x03(2); ... v01(1); v02(1); v03(1); v01(2); v02(2); v03(2)]; tend = 20*365*24*3600; %'InitialStep' and/or 'MaxStep' options = odeset ("RelTol", 1e-9, "AbsTol", 1e-9); [t, y] = ode45 (@threebodies, [0, tend], x0', options); h = figure; hax = axes (); hold (hax, "on"); grid (hax, "on"); % Bahnkurven plot (y(:, 1), y(:, 4), 'r-'); plot (y(:, 2), y(:, 5), 'g-'); plot (y(:, 3), y(:, 6), 'b-'); % Plotten der AWB Massen plot (x01(1), x01(2), 'r*'); plot (x02(1), x02(2), 'g*'); plot (x03(1), x03(2), 'b*'); xrange = [min(min(y(:, 1:3))), max(max(y(:, 1:3)))]; yrange = [min(min(y(:, 4:6))), max(max(y(:, 4:6)))]; xlim (xrange); ylim (yrange); %dz = 1/sqrt((xrange(2)-xrange(1))^2+(yrange(2)-yrange(1))^2); %daspect ([1, 1, dz]); axis equal; pause (3); for i = 1 : 4 : numel (t) tact = tend * (i - 1) / numel(t); x1act = interp1 (t, y(:,1), tact); y1act = interp1 (t, y(:,4), tact); x2act = interp1 (t, y(:,2), tact); y2act = interp1 (t, y(:,5), tact); x3act = interp1 (t, y(:,3), tact); y3act = interp1 (t, y(:,6), tact); cla (hax); scatter (hax, x1act, y1act, 30, "filled", "o", "markerfacecolor", "r"); scatter (hax, x2act, y2act, 30, "filled", "o", "markerfacecolor", "g"); scatter (hax, x3act, y3act, 30, "filled", "o", "markerfacecolor", "b"); drawnow; end end % dx = [x1', x2', x3', y1', y2', y3', v1', v2', v3', w1', w2', w3'] % x = [x1, x2, x3, y1, y2, y3, v1, v2, v3, w1, w2, w3] function dx = threebodies (t, x) global G m1 m2 m3 % Abstand dc12 = sqrt ((x(1)-x(2)).^2 + (x(4)-x(5)).^2).^3; dc13 = sqrt ((x(1)-x(3)).^2 + (x(4)-x(6)).^2).^3; dc23 = sqrt ((x(2)-x(3)).^2 + (x(5)-x(6)).^2).^3; % 12 Gleichungen erster Ordnung dx = zeros(12, 1); % x1' = v1; x2' = v2; x3' = v3 dx(1) = x(7); dx(2) = x(8); dx(3) = x(9); % y1' = w1; y2' = w2; y3' = w3 dx(4) = x(10); dx(5) = x(11); dx(6) = x(12); % v1' = G*m2*(x2-x1)/d12^3 + G*m3(x3-x1)/d13^3 dx(7) = G * m2 * (x(2)-x(1))./dc12 + G * m3 * (x(3)-x(1))./dc13; dx(8) = G * m1 * (x(1)-x(2))./dc12 + G * m3 * (x(3)-x(2))./dc23; dx(9) = G * m1 * (x(1)-x(3))./dc13 + G * m2 * (x(2)-x(3))./dc23; % w1' = G*m2*(y2-y1)/d12^3 + G*m3(y3-y1)/d13^3 dx(10) = G * m2 * (x(5)-x(4))./dc12 + G * m3 * (x(6)-x(4))./dc13; dx(11) = G * m1 * (x(4)-x(5))./dc12 + G * m3 * (x(6)-x(5))./dc23; dx(12) = G * m1 * (x(4)-x(6))./dc13 + G * m2 * (x(5)-x(6))./dc23; end