######################
##
##  Aufgabe 1 a)
##
#######################

# Anfangswerte und Konstanten
hP = 215;
hA = 939;
omega0 = deg2rad(340);
i0 = deg2rad(65.1);
w0 = deg2rad(58);
theta0 = deg2rad(332);
m = 100;
d = 1;
cD = 2.2;
wE = 72.9211E-6;

rE = 6378;
a0 = (hP + hA + 2*rE)/2;
e0 = 1 - (hP + rE)/a0;
p0 = a0*(1-e0^2);
r0 = p0/(1 + e0*cos(theta0));
h0 = r0 - rE;
j2 = 1082.6E-6;
mu = 398600;
n0 = sqrt(mu/(a0^3));
E0 = 2*atan2(sqrt(1-e0)*tan(theta0/2), sqrt(1+e0));
D = cD*pi*(d^2)/m;
rP = p0/(1 + e0*cos(0));

a = a0;
e = e0;
p = p0;
n = n0;
p = p0;

r = zeros(101,1);
dichte = zeros(101,1);
valueToIntegrate1 = zeros(101,1);
valueToIntegrate2 = zeros(101,1);


# Um die Änderung für a und e zu berechnen betrachte zuerst den Integralteil der Gleichungen.
# Simuliere eine Orbitperiode um die numerischen Werte innerhalb des Integrals zu bestimmen und löse dieses
# anschließend mithilfe die trapz funktion für einen Orbitdurclauf und berechne so a und e
# Update anschließend alle Werte mit den aktualisierten Werten von a und e und gehe analog für den
# nächsten Orbit vor bis die Perigäumshöhe auf 100 km abgefallen ist.

# Teil a) funktioniert leider nicht da die Änderung von a und e so klein berechnet wurden das sie keine
# Auswirkung auf die Rechnung haben. Programm läuft deswegen in Endlosschleife

while((rP - rE) > 100)
i = 1;
for k = 0:pi/50:2*pi
  theta = theta0 + k;
  r(i) = p/(1 + e*cos(theta));
  dichte(i) = rho(r(i) - rE);

  valueToIntegrate1(i) = dichte(i)*((1 + (e^2) + 2*e*cos(theta))^(3/2))/((1 + e*cos(theta))^(2));
  valueToIntegrate2(i) = dichte(i)*(((1 + (e^2) + 2*e*cos(theta))^(3/2))/((1 + e*cos(theta))^(2)))*(cos(theta) + e);

  i = i + 1;
endfor

adot = -D*(a^2)*(n/(2*pi))*trapz(0:pi/50:2*pi, valueToIntegrate1)
edot = -D*a*(1-(e^2))*(n/(2*pi))*trapz(0:pi/50:2*pi, valueToIntegrate2)

# scaling factor
o = 1;
a = a + o*adot
e = e + o*edot;

p = a*(1-e^2);
rP = p/(1 + e*cos(0));
n = sqrt(mu/(a^3));

end



# Plot der Dichte Funktion

##for i = 95:1100
##  s(i) = rho(i);
##endfor

##k = 600
##
##figure
##plot(k:1100, s(k:1100))
###ylim([0 1E-6])
##title('Plot Dichte')
###xlabel('Zeit [s]'); ylabel('v [km/s]')
##grid on



######################
##
##  Teil b)
##
#######################

# da teil a) nicht gelöst wurde wähle t und rechne damit
t= 24*60*60*7;

k = (3/2)*((j2*(rE^2))/(p0^2))*n0;
omega = omega0 - k*cos(i0)*t;
rad2deg(omega);

w = w0 + k*(2 - (5/2)*(sin(i0)^2))*t;
rad2deg(w);