with(student): 

 tanl := (x,y) -> (sqrt((-diff(x,t))^2+(-diff(y,t))^2)):
 norx := (x,y) -> ((-diff(x,t))/tanl(x,y)):

                                           
 nory := (x,y) -> ((-diff(y,t))/tanl(x,y)):

                  
 arc := (x,y,tstart) -> (simplify(evalf(Int(subs(t=n,(tanl(x,y))), 
n=tstart..t)))):

                                           
                                          
 involute := (x,y,tstart,tend) -> (simplify(x + (evalf(Int(subs(t=n,abs(sqrt((-diff(x,t))^2+(-diff(y,t))^2))),n=tstart..t)))*((-diff(x,t))/(sqrt((-diff(x,t))^2+(-diff(y,t))^2)))), simplify(y + (evalf(Int(subs(t=n,abs(sqrt((-diff(x,t))^2+(-diff(y,t))^2)))
,n=tstart..t)))*((-diff(y,t))/(sqrt((-diff(x,t))^2+(-diff(y,t))^2))))):
                        
 paraminvolute := (r,tstart,tend) -> involute(r*cos(t),r*sin(t),tstart,tend):


 inv := (x,y,tstart,tend) -> (simplify(x + arc(x,y,tstart)*norx(x,y)), 
simplify(y + arc(x,y,tstart)*nory(x,y))):

  
 paraminv := (r,tstart,tend) -> inv(r*cos(t),r*sin(t),tstart,tend):
  evolute2:= (x,y) ->[ simplify( x- ( diff(x,t)^2 + diff(y,t)^2 
)*diff(y,t)/(diff(y,t,t)*diff(x,t) - diff(x,t,t)*diff(y,t))), simplify(y 
+ ( diff(x,t)^2 + diff(y,t)^2 )*diff(x,t)/(diff(y,t,t)*diff(x,t) - diff(x,t,t)*diff(y,t)))]:
evolute3 := (x,y) -> ({evolute2(x,y),[x,y]}):
evolute4 := (x,y,t1,t2) -> ({[evolute(x,y),t=t1..t2],[x,y,t=t1..t2]}):
 evolute:= (x,y) ->( simplify( x- ( diff(x,t)^2 + diff(y,t)^2 
)*diff(y,t)/(diff(y,t,t)*diff(x,t) - diff(x,t,t)*diff(y,t))), simplify(y + ( diff(x,t)^2 + diff(y,t)^2 )*diff(x,t)/(diff(y,t,t)*diff(x,t) - diff(x,t,t)*diff(y,t)))):


Double_Evolute(x,y,t1,t2) -> ({[evolute(evolute(x,y)),t=t1..t2],[evolute(x,y),t=t1..t2],[x,y,t=t1..t2]}):

devolute(x,y) -> ({evolute2(evolute(x,y)),evolute2(x,y),[x,y]}): 
              
 paramevo := (x) -> (evolute(x*cos(t),x*sin(t))):
 involutef := (x,y,tstart,tend) -> (simplify(x + 
(evalf(simpson(subs(t=n,abs(sqrt((-diff(x,t))^2+(-diff(y,t))^2))),n=tstart..t)))*((-diff(x,t))/(sqrt((-diff(x,t))^2+(-diff(y,t))^2)))), simplify(y + (evalf(simpson(subs(t=n,abs(sqrt((-diff(x,t))^2+(-diff(y
,t))^2))),n=tstart..t)))*((-diff(y,t))/(sqrt((-diff(x,t))^2+(-diff(y,t))^2))))):

d0:= `evolute(x,y) - Returns the evolute of the parameterized function  `:
d1:= `               f(t) -> (x,y).                                     `:
d2:= `evolute2(x,y) - Same as evolute(x,y), only the result has         `:
d3:= `                brackets around it.                               `:
d4:= `evolute3(x,y) - Returns a the evolute and the original function   `:
d5:= `                in a form easily incorporated into the            `:
d6:= `                ParamPlot function:                               `: 
d7:= `                (ParamPlot(evolute3(x,y),t=0..2*Pi,color=green)). `:
e := `evolute4(x,y,tstart,tend) - Returns the evolute and the original  `:
e1:= `                function in a format easily incorporated into a   `:
e2:= `                plot function: plot(evolute(x,y,t0,t1),color=red).`:
#e3:= `Double_Evolute(x,y,tstart,tend) - Same as evolute4, except it     `:
#e4:= `                plots the function, the first evolute and the 2nd `:
#e5:= `                evolute.`:
d8:= `paramevo(t) - Returns the evolute of the polar function           `:
d9:= `              parameterized in terms of t (not theta).            `:
d10:=`involute(x,y,tstart) - Returns the involute of the parameterized  `:
d11:=`                       function f(t) -> (x,y) starting at t=tstart`:
d12:=`                       with parenthesis - uses the absolute value `:
d13:=`                       inside the integral.                       `:
d14:=`paraminvolute(t,tstart) - Uses the involute function on the polar `:
d15:=`                          function (in terms of t).               `:
d16:=`inv(x,y,tstart) - Same as the involute without the absolute value.`:
d17:=`paraminv(t,tstart) - Same as paraminvolute except it uses the inv `:
d18:=`                     function instead of the involute function.   `:

 
d0; d1; d2; d3; d4; d5; d6; d7;e; e1;e2; d8; d9; d10; d11; d12; 
d13; d14; d15; d16; d17; d18;

#########################
#
#  The function ParamPlot was taken from the following package:
#

##########################
# /math/mathlab/CalcP/Vike/MapleVR2/macro
###
# Convert from 5.2 to Release 3
##
##      CalcP, a package for calculus
##
##      Date    : October 21, 1993
##
##      Authors : Willaim Farr
##                Michael VanGeel
##                Mathematical Sciences
##                Worcester Polytechnic Institute
##                Worcester, MA 01564 
##
##                bfarr@wpi.wpi.edu
##

##########################
# /math/mathlab/CalcP/Vike/MapleVR2/ParamPlot.tst
###
# Convert from 5.2 to Release 3
ParamPlot := proc(flist:{vector(algebraic),list(algebraic),
              set(list(algebraic)),set(vector(algebraic))},
                  trange:name=range(constant))
 local xloc, td, tstart, plotlist, numpts, opargs, n, funclist,
       funcset,i,plotset;
 options `Copyright 1993 by William Farr and Michael VanGeel`;
   if type(flist,set) then 
     funcset := map(convert,flist,list)
   else
     funcset := {convert(flist,list)}
   fi;
   numpts := 100;
   opargs := [];
   td := op(1,trange);
   for i from 1 to nops(funcset) do 
     funclist := op(i,funcset);
     if nops(funclist)<>2 then
      ERROR(`first argument must be a two component vector or list.`) fi;
     if nops(indets(funclist,name) minus indets(funclist,constant))>1 then 
        ERROR(`the first argument has parameters that must be defined.`) fi;
     if not(member(td,indets(funclist,name))) then 
        ERROR(`second argument variable not present in the first.`) fi;
   od;
   if nargs > 2 then
      for n from 3 to nargs do
         if op(1,args[n])=numpoints then
            numpts := op(2,args[n])
          elif op(1,args[n])=frames then
            ERROR(`the frames option cannot be changed.`)
          else 
            opargs := [opargs[],args[n]];
         fi
       od
   fi;
   tstart := op(1,op(2,trange));
   plotset := {};
   for i from 1 to nops(funcset) do
     funclist := op(i,funcset);
     plotlist := subs(td = tstart + xloc*(td - tstart),funclist);
     plotlist := [op(plotlist),xloc=0..1];
     plotset := plotset union {plotlist};
   od;
   plots[animate](plotset,trange,frames=32,numpoints=numpts,opargs[]);
end: