%%
%close all, clear, clc
%figure
tic
% Vers.2.2 28.12.2013
% W. Schweizer
% Simulation physikalischer Systeme - Computational Physics mit MATLAB -

%% Parameter und Konstanten
G = 6.67384e-11; % Gravitationskonstante
m1 = 6e1;        % Masse des Satelliten in kg
m2 = 5.977e24;   % Masse des Erde in kg
%AE = 1.496e11;  % Mittlere Entfernung Erde-Sonne in m = 1 astronomische Einheit (1AE)
ae = 6.378140e06; % mittlerer Erdradius in m
J2 = 0.001082626; % Multipolmomente der Erde 
J3 = -0.254e-05;
J4 = -0.161e-05;
Vorf = G*m1*m2;
mu = m1*m2/(m1+m2); % Reduzierte Masse
%% Bahnparameter festlegen
h = input('Satellitenbahn: Hoehe ueber der Erdoberflaeche '); 
%rm = ae+4e05;                   % Abstand in m
rm = ae +h*1000;
epsilon = input('0 <= Exzentrizitaet der Bahn < 1 ');
%epsilon =0.95;%0.991;%0.5; % Exzentrizit?t (0 Kreis stets <1; Erde: 0.0167)
%% Umlaufzeit
Tn = sqrt(4*pi^2/(G*m2)*rm.^3)/3600; % in Stunden  --> n
Ausgabe = ['Wieviele Orbits? Die Umlaufzeit betraegt ',num2str(Tn), ' Stunden '];
n = input(Ausgabe);
%% Startwerte
r0 = (1+epsilon)*rm;   % Anfangswert: Aphel
theta0 = input('Winkel gegen die Polachse in Grad ');
theta0 = theta0*pi/180;;%pi/6;%0;%pi/2;%pi/5;%0;      
pr0=0;
%
E0 = -Vorf/(2*r0) * (1+epsilon); % Energie
pt0 = sqrt((E0*r0+Vorf)*2*mu*r0);
%% Differentialgleichung l??sen
y0 = [r0;theta0;pr0;pt0];
tmax=3600*Tn*n;
options = odeset('RelTol',1e-8,'AbsTol',1e-8,'MaxStep',3600*24);
tic
[t,y] = ode45(@(t,y) SatErdeDGL(t,y,mu,Vorf,ae,J2,J3,J4),[0,tmax],y0,options);
toc
%% Visualisieren
x1=y(:,1).*sin(y(:,2));
x2=y(:,1).*cos(y(:,2));
plot(x1,x2,'k'),hold on, plot(x1(1),x2(1),'pr'),plot(0,0,'yp'),axis equal, shg
hold on
toc