%%
close all, clear, clc

% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -
% Vers.3.0 23.12.2013
%% Parameter und Konstanten
G = 6.67384e-11;      % Gravitationskonstante
m1 = 5.977e24;        % Masse des Erde in kg
m2 = 1.983e30;        % Masse der Sonne in kg
AE = 1.496e11*0.999;  % Mittlere Entfernung Erde-Sonne in m = 1 astronomische Einheit (1AE)
%% Bahnparameter festlegen
epsilon =0.0167;%0.991;%0.5; % Exzentrizit?t (0 Kreis stets <1; Erde: 0.0167)

% wp, isnumeric(wp)
% handles
char('Merkur  1','Venus   2','Erde    3','Mars    4',  'Jupiter 5',...
     'Saturn  6', 'Uranus  7','Neptun  8','Pluto   9','Eig.Werte 10')   
wp = input('Welcher Himmelskoerper soll berechnet werden? Bitte Nr. aus Liste eingeben:   ')
switch wp
    case 1  % Merkur
        m1 = 0.037*m1;
        AE = 57.8e09;
        epsilon = 0.2056;
    case 2   % Venus
        m1 = 0.826*m1;
        AE = 108.1e09;
        epsilon = 0.0068;
    case 3 % Erde
        % Defaultwerte 'oben'
    case 4    % Mars
        m1 = 0.108*m1;
        AE = 227.7e09;
        epsilon = 0.0934;
    case 5    % Jupiter
        m1 = 318.36*m1;
        AE = 777.8e09;
        epsilon = 0.0484;
    case 6    % Saturn
        m1 = 95.22*m1;
        AE = 1425.6e09;
        epsilon = 0.0558;
    case 7    % Uranus
        m1 = 14.58*m1;
        AE = 2868.1e09;
        epsilon = 0.0471;
    case 8    % Neptun
        m1 = 17.27*m1;
        AE = 4494.1e09;
        epsilon = 0.0086;
    case 9    % Pluto
        m1 = 0.0021*m1;
        AE = 5906e09;
        epsilon = 0.2468;
    case 10    % eigen
        m1 = input('Masse Planet m1 ist Erdmasse; Bsp: 0.8*m1   ')
        m2 = input('Masse Sonne m2 ist Sonenmasse; Bsp: 2*m2    ')
        AE = input('Abstand des Himmelsk??rpers von der Sonne; Erde-Sonne AE; Bsp: 0.5*AE  ')
        epsilon = input('Exzentrizit??t  Werte zwischen 0 und <1:  ')  
end

Vorf = G*m1*m2;
mu = m1*m2/(m1+m2); % Reduzierte Masse
Periode =  2*pi*sqrt(1/(G*m2))*AE^(3/2)/3600/24, 'Tage';
%% Startwerte
%r0 = (1+epsilon)*AE;   % Anfangswert: Aphel
r0 = (1-epsilon)*AE;   % Anfangswert: Perihel
theta0 = 0;      
pr0=0;
%
% E0 = -Vorf/(2*r0) * (1+epsilon); % Energie - Aphel
E0 = -Vorf/(2*r0) * (1-epsilon); % Energie - Perihel
pt0 = sqrt((E0*r0+Vorf)*2*mu*r0);
%% Differentialgleichung l??sen
y0 = [r0;theta0;pr0;pt0];
n = input('Berechnungsdauer in Tage:  ')
%n = 365.75;%*10;%*2; % Anzahl der Tage
tmax=24*3600*n;
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',3600*24);
tic
[t,y] = ode45(@(t,y) keplerDGL(t,y,mu,Vorf),[0,tmax],y0,options);
toc
%% Visualisieren
x1=y(:,1).*sin(y(:,2));
x2=y(:,1).*cos(y(:,2));
plot(x1,x2),hold on, plot(x1(1),x2(1),'pr'),plot(0,0,'yp'),axis equal, shg
