%Pflichtuebung 11
%Daniel Jung 1411014

%Alle Variablen aus dem Speicher löschen
clear all;
clear fig;

%Parameter der Differentialgleichung festlegen:
global m g c p d A
m=0.145;
g=9.81;
c=0.5;
p=1.2;
d=0.5;
A=pi/(4*d^2);

%Anfangs- und Endzeitpunkt festlegen:
t0=0;               %Anfangszeitpunkt
tE=0.3;               %Endzeitpunkt
N=10000;            %Anzahl der gespeicherten Datenpunkte
deltatau=(tE-t0)/(N-1);
tauspan=t0:deltatau:tE;

%Anfangsparameter abfragen
alfarad=input('Abwurfwinkel in Grad: ');
v0=input('Betrag der Anfangsgeschwindigkeit in m/s: ');

alfa=alfarad*(pi/180);  %Grad in rad rechnen

x(1,1)=0;               %Wert in x-Richtung
x(2,1)=0;               %Wert in y-Richtung
x(3,1)=v0*cos(alfa);    %Geschwindigkeit in x-Richtung
x(4,1)=v0*sin(alfa);    %Geschwindigkeit in y-Richtung

xstart=[x(1,1); x(2,1); x(3,1); x(4,1)]; %Vektor bilden

[t,X]=ode45(@wurf,tauspan,xstart);       %Aufruf der Funktion

ta=t0;
for n=1:N
    sx(n)=x(3,1)*ta(n);
    sy(n)=-1/2*g*ta(n)^2+x(4,1)*ta(n);
    ta(n+1)=ta(n)+deltatau;
end

%Ergebnisse plotten
figure(1)
plot(X(:,1),X(:,2))
hold on;
plot(sx(1,:),sy(1,:))
xlabel('Weite in m');
ylabel('Höhe in m');

function f=wurf(t,x)            %Funktion
    global m g c p A
    f(1,1)=x(3);
    f(2,1)=x(4);
    f(3,1)=-(c*p*A*sqrt(x(3)^2+x(4)^2)*x(3))/(2*m);
    f(4,1)=-((c*p*A*sqrt(x(3)^2+x(4)^2)*x(4))/(2*m))-g;
end