%%  SwingBy Technik am Beispiel Jupiter
% Ausgangsbahn des Satelliten Parabel oder Hyperbel
% close all, 
% figure,clear, clc
% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
tic
%% Parameter und Konstanten
G = 6.6726e-11; % Gravitationskonstante
m1 = 1e6;       % Masse des Satelliten in kg.
% Jupiterdaten
% m2 = 1.193e30;  % Zum Vergleich: Masse der Sonne in kg
% AE = 1.496e11;  % Zum Vergleich: Entfernung Erde-Sonne in m = 1 astronomische Einheit (1AE)
m2 = 5.977e24;   % Masse des Erde in kg
m2 = 318.36*m2;  % Masse des Jupiters in kg
AE = 777.8e09;   % Abstand Jupiter Sonne
epsilon = 0.0484;% wird zu Null gesetzt; Annahme ideale Kreisbahn 
%
Vorf = G*m1*m2;
mu = m1*m2/(m1+m2); % Reduzierte Masse
%% Bahnparameter festlegen
rmin = 1e09; % Mindestabstand von Jupiter in m 
rmin = input('Mindestabstand zu Jupiter in km; Bsp. 1e06  ');
rmin = rmin*1000;
pt0 = sqrt(2*mu*Vorf*rmin);  % Drehimpuls zu E=0
%% Startwerte bestimmen
r0 = 150*rmin;   % Abstand fuer Start der Rechnung
was = input('Parabel ''P'' oder Hyperbel ''H'' ?');
thetastart = input('Winkel gegen die Jupiterbahn 0 ... 360  ');
thetastart = thetastart*pi/180;
if 'P' == was(1) 
disp('Parabel')
Er = Vorf/r0-pt0^2/(2*mu*r0^2);
pr0 = sqrt(2*mu*Er);
C = mu*Vorf/pt0^2;
theta0=acos((1/(r0*C)-1))+thetastart;
else
disp('Nehme Hyperbel')
epsilon = input('Exzentrizit??t > 1  ')
if epsilon < 1
    disp('falscher Wert, setze')
    epsilon = 1.05
end
nullst=[1./(2*mu*rmin^2),-Vorf/rmin,-(epsilon^2-1)*mu*Vorf^2/2];
pt0=sqrt(max(roots(nullst)));
C = mu*Vorf/pt0^2;
E = pt0^2/(2*mu*rmin^2)-Vorf/rmin;
Er = E + Vorf/r0-pt0^2/(2*mu*r0^2);
pr0 = sqrt(2*mu*Er);
theta0 = acos((1/(r0*C)-1)/epsilon) + thetastart;
end

%%
%y0 = [r0;theta0;-pr0;-pt0];
y0 = [r0,theta0,-pr0,-pt0];
%y0 = [y(3500,1),y(3500,2),-y(3500,3),y(3500,4)];
clear y
n = 3000; % Anzahl der Tage
tmax=24*3600*n;
%options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',3600*24);
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',3600*24,'Events', @(t,y) swingByevent(t,y,r0));
[t,y] = ode45(@(t,y) kometenDGL(t,y,mu,Vorf),[0,tmax],y0,options);
%%
x1=y(:,1).*sin(y(:,2));
x2=y(:,1).*cos(y(:,2));
plot(x1,x2),hold on, plot(x1(1),x2(1),'pr'),shg,axis equal,grid on
title('Bahn im ruhenden Jupitersystem')
% E = y(:,3).^2/(2*mu) + y(:,4).^2./(2*mu*y(:,1).^2) - Vorf./y(:,1);
% figure;
% subplot(2,1,1)
% plot(y(:,1),E./V0),ylabel('rel.Energie')
% subplot(2,1,2)
% epsilon = sqrt(1+2*E*pt0^2/(mu*Vorf^2));
% plot(y(:,1),epsilon-1),ylabel('epsilon-1')
toc
% x1pl=x1;  Fuer Plotzwecke -> Buch
% x2pl=x2;