restart; with(plots): Su2 := K=tan(1/2*psi), psi=146/180*Pi, phi=108/180*Pi, Xv=2, Yv=-0.18, s=[0.1, 0.3, 1.085], S=1.3, h=0.4; HP := x^2/a^2-y^2/b^2=1: HP := simplify(map(u->(u-1)*b^2,HP)): HP := simplify(subs(b=K*a,HP)): HPR := subs([x=x*cos(phi)+y*sin(phi), y=-x*sin(phi)+y*cos(phi)],HP): HPRS := subs(x=x-Xv,y=y-Yv,HPR); assign(Su2); F := expand(evalf(lhs(HPRS))); x1 := Xv-1; x2 := Xv+0.15; y1 := Yv-0.02; y2 := 0; AP[12] := implicitplot([seq(F, a=1/100*[$1..10])], x=x1..x2, y=y1..y2, labels=["x [m]","y [m]"], grid=[60, 60], thickness=[4, seq(1, i=2..10)]): %; read "drive.sav": a := 0.01; dT := a/v; N := ceil(Tf/dT*2); T := Tf/N*[$0..N]: N := nops(T); Gx := map(u->evalf(subs(t=u, Gx)), T): Gy := map(u->evalf(subs(t=u, Gy)), T): Gxt := map(u->evalf(subs(t=u, Gxt)),T): Gyt := map(u->evalf(subs(t=u, Gyt)), T): Gxtt := map(u->evalf(subs(t=u, Gxtt)), T): Gytt := map(u->evalf(subs(t=u, Gytt)), T): Fx := unapply(diff(F, x), x, y); Fy := unapply(diff(F, y), x, y); Fxx := unapply(diff(F, x, x), x, y); Fyy := unapply(diff(F, y, y), x, y); Fxy := unapply(diff(F, x, y), x, y); read "Functions.sav": Dxy:=(Gx[2]-Gx[1])*2; position := subs(fsolve({(Gx[1]-x)^2+(Gy[1]-y)^2=S^2, F}, {x, y}, {x=x1..x2, y=y1..y2}), [x, y]); Hx[1] := position[1]: Hy[1] := position[2]: for i from 2 to N do; position := subs(fsolve({(Gx[i]-x)^2+(Gy[i]-y)^2=S^2, F}, {x, y}, {x=Hx[i-1]-Dxy..Hx[i-1]+Dxy, y=Hy[i-1]-Dxy..Hy[i-1]+Dxy}), [x, y]); Hx[i] := position[1]: Hy[i] := position[2]: if i mod 25 = 0 then print(i) end if; end do: plot([seq([Hx[i], Hy[i]], i=1..N)]); read "plots.sav": Hxt := [seq(Vx(Gx[i], Gy[i], Gxt[i], Gyt[i], Hx[i], Hy[i]), i=1..N)]: Hyt := [seq(Vy(Gx[i], Gy[i], Gxt[i], Gyt[i], Hx[i], Hy[i]), i=1..N)]: pl1 := plot(zip((u, v)->[u, v], Hxt, Hyt), color=grey, thickness=3): AP[13] := display({PL1, pl1}, labels=["Vx [m/s]", "Vy [m/s]"]): %; Hxtt := [seq(Ax(Gx[i], Gy[i], Gxt[i], Gyt[i], Gxtt[i], Gytt[i], Hx[i], Hy[i], Hxt[i], Hyt[i]), i=1..N)]: Hytt := [seq(Ay(Gx[i], Gy[i], Gxt[i], Gyt[i], Gxtt[i], Gytt[i], Hx[i], Hy[i], Hxt[i], Hyt[i]), i=1..N)]: pl2 := plot(zip((u,v)->[u, v], Hxtt, Hytt), color=grey, thickness=3): AP[14] := display({PL2, pl2}, labels=["Ay [m/s^2]", "Ay [m/s^2]"]): %; Ex := (x(t)-X(t))/S; Ey := (y(t)-Y(t))/S; for i from 1 to 3 do; Px[i] := X(t)+Ex*s[i]+Ey*h; Py[i] := Y(t)+Ey*s[i]-Ex*h; Vx[i] := diff(Px[i], t); Vy[i] := diff(Py[i], t); Ax[i] := diff(Vx[i], t); Ay[i] := diff(Vy[i], t); V[i] := sqrt(Vx[i]^2+Vy[i]^2); A[i] := sqrt(Ax[i]^2+Ay[i]^2); At[i] := diff(V[i],t); An[i] := sqrt(A[i]^2-At[i]^2); end do: i := 'i': SU := [diff(X(t),t,t)=Gxtt[i], diff(X(t),t)=Gxt[i], X(t)=Gx[i], diff(Y(t),t,t)=Gytt[i], diff(Y(t),t)=Gyt[i], Y(t)=Gy[i], diff(x(t),t,t)=Hxtt[i], diff(x(t),t)=Hxt[i], x(t)=Hx[i], diff(y(t),t,t)=Hytt[i], diff(y(t),t)=Hyt[i], y(t)=Hy[i]]: AP[15] := plot([seq([seq(subs(SU, [T[i], At[j]]), i=1..N)], j=1..3)], color=black, thickness=[3, 2, 1], labels=["t [s]", "At [m/s^2]"]): %; AP[16] := plot([seq([seq(subs(SU, [T[i], An[j]]), i=1..N)], j=1..3)], color=black, thickness=[3, 2, 1], labels=["t [s]", "An [m/s^2]"]): %; for j from 12 to 16 do; plotsetup(ps, plotoutput=cat(ap, j, `.ps`), plotoptions=`landscape, noborder`); AP[j]; end do: >