#==============================================================================
# Curves package
#==============================================================================
# Marc A. Murison
# Astronomical Applications Dept.
# U.S. Naval Observatory
# 3450 Massachusetts Ave., NW
# Washignton, DC 20392
# murison@aa.usno.navy.mil
# http://aa.usno.navy.mil/murison/
#==============================================================================

macro(circle=curves['circle']):
macro(ellipse=curves['ellipse']):
macro(parabola=curves['parabola']):
macro(hyperbola=curves['hyperbola']):
macro(lemniscate=curves['lemniscate']):
macro(cissoid=curves['cissoid']):
macro(tractrix=curves['tractrix']):
macro(clotharg=curves['clotharg']):
macro(clothoid=curves['clothoid']):
macro(folium=curves['folium']):
macro(cassinian=curves['cassinian']):
macro(logspiral=curves['logspiral']):
macro(catenary=curves['catenary']):
macro(astroid=curves['astroid']):
macro(cardioid=curves['cardioid']):
macro(cycloid=curves['cycloid']):
macro(deltoid=curves['deltoid']):
macro(lissajous=curves['lissajous']):
macro(torusknot=curves['mtorusknot']):
macro(twicubic=curves['twicubic']):
macro(helix=curves['helix']):
macro(ParametrizedToImplicit=curves['ParametrizedToImplicit']):
macro(tangent=curves['tangent']):
macro(unitnormal=curves['unitnormal']):
macro(binormal=curves['binormal']):
macro(curvature=curves['curvature']):
macro(torsion=curves['torsion']):
macro(arclength=curves['arclength']):
macro(evolute=curves['evolute']):
macro(involute=curves['involute']):
macro(radial=curves['radial']):
macro(pedal=curves['pedal']):
macro(inversion=curves['inversion']):

#------------------------------------------------------------------------------
# Plane curves (parametric forms, except for Cassinian Oval)
#------------------------------------------------------------------------------
`curves/circle`     := (a,t)     -> vec(a*cos(t), a*sin(t));
`curves/ellipse`    := (a,b,t)   -> vec(a*cos(t), b*sin(t));
`curves/parabola`   := (f,t)     -> vec( 2*f*t, f*t^2 );
`curves/hyperbola`  := (a,b,t)   -> vec( a*cosh(t), b*sinh(t) );
`curves/lemniscate` := (a,t)     -> vec( a*cos(t)/(1+sin(t)^2), a*sin(t)*cos(t)/(1+sin(t)^2) );
`curves/cissoid`    := (a,t)     -> vec( 2*a*t^2/(1+t^2), 2*a*t^3/(1+t^2) );
`curves/tractrix`   := (a,t)     -> vec( a*sin(t), a*(cos(t)+log(tan(t/2))) );
`curves/clotharg`   := (n,a,t)   -> vec( a*sin(t^(n+1)/(n+1)), a*cos(t^(n+1)/(n+1)) );
`curves/clothoid`   := (n,a,t)   -> map( int, clotharg(n,a,u), u=0..t );
`curves/folium`     := (t)       -> vec( 3*t/(1+t^3), 3*t^2/(1+t^3) );
#curves[cassinian]  := (a,b,x,y) -> (x^2+y^2+a^2)^2 - (2*a*x)^2 - b^4;
`curves/cassinian`  := (a,b,t)   -> evalm( (a^2*cos(2*t) + sqrt(b^4 - a^4*sin(2*t)^2))*vec(cos(t),sin(t)) );
`curves/logspiral`  := (a,b,t)   -> vec( a*exp(b*t)*cos(t), a*exp(b*t)*sin(t) );
`curves/catenary`   := (a,t)     -> vec( t, a*cosh(t/a) );
`curves/astroid`    := (n,a,b,t) -> vec( (a*cos(t))^n, (b*sin(t))^n );
`curves/cardioid`   := (a,t)     -> vec( 2*a*cos(t)*(1+cos(t)), 2*a*sin(t)*(1+cos(t)) );
`curves/cycloid`    := (a,b,t)   -> vec( a*t-b*sin(t), a-b*cos(t) );
`curves/deltoid`    := (a,t)     -> vec( 2*a*cos(t)*(1+cos(t))-a, 2*a*sin(t)*(1-cos(t)) );
`curves/lissajous`  := (n,m,a,b,phi,t) -> vec( a*sin(n*t), b*sin(m*t+phi) );

#------------------------------------------------------------------------------
# Space curves (parametric forms)
#------------------------------------------------------------------------------
`curves/torusknot` := (a,b,c,p,q,t) -> vec( (a+b*cos(q*t))*cos(p*t),
                                            (a+b*cos(q*t))*sin(p*t),
                                            c*sin(q*t) );
`curves/twicubic`  := (t) -> vec( t, t^2, t^3 );
`curves/helix`     := (curve2D,a,t) -> vec( curve2D[1], curve2D[2], a*t );


#------------------------------------------------------------------------------
# Given a parameterized plane curve, convert it to an implicit curve in (x,y)
#------------------------------------------------------------------------------
`curves/ParametrizedToImplicit` := proc(curve::{vector(2),vector(3)}, t::name, vars::list(name))
    local S, k, c;
    c := convert(curve,list);
    S := { };
    for k from 1 to nops(c) do
        S := S union { vars[k]=c[k] };
    od;
    eliminate( S, t ); 
    op(%[2]) = 0;
end proc;

#------------------------------------------------------------------------------
# Unit tangent vector to a curve in 3D Euclidean space
#
# The user will have to evalm() the result in order to perform
# operations (like simplify) on it.
#------------------------------------------------------------------------------
`curves/tangent` := proc( curve::{vector(2),vector(3)}, t::name )
    local d;
    d := map( diff, curve, t );
    evalm(d/mag(d));
end proc;

#------------------------------------------------------------------------------
# Unit normal vector to a curve in 3D Euclidean space
#
# The user will have to evalm() the result in order to perform
# operations (like simplify) on it.
#------------------------------------------------------------------------------
`curves/unitnormal` := proc( curve::{vector(2),vector(3)}, t::name )
    local d1, d2, cp, un;
    d1 := map( diff, curve, t );
    d2 := map( diff, curve, t$2 );
    cp := cross(d2,d1);
    cross( d1, eval(cp))/mag(eval(cp))/mag(d1);
    un := evalm( cross(d1,cp)/(mag(cp)*mag(d1)) );
    if size(curve)=2 then
        vec( un[1], un[2] );
    else
        evalm(un);
    end if;
end proc;

#------------------------------------------------------------------------------
# Unit binormal vector to a curve in 3D Euclidean space
#
# The user will have to evalm() the result in order to perform
# operations (like simplify) on it.
#------------------------------------------------------------------------------
`curves/binormal` := proc( curve::{vector(2),vector(3)}, t::name )
    local d1, d2, cp;
    d1 := map( diff, curve, t );
    d2 := map( diff, curve, t$2 );
    cp := cross(d1,d2);
    evalm(cp/mag(cp));
end proc;

#------------------------------------------------------------------------------
# Curvature of 2D or 3D curve
#------------------------------------------------------------------------------
`curves/curvature` := proc( curve::{vector(2),vector(3)}, t::name )
    local d1, d2;
    d1 := map( diff, curve, t );
    d2 := map( diff, d1, t );
    mag(cross(d1,d2))/mag(d1)^3;
end proc;

#------------------------------------------------------------------------------
# Torsion of 3D curve
#------------------------------------------------------------------------------
`curves/torsion` := proc( curve::{vector(2),vector(3)}, t::name )
    local d1, d2, d3, cp;
    d1 := map( diff, curve, t );
    d2 := map( diff, curve, t$2 );
    d3 := map( diff, curve, t$3 );
    cp := cross(d1,d2);
    dot(eval(cp),d3);
    dot(eval(cp),d3)/mag(eval(cp))^2;
end proc;

#------------------------------------------------------------------------------
# arc length of a 2D or 3D curve (parametric form)
#
# usage:
#   arclength( curve, t );
#   arclength( curve, t=a..b );
#   arclength( curve, t=f );
#   arclength( curve, t=f..g );
# where a and b are Maple names, and f and g are numeric
#
# Warning: simplify( , assume=positive) is applied to symbolic results
#------------------------------------------------------------------------------
`curves/arclength` := proc( curve::{vector(2),vector(3)}, t::{name,name=range,name=numeric})
    local d, s, z;

    #grab independent var and differentiate the curve
    if type(t,name) then s := t; else s := lhs(t); end if;
    d := map( diff, curve, s );
     
    if type(t,name=numeric) then
        evalf( int( mag(d), s=0..rhs(t) ) );
    elif type(t,name=range(numeric)) then
        evalf( int( mag(d), t ) );
    elif type(t,name) then
        z := int( mag(d), t );      #try to integrate symbolically
        if evalb(op(0,z)=int) then  #failed -- return unevaluated
            Int( mag(d), t );
        else                        #succeeded
            simplify( z, assume=positive );
        end if;
    elif type(t,name=range(name)) then
        z := int( mag(d), s );      #try to integrate symbolically
        if evalb(op(0,z)=int) then  #failed -- return unevaluated
            Int( mag(d), t );
        else                        #succeeded
            z := simplify( z, assume=positive );
            subs( s=op(2,rhs(t)), z ) - subs( s=op(1,rhs(t)), z );
        end if;
    end if;
 
end proc;

#------------------------------------------------------------------------------
# evolute
#
# Warning: simplify( ..., assume=positive) is applied to symbolic results
#------------------------------------------------------------------------------
`curves/evolute` := proc ( curve::{vector(2),vector(3)}, t::name )
    local r, N;
    r := 1/curvature(curve,t);
    N := unitnormal(curve,t);
    if size(curve)=2 then
        N := vec( N[1], N[2] );
    end if;
    map( simplify, evalm(curve + r*N), assume=positive );
end proc;


#------------------------------------------------------------------------------
# involute
#
# Warning: simplify( ..., assume=positive) is applied to symbolic results
#------------------------------------------------------------------------------
`curves/involute` := proc( curve::{vector(2),vector(3)}, t::name )
    curve - arclength(curve,t)*tangent(curve,t);
    map( simplify, evalm(%), assume=positive );
end proc;

#------------------------------------------------------------------------------
# radial
#
# Warning: simplify( ..., assume=positive) is applied to symbolic results
#------------------------------------------------------------------------------
`curves/radial` := proc ( curve::{vector(2),vector(3)}, t::name )
    local r, N;
    r := 1/curvature(curve,t);
    N := unitnormal(curve,t);
    if size(curve)=2 then
        N := vec( N[1], N[2] );
    end if;
    map( simplify, evalm(r*N), assume=positive );
end proc;

#------------------------------------------------------------------------------
# pedal
#
# Warning: simplify( ..., assume=positive) is applied to symbolic results
#------------------------------------------------------------------------------
`curves/pedal` := proc ( curve::{vector(2),vector(3)}, t::name, refpoint::{vector(2),vector(3)} )
    local T, N;
    T := tangent(curve,t);
    N := unitnormal(curve,t);
    if size(curve)=2 then
        N := vec( N[1], N[2] );
    end if;
    map( simplify, evalm(refpoint + dot(curve-refpoint,N)*N), assume=positive );
end proc;

#------------------------------------------------------------------------------
# inverse
#
# r = inversion radius
#
# Warning: simplify( ..., assume=positive) is applied to symbolic results
#------------------------------------------------------------------------------
`curves/inversion` := proc ( curve::{vector(2),vector(3)}, t::name, 
                             refpoint::{vector(2),vector(3)}, r::algebraic )
    refpoint + r^2*(curve - refpoint)/dot(evalm(curve - refpoint),evalm(curve - refpoint));
    map( simplify, evalm(%), assume=positive );
end proc;


#==============================================================================
savelibname := CURVES_PATH;

savelib(`curves/circle`,"curves/circle.m"):
savelib(`curves/ellipse`,"curves/ellipse.m"):
savelib(`curves/parabola`,"curves/parabola.m"):
savelib(`curves/hyperbola`,"curves/hyperbola.m"):
savelib(`curves/lemniscate`,"curves/lemniscate.m"):
savelib(`curves/cissoid`,"curves/cissoid.m"):
savelib(`curves/tractrix`,"curves/tractrix.m"):
savelib(`curves/clotharg`,"curves/clotharg.m"):
savelib(`curves/clothoid`,"curves/clothoid.m"):
savelib(`curves/folium`,"curves/folium.m"):
savelib(`curves/cassinian`,"curves/cassinian.m"):
savelib(`curves/logspiral`,"curves/logspiral.m"):
savelib(`curves/catenary`,"curves/catenary.m"):
savelib(`curves/astroid`,"curves/astroid.m"):
savelib(`curves/cardioid`,"curves/cardioid.m"):
savelib(`curves/cycloid`,"curves/cycloid.m"):
savelib(`curves/deltoid`,"curves/deltoid.m"):
savelib(`curves/lissajous`,"curves/lissajous.m"):
savelib(`curves/torusknot`,"curves/mtorusknot.m"):
savelib(`curves/twicubic`,"curves/twicubic.m"):
savelib(`curves/helix`,"curves/helix.m"):
savelib(`curves/ParametrizedToImplicit`,"curves/ParametrizedToImplicit.m"):
savelib(`curves/tangent`,"curves/tangent.m"):
savelib(`curves/unitnormal`,"curves/unitnormal.m"):
savelib(`curves/binormal`,"curves/binormal.m"):
savelib(`curves/curvature`,"curves/curvature.m"):
savelib(`curves/torsion`,"curves/torsion.m"):
savelib(`curves/arclength`,"curves/arclength.m"):
savelib(`curves/evolute`,"curves/evolute.m"):
savelib(`curves/involute`,"curves/involute.m"):
savelib(`curves/radial`,"curves/radial.m"):
savelib(`curves/pedal`,"curves/pedal.m"):
savelib(`curves/inversion`,"curves/inversion.m"):

macro(circle=circle):
macro(ellipse=ellipse):
macro(parabola=parabola):
macro(hyperbola=hyperbola):
macro(lemniscate=lemniscate):
macro(cissoid=cissoid):
macro(tractrix=tractrix):
macro(clotharg=clotharg):
macro(clothoid=clothoid):
macro(folium=folium):
macro(cassinian=cassinian):
macro(logspiral=logspiral):
macro(catenary=catenary):
macro(astroid=astroid):
macro(cardioid=cardioid):
macro(cycloid=cycloid):
macro(deltoid=deltoid):
macro(lissajous=lissajous):
macro(torusknot=mtorusknot):
macro(twicubic=twicubic):
macro(helix=helix):
macro(ParametrizedToImplicit=ParametrizedToImplicit):
macro(tangent=tangent):
macro(unitnormal=unitnormal):
macro(binormal=binormal):
macro(curvature=curvature):
macro(torsion=torsion):
macro(arclength=arclength):
macro(evolute=evolute):
macro(involute=involute):
macro(radial=radial):
macro(pedal=pedal):
macro(inversion=inversion):

curves[circle]     := 'readlib('`curves/circle`')':
curves[ellipse]    := 'readlib('`curves/ellipse`')':
curves[parabola]   := 'readlib('`curves/parabola`')':
curves[hyperbola]  := 'readlib('`curves/hyperbola`')':
curves[lemniscate] := 'readlib('`curves/lemniscate`')':
curves[cissoid]    := 'readlib('`curves/cissoid`')':
curves[tractrix]   := 'readlib('`curves/tractrix`')':
curves[clotharg]   := 'readlib('`curves/clotharg`')':
curves[clothoid]   := 'readlib('`curves/clothoid`')':
curves[folium]     := 'readlib('`curves/folium`')':
curves[cassinian]  := 'readlib('`curves/cassinian`')':
curves[logspiral]  := 'readlib('`curves/logspiral`')':
curves[catenary]   := 'readlib('`curves/catenary`')':
curves[astroid]    := 'readlib('`curves/astroid`')':
curves[cardioid]   := 'readlib('`curves/cardioid`')':
curves[cycloid]    := 'readlib('`curves/cycloid`')':
curves[deltoid]    := 'readlib('`curves/deltoid`')':
curves[lissajous]  := 'readlib('`curves/lissajous`')':
curves[torusknot]  := 'readlib('`curves/mtorusknot`')':
curves[twicubic]   := 'readlib('`curves/twicubic`')':
curves[helix]      := 'readlib('`curves/helix`')':
curves[ParametrizedToImplicit] := 'readlib('`curves/ParametrizedToImplicit`')':
curves[tangent]    := 'readlib('`curves/tangent`')':
curves[unitnormal] := 'readlib('`curves/unitnormal`')':
curves[binormal]   := 'readlib('`curves/binormal`')':
curves[curvature]  := 'readlib('`curves/curvature`')':
curves[torsion]    := 'readlib('`curves/torsion`')':
curves[arclength]  := 'readlib('`curves/arclength`')':
curves[evolute]    := 'readlib('`curves/evolute`')':
curves[involute]   := 'readlib('`curves/involute`')':
curves[radial]     := 'readlib('`curves/radial`')':
curves[pedal]      := 'readlib('`curves/pedal`')':
curves[inversion]  := 'readlib('`curves/inversion`')':
savelib(curves,"curves.m");

savelibname := 'savelibname';
#==============================================================================