# Kapitel 3_a: Der schiefe Wurf ohne Luftwiderstand
# © 2015 Friedrich U. Mathiak,
# mathiak@mechanik-info.de
# 
> restart: with(LinearAlgebra): with(plottools): with(plots):
# In den Impulssatz gehen als äußere Kraft die Gewichtskraft G, die Geschwindigkeit v und die Anfangsgeschwindigkeit v_0 ein.
> G:=Vector([0,-m*g]);
> v:=Vector([diff(x1(t),t),diff(x2(t),t)]);
> v_0:=Vector([v0*cos(alpha0),v0*sin(alpha0)]); 
# Wir notieren den Impulssatz
> IS:=map(int,G,t=0..t) - m*v + m*v_0;
# Damit liegt ein entkoppeltes Differenzialgleichungssystem für die beiden unbekannten Funktionen x1(t) und x2(t) vor
> sys_ode:= IS[1], IS[2];
# Die allgemeine Lösung ist:
> dsolve([sys_ode]);
# Wir lösen nun das Differenzialgleichungssystem unter den vorgegebenen  Anfangswerten
> AW:=x1(0)=0, x2(0)=0;
> loe:=dsolve({sys_ode,AW});assign(%);
# Die Geschwindigkeit ist dann
> v:=map(eval,v);
# und für den Betrag der Geschwindigkeit errechnen wir
> vNt:=unapply(VectorNorm(v,2,conjugate = false),v0,alpha0,t);
# oder durch Elimination der Zeit t
> vNx1:=unapply(vNt(v0,alpha0,x1/v0/cos(alpha0)),v0,alpha0,x1);
# und für die Gleichung der Wurfparabel erhalten wir
> x2:=unapply(subs(t=x1/v0/cos(alpha0),x2(t)),x1,v0,alpha0);
# Für die Einhüllende der Wurfparabeln gilt (s.h. Buch TM III)
> EH:=v0^2/2/g*(1-(g*x1/v0^2)^2);
# Die Wurfhöhe  H erhalten wir aus der Forderung, dass die Neigung der Tangente an die Wurfparabel x2(x1) im Scheitelpunkt verschwinden muss. Wir bilden zunächst die Ableitung 
# x2'  = dx2/dx1 und erhalten
> x2str:=diff(x2(x1,v0,alpha0),x1);
# Lösen wir diese Beziehung nach x1 = x1H auf, dann folgt 
> x1H:=unapply(solve(x2str,x1),v0,alpha0);
# und für die Wurfhöhe errechnen wir
> H:=unapply(x2(x1H(v0,alpha0),v0,alpha0),v0,alpha0);
# Die größte Wurfhöhe erreichen wir mit dem senkrechten Wurf (± = À/2) zu
> Hmax:=H(v0,Pi/2);
# Wir berechnen die Wurfweite und setzen dazu  x1 = x1w in die Gleichung für x2 ein:
> gl:=-h=x2(x1w, v0, alpha0);
> loe:=solve(gl,x1w);
# Wir wählen für x1w die positive Lösung (Wurf nach rechts)
> x1w:=unapply(loe[1],v0,alpha0,h);
# Dazu gehört die Wurfzeit
> tw:=unapply(x1w(v0,alpha0,h)/v0/cos(alpha0),v0,alpha0,h);
# Speziell für den Fall h = 0 erhalten wir 
> x1w0:=x1w(v0,alpha0,0); tw0:=tw(v0,alpha0,0);
# Am Ende der Wurfbahn hat die Punktmasse die Geschwindigkeit
> vN_E:= simplify(vNx1(v0,alpha0,x1w(v0,alpha0,h)));
# Zur Berechnung der größten Wurfweite benötigen wir die Ableitung von x1w nach ±:
> x1wstr:=diff(x1w(v0,alpha0,h),alpha0);
# Die größte Wurfweite ergibt sich für den Abwurfwinkel
> alpha0max:=solve(x1wstr,alpha0);
# Beispiel 3-2:
# Im Folgenden wird mit den Werten des Beispiels eine nummerische Rechnung durchgeführt.
> v0:=10.; h:=3.; g:=9.81; x0:=10.;
> plot([x0,x1w(v0,alpha0,h)],alpha0=-Pi/2..Pi/2,axes = boxed,gridlines=true,color=[blue,black],title = "\nWurfweite in Abhängigkeit vom Abwurfwinkel für h = 3 m\n",titlefont = ["ARIAL", 15], labels = ["Abwurfwinkel [rad]", "Wurfweite [m]"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8]);
# Der obigen Abbildung entnehmen wir, dass dieselbe Wurfweite (von beispielsweise x0 = 10 m) durch zwei verschiedene Abwurfwinkel erreicht wird:
> alpha01:=fsolve(x0-x1w(v0,alpha0,h),alpha0=0..Pi/4);alpha01deg:=evalf(convert(alpha01,degrees));
> alpha02:=fsolve(x0-x1w(v0,alpha0,h),alpha0=Pi/4..Pi/2);alpha02deg:=evalf(convert(alpha02,degrees));
# Die zugehörigen Wurfzeiten sind
> t1w:=tw(v0,alpha01,h);
> t2w:=tw(v0,alpha02,h);
# Wir prüfen das Ergebnis, indem wir für beide Abwurfwinkel die Wurfweite berechnen.
> l1:=x1w(v0,alpha01,h);l2:=x1w(v0,alpha02,h);
# und geben die Wurfbahnen grafisch aus
> p1:=plot([x2(x1,v0,alpha01),x2(x1,v0,alpha02)],x1=0..l1,axes=boxed,gridlines=true,color=[blue,red]):
> p2:=line([0, 0], [0, -h],color = black, thickness=5):
> display(p1,p2,scaling=constrained,title="\nWurfbahnen für den Flach- u. Steilwurf", titlefont = ["ARIAL", 15], labels = ["x1 [m]", "x2 [m]"],labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10]);
# Beim  Steilwurf ergibt sich eine Wurfhöhe von 
> H(v0,alpha02);
# Die größte Wurfweite L wird erreicht für den Abwurfwinkel
> alpha03:=alpha0max[1];evalf(convert(alpha03,degrees));
> L:=x1w(v0,alpha03,h);
# Die zugehörige Wurfzeit ist
> t3w:=tw(v0,alpha03,h);
# Für den senkrechten Wurf erhalten wir die Wurfzeit
> tsw:=tw(v0,Pi/2,h);
> plot([x1w(v0,alpha0,0),x1w(v0,alpha0,1),x1w(v0,alpha0,h)],alpha0=-Pi/2..Pi/2,gridlines=true,axes=boxed,color=[blue,red,green],title = "\nWurfweiten in Abhängigkeit vom Abwurfwinkel\n", legend = ["h=0 m", "h=1 m","h=3 m"], titlefont = ["ARIAL", 15], labels = ["Abwurfwinkel [rad]", "Wurfweite [m]"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom]);
# Die zugehörigen Wurfzeiten in Abhängigkeit vom Abwurfwinkel ± können der folgenden Grafik entnommen werden
> plot([tw(v0,alpha0,0),tw(v0,alpha0,1),tw(v0,alpha0,h)],alpha0=-Pi/2..Pi/2,y=0..tsw,gridlines=true,axes=boxed,color=[blue,red,green],title = "\nWurfzeiten in Abhängigkeit vom Abwurfwinkel\n", legend = ["h=0 m", "h=1 m","h=3 m"], titlefont = ["ARIAL", 15], labels = ["Abwurfwinkel [rad]", "Wurfzeit [s]"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom]);
# Wir stellen die Wurfparabeln mit ihrer Einhüllenden grafisch dar
> plot([x2(x1,v0,alpha01),x2(x1,v0,alpha02),x2(x1,v0,alpha03),EH],x1=0..L,y=-h..Hmax,gridlines=true,axes=boxed,color=[blue,red,green,black],title = "\nWurfbahnen in Abhängigkeit vom Abwurfwinkel\n", legend = ["Flachwurf", "Steilwurf","größte Wurfweite","Einhüllende"], titlefont = ["ARIAL", 15], labels = ["x1 [m]", "x2 [m]"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8], legendstyle = [font = ["HELVETICA", 9], location = bottom],linestyle=[solid,solid,solid,dash]);
# Die Koordinaten des Scheitelpunktes für den Flachwurf sind
> x1S:=x1H(v0,alpha01);x2S:=H(v0,alpha01);
# Die Koordinaten des Scheitelpunktes für den Steilwurfwurf sind
> x1S:=x1H(v0,alpha02);x2S:=H(v0,alpha02);
# Die Koordinaten des Scheitelpunktes für die größte Wurfweite sind
> x1S:=x1H(v0,alpha03);x2S:=H(v0,alpha03);
# Geschwindigkeiten als Funktion der Zeit:
> v1:=plot(vNt(v0,alpha01,t),t=0..t1w,axes=boxed,color= blue):
> v2:=plot(vNt(v0,alpha02,t),t=0..t2w,axes=boxed,color=  red):
> v3:=plot(vNt(v0,alpha03,t),t=0..t3w,axes=boxed,color=green):
> display(v1,v2,v3,gridlines=true,title = "\nGeschwindigkeiten v(t)\n", titlefont = ["ARIAL", 15], labels = ["Zeit t in [s]", "Geschwindigkeit v [m/s]"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8]);
# Geschwindigkeiten als Funktion des Weges:
> v1:=plot(vNx1(v0,alpha01,x1),x1=0..x1w(v0, alpha01, h),color= blue):
> v2:=plot(vNx1(v0,alpha02,x1),x1=0..x1w(v0, alpha02, h),color=  red):
> v3:=plot(vNx1(v0,alpha03,x1),x1=0..x1w(v0, alpha03, h),color=green):
> display(v1,v2,v3,gridlines=true,title = "\nGeschwindigkeiten v(x1) \n", titlefont = ["ARIAL", 15], labels = ["x1 in [m]", "Geschwindigkeit v [m/s]"], labeldirections = ["horizontal", "vertical"], labelfont = ["HELVETICA", 10], axesfont = ["HELVETICA", "ROMAN", 8],axes=boxed);
# Am Ende des Wurfes beträgt die Geschwindigkeit und zwar unabhängig vom Abwurfwinkel ±, 
> vN_E;
# Es folgt abschließend eine Animation der Wurfparabeln.  Der Senkrechtwurf mit alpha02 = Pi/2  ist bei der folgenden Animation nicht zulässig.
> nframes:=50:
> pl1:=animatecurve(x2(x1,v0,alpha01),x1=0..l1,frames=nframes,color= blue):
> pl2:=animatecurve(x2(x1,v0,alpha02),x1=0..l1,frames=nframes,color=  red):
> pl3:=animatecurve(x2(x1,v0,alpha03),x1=0.. L,frames=nframes,color=green):
> display(pl1,pl2,pl3,axes=boxed,gridlines=true,view=-h..Hmax,scaling = constrained,title="\nWurfparabeln verschiedener Abwurfwinkel",titlefont=["ROMAN", 15],labels=["x1 [m]","x2 [m]"],labeldirections=["horizontal", "vertical"], labelfont = ["HELVETICA", 10]);
> 
;
> 
;
