%% Problem Set C: Problem 14
%% (a) parametric equations

syms a b t real
assume(a>0 & b>0)
r = [(a+b)*cos(t) - b*cos((a+b)*t/b), (a+b)*sin(t) - b*sin((a+b)*t/b)];

%
%% (b) graphs
% Let's fix a = 1 (no loss of generality, since we can rescale) and try b = 0.25, 0.5, and 1.

r1 = subs(r, [a, b], [1, 1/4]);
r2 = subs(r, [a, b], [1, 1/2]);
r3 = subs(r, [a, b], [1, 1]);
fplot(r1(1), r1(2), [0, 4*pi], 'r'), hold on
fplot(r2(1), r2(2), [0, 4*pi], 'b')
fplot(r3(1), r3(2), [0, 4*pi], 'g'), hold off
axis equal

%
%% (c) speed and curvature

veclength = @(x) sqrt(x*transpose(x));
vel = diff(r,t);
speed = veclength(vel);
accel = diff(vel, t);
vta = simplify(cross([vel,0], [accel,0]));
curvature = simplify(veclength(vta)/speed^3);
sqspeed = simplify(speed^2, 'Steps', 50);
pretty(sqspeed)
sqcurv = simplify(curvature^2, 'Steps', 50);
pretty(sqcurv)

% We can see from this that the curvature vanishes only when ``#1" vanishes. The curvature blows up (goes to infinity) and the speed vanishes when ``#2" and ``#3" vanish simultaneously, i.e., when sin(t*(a+b)/b) = sin(t) and cos(t*(a+b)/b) = cos(t). This always happens at t=0 (the cusp on the right on the figure) but whether it happens at other points or not depends on whether a/b is rational or not and what it's exact value is.
%
%% (d) the evolute

UT = [vel,0]/speed;
UB = simplify(vta/veclength(vta));
UN = cross(UB, UT); UN = UN(1:2);
evolute = simplify(r + UN/curvature);
pretty(evolute)

%
% Now let's specialize the case a = 3, b = 1. We will plot the fixed circle in green, the epicycloid curve in red, and the evolute in blue.

r = simplify(subs(r , [a, b], [3, 1]));
evolute = simplify(subs(evolute , [a, b], [3, 1]));
fplot(r(1), r(2), [0, 4*pi], 'r'), hold on
fplot(3*cos(t), 3*sin(t), [0, 2*pi], 'g')
fplot(evolute(1), evolute(2), [0, 4*pi], 'b'), hold off
axis equal

% A very pretty picture! The evolute looks like a shrunken and rotated version of the same curve.