# The Classical Tractrix / int
restart;
y := -Int(sqrt(a^2-x^2)/x, x)+c;
y := simplify(value(y), symbolic);
c := solve(subs(x=a, y), c);
y := combine(y, ln, symbolic);

# The Classical Tractrix / dsolve
unassign('y');
eq := (y(x)/diff(y(x), x))^2 + y(x)^2 = a^2;
p:=[dsolve(eq,y(x))]: 
p1:=simplify(p,symbolic);
p2:=subs({x=0,y(x)=a},p1);
p3:=map(u->solve(u,_C1),p2);
zip((u,v)->`if`(u=c,subs(_C1=c,v),NULL),p3,p1)[];

# Jogger with Constant Velocity
restart;
Xe := m1 + a*cos(s(t)): Ye := m2 + b*sin(s(t)):
Ve2 := diff(Xe, t)^2 +diff(Ye, t)^2 =V^2;
diff(s(t),t)={solve(Ve2, diff(s(t), t))};
T := 2*Pi:
L := evalf(4*int(sqrt(20^2*sin(f)^2 + 15^2*cos(f)^2), 
     f = 0..Pi/2));
V := evalf(L/T);


# Using a Moving Coordinate System
restart;
S(t) := sqrt((X(t) - x(t))^2 + (Y(t) - y(t))^2):
rx := diff(x(t), t) = W*(X(t) - x(t))/S(t);             
ry := diff(y(t), t) = W*(Y(t) - y(t))/S(t); 
Rx := x(t) = X(t) + rho(t)*cos(phi(t)): 
Ry := y(t) = Y(t) + rho(t)*sin(phi(t)):
Vx := diff(Rx, t);  Vy := diff(Ry, t);
qx := subs(Vx, rx):  qy := subs(Vy, ry):
qx := simplify(subs(Rx, Ry, qx), symbolic);
qy := simplify(subs(Rx, Ry, qy), symbolic);
Dsys := solve({qx, qy}, {diff(rho(t), t), diff(phi(t), t)});
Dradial:=select(has,Dsys,diff(rho(t),t))[];
Daxial:=select(has,Dsys,diff(phi(t),t))[];

#Transformation for Child/Toy
restart; 
Rx := x(t) = X(t) + a*cos(phi(t)):
Ry := y(t) = Y(t) + a*sin(phi(t)):
Vx := diff(Rx, t): Vy := diff(Ry, t):
RW:=W(t)=subs(Vx,Vy,sqrt(diff(x(t),t)^2+diff(y(t),t)^2));
qx := diff(x(t), t) = W(t)/a*(X(t) - x(t));
qy := diff(y(t), t) = W(t)/a*(Y(t) - y(t));
qx:=subs(RW, Vx, Rx, map(u -> u^2, qx));
qy:=subs(RW, Vy, Ry, map(u -> u^2, qy));
Subst := [rhs(Vx) = A, rhs(Vy) = B, phi(t) = F];
q1x := expand(subs(Subst, qx)); 
q1y := expand(subs(Subst, qy));
[map(u->simplify(sqrt(u-cos(F)^2*A^2),symbolic),q1x),
map(u->simplify(sqrt(-u+sin(F)^2*A^2+B^2),symbolic),q1y)];
Q1:=subs(sqrt(1-cos(F)^2)=sin(F),%[1]);
BackSubst:=map(u->rhs(u)=lhs(u),Subst);
Q2 := subs(BackSubst, Q1);
Daxial := diff(phi(t), t) = 
            simplify(solve(Q2, diff(phi(t), t)));

# Examples
Child_Subst := X(t) = R*cos(omega*t), Y(t) = R*sin(omega*t);
Das := subs(Child_Subst, Daxial);
Das:=combine(Das,trig);
Sol1 := dsolve(Das, phi(t))  assuming a>R;
Sol2 := dsolve(Das, phi(t))  assuming R>a;
SuC1:=_C1 = solve(subs(t = 0, phi(0) = alpha, Sol1),_C1);
SuC2:=_C1 = solve(subs(t = 0, phi(0) = alpha, Sol2),_C1);
Sol1:=subs(SuC1,Sol1):
Sol2:=subs(SuC2,Sol2):
Sol10:=eval(subs(alpha=0,Sol1));
Sol20:=eval(subs(alpha=0,Sol2));
limit(rhs(Sol10), a = R),limit(rhs(Sol20), a = R);
Sol0 :=phi(t)= %[1];
Toy := [x(t), y(t)]:
Toy_Plot := subs(Rx, Ry, Child_Subst, Sol0, omega = 1, 
                 a = R, R = 5, [Toy[], t = 0..40]);
plot(Toy_Plot, scaling = constrained, color = black);
Toy_Plot := subs(Rx, Ry, Child_Subst, Sol1, omega = 1, alpha = arctan(-2)+Pi,
                  a = sqrt(125), R = 5, [Toy[], t  = 0..200]):
plot(Toy_Plot, scaling = constrained, color = black);
Toy_Plot := subs(Rx, Ry, Child_Subst, Sol2, omega = 1, alpha = 0,
                 a = 4.95, R = 5, [Toy[], t = 0..40]):
plot(Toy_Plot, scaling = constrained, color = black);
a := 'a': Vx := 'Vx': Vy := 'Vy':
Child_Subst := X(t) = Xo + Vx*t, Y(t) = Yo + Vy*t:
Das := eval(subs(Child_Subst, Daxial));
Sol := dsolve(Das, phi(t));
X := t -> t:  Y := t -> 5*sin(t):  a:= 10:
Daxial;
F := dsolve({Daxial, phi(0) = Pi/2}, phi(t),numeric);
N := 200: L := a:
T:=evalf(6*Pi/N*[$0..N]):
SF := [seq(rhs(F(t)[2]), t = T)]:
plot(zip((u,v)->[X(u) + L*cos(v), Y(u) + L*sin(v)],T,SF));
with(plots): 
TS := display(seq(plot(zip((u,v)->[X(u) + L*cos(v), 
      Y(u) + L*sin(v)], T[1..j],SF[1..j]),color=blue,thickness=2), j = 1..N+1),insequence=true):
BS := display(zip((u,v)->plot([[X(u), Y(u)], 
      [X(u) + L*cos(v), Y(u) + L*sin(v)]],color=red),T,SF),insequence=true):
PS := display(seq(plot([seq([X(u), Y(u)],u=T[1..j])],color=black),j=1..N+1),insequence=true):
display({TS,BS,PS},view=[-2..20,-5..15]);