(* Copyright 1993 Alfred Gray *)
BeginPackage["CSPROGS`"]

(* Arc Length Functions *)

    arclengthprime::usage="t->arclengthprime[alpha][t] is the derivative of
the arclength function of the curve alpha:(a,b)->R^n."

    arclength::usage="t->arclength[alpha][t] is the arclength function of
the curve alpha:(a,b)->R^n computed using Integrate."

    length::usage="length[a,b][alpha] is the length of the curve
alpha:(a,b)->R^n computed symbolically."

    ndslength::usage="ndlsength[a,b][alpha][t] is the length of the curve
alpha:(a,b)->R^n over the interval [a,t] computed using NDSolve over
the interval [a,b]."

    nlength::usage="nlength[a,b][alpha] is the length of the curve
alpha:(a,b)->R^n computed using NIntegrate."

(* Brioschi Curvature *)

    brioschicurvature::usage="{u,v}->brioschicurvature[x][u,v] is the
Gaussian curvature of the patch x:U->R^n computed using
Brioschi's formula.  {u,v}->brioschicurvature[ee,ff,gg][u,v] is the
Gaussian curvature of the metric \n
ds^2 = ee*du^2 + 2*ff*du dv + gg*dv^2."

(* Christoffel Symbols *)

    christoffel::usage="{u,v}->christoffel[i,j,k][x][u,v] is a Christoffel
symbol of the metric tensor of the patch x:U->R^n.
{u,v}->christoffel[i,j,k][ee,ff,gg][u,v] is a Christoffel symbol
of the metric \n
ds^2 = ee*du^2 + 2*ff*du dv + gg*dv^2."

(* Complex to Real to Complex\label{ComplextoRealtoComplex *)

    ctor::usage="ctor[z] converts a complex number z to an ordered pair
of real numbers whose components are the real and imaginary parts of z;
also works for lists of complex numbers."

    rtoc::usage="rtoc[{x,y}] converts an ordered pair {x,y} of
real numbers to a complex number whose real part is x and
whose imaginary part is y;
also works for lists of ordered pairs."

    cdot::usage="cdot[z,w] is the dot product of complex vectors
z and w."

(* Complex Structure *)

    J::usage="J is the linear map J:R^2->R^2 defined by
that maps {p1,p2} onto {-p2,p1}."

(* Cross Product *)

    Cross::usage="Cross[u,v] is the cross product of the vectors
v={v1,v2,v3} and w={w1,w2,w3}."

(* Construction of Plane Curves from
Plane Curves  *)

    circle3point::usage="circle3point[{p1,q1},{p2,q2},{p3,q3}] is
the circle in the plane passing through the points
{p1,q1},{p2,q2},{p3,q3}."

    evolute::usage="t->evolute[alpha][t] is the
evolute of the plane curve alpha.  evolute[n,alpha][t] is the
nth evolute of the plane curve alpha."

     involute::usage="t->involute[alpha,c][t] is the involute of alpha
starting at alpha[c], that is, the curve traced by a point on a thread
attached to alpha at alpha[c] that is kept taut as it is unwound from
alpha. t->involute[alpha][t] is t->involute[alpha,0][t]."

     ndsinvolute::usage="t->ndsinvolute[a,b][alpha][t] is the involute of alpha starting at alpha[a] computed over the interval [a,b] using NDSolve.
It is the curve traced by a point on a thread attached to alpha at alpha[a] that is kept taut as it is unwound from alpha."

    normalline::usage="t->normalline[alpha][s,t] is the normal line of
length s emanating from a plane curve alpha at alpha[t]."

    osculatingcircle::usage="t->osculatingcircle[alpha][t] is the osculating
circle to the plane curve alpha at alpha[t]."

    tangentline::usage="t->tangentline[alpha][s,t] is the tangent line of
length s emanating from the plane curve alpha at alpha[t]."

    plottl::usage="plottl[curve,t0,{tmin,tmax},opts]} plots
the tangent line emanating from curve."

    parcurve::usage="t->parcurve[alpha][s][t] is one component of the
curve parallel to a plane curve alpha at a distance Abs[s],
where Abs[s] is small. The other component is
t->parcurve[alpha][-s][t]."

    pedal::usage="t->pedal[alpha][t] is the pedal curve of the
plane curve alpha with respect to the origin.
pedal[{a,b},alpha][t] is the pedal curve of the plane curve alpha
with respect to the point {a,b}."

(* Construction of Space Curves from
Plane Curves  *)

    helical::usage="t->helical[alpha,a][t] is the helical curve
in R^3 with slant a over the plane curve alpha.
helical[alpha][f][t] is the helical curve in R^3 with
slant function f over the plane curve alpha."

    inscurv::usage="t->inscurv[n,a,alpha][t] is the curve that results when
the value a is inserted as the nth coordinate into the parametrization of
the curve alpha."

    kappa23d::usage="t->kappa23d[alpha][t] is the space curve whose
first two components are those of the plane curve alpha, and whose third component is the signed curvature of alpha.
t-> kappa23d[mag,alpha][t]  is the space curve whose
first two components are those of the plane curve alpha, and whose third component is the signed curvature of alpha multiplied by mag."

    writhe::usage="t->writhe[alpha][t] is a curve in R^3 that is a
twisted generalization of the plane curve t->alpha[t]."

(* Construction of Surfaces from
Plane Curves  *)

    cone::usage="{u,v}->cone[{p1,p2,p3},gamma][u,v] is the cone
with vertex {p1,p2,p3} over the curve gamma:(a,b)->R^3."

    cylinder::usage="cylinder[{q1,q2,q3},gamma][u,v] is the cylinder
in the direction {q1,q2,q3} over the curve gamma:(a,b)->R^3."

    genhel::usage="{u,v}->genhel[alpha][sl][u,v] is the generalized helicoid
with slant  sl  generated by the plane curve alpha."

    gensurfrev::usage="{u,v}->gensurfrev[alpha,gamma][u,v] is the
generalized surface of revolution generated by the plane curves alpha
and gamma."

    ruled::usage="{u,v}->ruled[alpha,gamma][u,v] is the ruled patch 
generated by the curves alpha and gamma."

    surfast::usage="{u,v}->surfast[alpha][u,v] is the astroidial patch 
generated by the plane curve alpha. surfast[n,alpha][u,v] is the
astroidial patch  of degree n generated by the plane curve alpha."

    surfrev::usage="{u,v}->surfrev[alpha][u,v] is the
surface of revolution generated by the plane curve alpha."

    twist::usage="{u,v}->twist[a,b,gamma][u,v] is the twisted patch 
generated by a curve gamma:(c,d)->R^2. The Moebius strip and Klein bottle
are special cases."

(* Construction of Surfaces from
Space Curves  *)

    binormaldev::usage="{u,v}->binormaldev[alpha][u,v] is the
developable surface generated by the curve alpha:(a,b)->R^3 and
its binormal vector field."

    normaldev::usage="{u,v}->normaldev[alpha][u,v] is the
developable surface generated from the curve alpha:(a,b)->R^3 and
its normal vector field."

   perpdev::usage="{u,v}->perpdev[phi,alpha][u,v] is a 1-parameter
family of developable surfaces generated by the curve alpha:(a,b)->R^3
and some normal vector field to alpha.
Special cases are {u,v}->perpdev[0,alpha][u,v], which is
the normal developable of alpha, and
{u,v}->perpdev[Pi/2,alpha][u,v], which is
the binormal developable of alpha."

    seashell::usage="{t,theta}->seashell[gamma][r][t,theta] is the
tube of variable radius r t about the curve gamma:(a,b)->R^3."

    tandev::usage="{u,v}->tandev[alpha][u,v] is the
developable surface generated by the curve alpha:(a,b)->R^n and its
tangent vector field."

    tubecurve::usage="{t,theta}->tubecurve[gamma][r][t,theta] is the
tube of radius r about the curve gamma:(a,b)->R^3."

    vartube::usage="{t,theta}->vartube[gamma][r][t,theta] is the
tube of variable radius r[t] about the curve gamma:(a,b)->R^3."

(* Construction of Surfaces from
Surfaces  *)

    parsurf::usage="{u,v}->parsurf[x][r][u,v] is the
parallel patch  at a distance r from a patch  x:U->R^3."

    pedsurf::usage="{u,v}->pedsurf[{a,b,c},x][u,v] is the pedal patch
of the patch x with respect to the point {a,b,c}.
{u,v}->pedsurf[x][u,v] is the pedal patch  of the
patch x with respect to the origin."

    stretch::usage="{u,v}->stretch[a,b,c,surf][u,v] is the
patch surf, but stretched by a,b,c."

    tangentplane::usage="{u,v}->tangentplane[x][u0,v0][u,v] is the
tangent plane to the patch  x at x[u0,v0]."

    plottp::usage="plottp[x,{u0,v0},{umin,umax},{vmin,vmax},opts]
plots a patch  x and its tangent plane at x[u0,v0]."

(* Curvature of Surfaces (Fast) *)

    gcurvature::usage="{u,v}->gcurvature[x][u,v] is the Gaussian curvature
of the patch x:U->R^3."

    mcurvature::usage="{u,v}->mcurvature[x][u,v] is the mean curvature of the
patch x:U->R^3."

(* Curvature of Surfaces (Slow) *)

    gaussiancurvature::usage="{u,v}->gaussiancurvature[x][u,v] is the
Gaussian curvature of the patch x:U->R^3."

    meancurvature::usage="{u,v}->meancurvature[x][u,v] is the
mean curvature of the patch x:U->R^3."

(* Curvature and Torsion of Curves *)

    kappa::usage="t->kappa[alpha][t] is the
unsigned curvature of the curve alpha:(a,b)->R^3."

    kappan::usage="t->kappan[alpha][t] is the
unsigned curvature of the curve alpha:(a,b)->R^n."

    kappa2::usage="t->kappa2[alpha][t] is the signed curvature of
the plane curve alpha."

    kappa2c::usage="t->kappa2c[alpha][t] is the signed curvature of the
complex plane curve alpha."

    kappa2implicit::usage="{x,y}->kappa2implicit[g][x,y] is the signed
curvature of the plane curve g[x,y]==0."

    kappa2polar::usage="t->kappa2polar[rho][t] is the signed
curvature of the plane curve whose radius function is rho."

    tau::usage="t->tau[alpha][t] is the
torsion of the curve alpha:(a,b)->R^3."

    totalacsalrcurvature::usage="totalacsalrcurvature[a,b,alpha] is the
total curvature of the curve alpha:(a,b)->R^2."

    turningangle::usage="t->turningangle[alpha][t] is the
turning angle of the curve alpha:(a,b)->R^2."

    turningnumber::usage="turningnumber[alpha,{a,b},opts] is the
turning number of the closed curve alpha:(a,b)->R^2."

(* Euclidean Motions Acting on
Plane Curves *)

    tl::usage="t->tl[p,alpha][t] is the curve alpha:(a,b)->R^n
translated by p."

    rt::usage="t->rt[theta,alpha][t] is the curve alpha:(a,b)->R^2
rotated by theta."

    rf::usage="t->rf[p,alpha][t] is the reflection about the
horizontal axis of the curve alpha:(a,b)->R^2."

(* Focal Set Programs *)

    focalkmu::usage="{u,v}->focalkmu[x][u,v] is a parametrization of the
focal set of a principal patch x corresponding to the
principal curvature kmu."

    focalkpi::usage="{u,v}->focalkpi[x][u,v] is a parametrization of the
focal set of a principal patch x corresponding to the
principal curvature kpi."

    focalk1::usage="{u,v}->focalk1[x][u,v] is a parametrization of the
focal set of a patch x corresponding to the principal curvature k1."

    focalk2::usage="{u,v}->focalk2[x][u,v] is a parametrization of the
focal set of a patch x corresponding to the principal curvature k2."

    tangentkmu::usage="{u,v}->tangentkmu[x][u,v] is an eigenvector of
unit length of a principal patch x corresponding to the
principal curvature kmu."

    tangentkpi::usage="{u,v}->tangentkpi[x][u,v] is an eigenvector of
unit length of a principal patch x corresponding to the
principal curvature kpi."

(* Frenet Frame Programs *)

    darbouxvector::usage="t->darbouxvector[gamma][t] is the Darboux vector
field of the curve gamma:(a,b)->R^3."

    binormal::usage="t->binormal[alpha][t] is the
binormal of the curve alpha:(a,b)->R^3."

    normal::usage="t->normal[alpha][t] is the
normal of the curve alpha:(a,b)->R^3."

    tangent::usage="t->tangent[alpha][t] is the
tangent of the curve alpha:(a,b)->R^3."

(* Geodesics *)

     geo::usage="{u,v}->geo[x][p,q,u,v] gives the system of first order
differential equations for a patch x:U->R^n. geo[ee,ff,gg][p,q,u,v] gives
the system of first order differential equations for the geodesics of a
metric with coefficients ee,ff,gg."

    solvegeoeqs::usage="solvegeoeqs[x,a:0,{u0:0,v0:0},tfin,div,optsnd___]
solves the geodesic equations of the path x using NDSolve."

(* Intrinsic Curvature of Curves *)

   intrinsic::usage="intrinsic[fun,a,{c,d,theta0},optsnd,{smin,smax}] is a
curve alpha starting at {c,d} and initial turning angle theta0 such that
the  signed curvature of alpha is fun."

    plotintrinsic::usage="plotintrinsic[fun,a,{c,d,theta0},optsnd,
{smin,smax}, optspp] plots a curve alpha starting at {c,d} and initial
turning angle theta0 such that the  signed curvature of alpha is fun."

    plotintrinsic3d::usage="plotintrinsic3d[{kk,tt},{a,{p1,p2,p3},{q1,q2,q3},
{r1,r2,r3}},{smin,smax},opts] plots a unit-speed space curve alpha with curvature kk and torsion tt such that alpha[a:0]={p1:0,p2:0,p3:0}, alpha'[a:0]={q1:1,q2:0,q3:0} and alpha''[a:0]=kk[0]*{r1:0,r2:1,r3:0}."

(* Inversion Programs*)

        inversion::usage="x->inversion[p,rho,x] is the inversion of a point
x with respect to a sphere of radius rho centered at p.  inversion[rho,x]
is the inversion of a point x with respect to a sphere of radius rho
centered at the origin."

        inversioncurve::usage="t->inversioncurve[p,rho,alpha][t] is the
inversion of a curve alpha with respect to a sphere of radius rho
centered at p.  t->inversioncurve[rho,alpha][t] is the inversion of a
curve alpha with respect to a sphere of radius rho centered at
the origin."

        inversionsurf::usage="{u,v}->inversionsurf[p,rho,surf][u,v]
is the inversion of a patch  surf with respect to a sphere of radius rho
centered at p.  {u,v}->inversion[rho,surf][u,v] is the inversion of a patch 
surf with respect to a sphere of radius rho centered at the origin."

(* Jacobians *)

    surfacejacobianmatrix::usage="{u,v}->surfacejacobianmatrix[x][u,v]
is the jacobian matrix of the patch x:U->R^n."

    surfacejacobian::usage="{u,v}->surfacejacobian[x][u,v] is
the determinant of the jacobian of the patch x:U->R^n."

(* Metric (First Fundamendal Form) *)

    ee::usage="{u,v}->ee[x][u,v] is the first coefficient of the metric
tensor (or first fundamental form) of the patch x:U->R^n.  The usual
mathematics symbol is E."

    ff::usage="{u,v}->ff[x][u,v] is the second coefficient of the metric
tensor of the patch x:U->R^n  The usual mathematics symbol is F."

    gg::usage="{u,v}->gg[x][u,v] is the third coefficient of the metric
tensor of the patch x:U->R^n.  The usual mathematics symbol is G."

    infarea::usage="{u,v}->infarea[x][u,v] is the infinitesimal area of
the metric tensor (or first fundamental form) of the patch x:U->R^n."

(* Minimal Surface Programs *)

    analytictominimal::usage="{u,v}->analytictominimal[Psi][t][u,v]
is the minimal patch  corresponding to the analytic n-tuple
z->E^(-I*t)*Psi[z]."

     geominsurf::usage="z->geominsurf[alpha][z] is the minimal curve
in C^3 that contains a plane curve alpha:(a,b)->R^3 as a geodesic."

    harmonictoanalytic::usage="z->harmonictoanalytic[x][z] is the
analytic function whose real part is the harmonic n-tuple {u,v}->x[u,v]
and whose imaginary part vanishes at {0,0}.
harmonictoanalytic[x][x0,y0][z] is the analytic function whose real part
is {u,v}->x[u,v] and whose imaginary part vanishes at {x0,y0}."

    wei::usage="z->wei[f,g][z] is the triple of meromorphic functions
constructed from meromorphic functions f and g via the Weierstrass
representation."

    weicurvature::usage="z->weicurvature[f,g][z] is the
Gaussian curvature of the Weierstrass patch determined by
meromorphic functions f and g."

(* Miscellaneous Programs *)

    asymeq1::usage="asymeq1[x][u,v] is the differential equation
eee[x][u,v]u'[v]^2 + 2*fff[x][u,v]u'[v] + ggg[x][u,v]==0 of the
patch x:U->R^3."

    asymeq2::usage="asymeq2[x][u,v] is the differential equation
eee[x][u,v] + 2fff[x][u,v]v'[u] + ggg[x][u,v]v'[u]^2==0 of the
patch x:U->R^3."

    monge::usage="{u,v}->monge[h][u,v] is the parametrization
{u,v}->{u,v,h[u,v]} of the graph of h:U->R^3,
where U is an open subset of R^2."

    nateq::usage="nateq[alpha][k2,s] sometimes can be used to
find the natural equations of a plane curve alpha."

    paramtoimp::usage="paramtoimp[alpha,{x,y}] attempts to convert the
parametric form of a curve alpha to the implicit form.  The result is a
number of equations, one of which is the implicit form of the curve.  The
command will usually not succeed if the components of alpha are not rational
functions."

    polarize::usage="polarize[func][r,theta] is the composition of the
maps {u,v}->func[u,v] and {r,theta}->{r*Cos[theta],r*Sin[theta]},
where func[u,v] is an n-tuple."

    umeqs::usage="umeqs[k][x][u,v] attempts to find the umbilic points of
the patch x:U->R^3."

(* Principal Curvatures *)

    kmu::usage="{u,v}->kmu[x][u,v] is a principal curvature of a
patch x for which ff[x] = fff[x] = 0."

    kpi::usage="{u,v}->kpi[x][u,v] is a principal curvature of a
patch x for which ff[x] = fff[x] = 0."

    k1::usage="{u,v}->k1[x][u,v] is a principal curvature of a
patch x (k1 is greater than or equal to k2)."

    k2::usage="{u,v}->k2[x][u,v] is a principal curvature of a
patch x (k1 is greater than or equal to k2)."

(* Second Fundamental Form and
Gauss Map Programs *)

    eee::usage="{u,v}->eee[x][u,v] is the first coefficient of the
second fundamental form of the patch x:U->R^3. The usual
mathematics symbol is e."

    eeee::usage="{u,v}->eeee[x][u,v] is the first coefficient of the
second fundamental form of the patch x:U->R^3 multiplied by
the infinitesimal area."

    fff::usage="{u,v}->fff[x][u,v] is the second coefficient of the
second fundamental form of the patch x:U->R^3.  The usual
mathematics symbol is f."

    ffff::usage="{u,v}->ffff[x][u,v] is the second coefficient of the
second fundamental form of the patch x:U->R^3 multiplied by
the infinitesimal area."

    ggg::usage="{u,v}->ggg[x][u,v] is the third coefficient of the
second fundamental form of the patch x:U->R^3.  The usual
mathematics symbol is g."

    gggg::usage="{u,v}->gggg[x][u,v] is the third coefficient of the
second fundamental form of the patch x:U->R^3 multiplied by
the infinitesimal area."

    jacobian2::usage="{u,v}->jacobian2[x][u,v] is the square of the
determinant of the jacobian matrix of the patch x:U->R^3."

    unitnormal::usage="{u,v}->unitnormal[x][u,v] is the
unitnormal to the patch x:U->R^3. It is the Gauss map of x."

    weingarten::usage="{u,v}->weingarten[x][u,v] is the matrix of the
shape operator of the patch x:U->R^3 with respect to xu and xv."

(* Vector Analysis Programs *)

    der::usage="der[v][ww][x,y] is the derivative of the vector field
ww on R^2 in the direction v.  der[v][ww][x,y,z] is the derivative
of the vector field ww on R^3 in the direction v."

    divergence::usage="divergence[ww][x,y] is the divergence of a
vector field ww on R^2.  divergence[ww][x,y,z] is the divergence
of the vector field ww on R^3."

    grad::usage="grad[g][x,y] is the gradient of a function g:R^2->R.
grad[g][x,y,z] is the gradient of a function g:R^3->R."

    laplace::usage="laplace[g][x,y] is the Laplacian of a function g:R^2->R.
laplace[g][x,y,z] is the Laplacian of a function g:R^3->R."

Begin["`Private`"]
(* Arc Length Functions *)

    arclengthprime[alpha_][t_]:=
        Sqrt[Simplify[D[alpha[tt],tt].
                D[alpha[tt],tt]]] /. tt->t

    arclength[alpha_][t_]:=
        Integrate[arclengthprime[alpha][tt],tt]  /. tt->t

    length[a_,b_][alpha_]:=
        Integrate[arclengthprime[alpha][u],{u,a,b}]

    nlength[a_,b_][alpha_]:=
        NIntegrate[arclengthprime[alpha][u],{u,a,b}]

    ndslength[a_,b_][alpha_][t_]:=
        Module[{aaa,bbb,s,tt},
            aaa=NDSolve[{s'[tt]==
                arclengthprime[alpha][tt],
                s[a]==0},s,{tt,a,b}];
            bbb = s[tt] /. aaa;
            bbb /. tt->t][[1]]//N

(* Brioschi Curvature *)
             
    brioschicurvature[x_][u_,v_]:=
        (Det[{
           {-D[ee[x][uu,vv],{vv,2}]/2 + 
             D[ff[x][uu,vv],uu,vv] -
             D[gg[x][uu,vv],{uu,2}]/2,
             D[ee[x][uu,vv],uu]/2,
             D[ff[x][uu,vv],uu] -
             D[ee[x][uu,vv],vv]/2},
           {D[ff[x][uu,vv],vv] -
            D[gg[x][uu,vv],uu]/2,
            ee[x][uu,vv],ff[x][uu,vv]},
          {D[gg[x][uu,vv],vv]/2,
          ff[x][uu,vv],gg[x][uu,vv]}}]
       -Det[{
          {0,D[ee[x][uu,vv],vv]/2,
             D[gg[x][uu,vv],uu]/2},
          {D[ee[x][uu,vv],vv]/2,
             ee[x][uu,vv],ff[x][uu,vv]},
          {D[gg[x][uu,vv],uu]/2,
             ff[x][uu,vv],gg[x][uu,vv]}}])/
       (ee[x][uu,vv]*gg[x][uu,vv] - ff[x][uu,vv]^2)^2 /.
                  {uu->u,vv->v}

    brioschicurvature[ee_,ff_,gg_][u_,v_]:=
        (Det[{{-D[ee[uu,vv],{vv,2}]/2 + D[ff[uu,vv],uu,vv] -
             D[gg[uu,vv],{uu,2}]/2,
             D[ee[uu,vv],uu]/2,
             D[ff[uu,vv],uu] - D[ee[uu,vv],vv]/2},
        {D[ff[uu,vv],vv] - D[gg[uu,vv],uu]/2,ee[uu,vv],ff[uu,vv]},
        {D[gg[uu,vv],vv]/2,ff[uu,vv],gg[uu,vv]}}] -
         Det[{{0,D[ee[uu,vv],vv]/2,D[gg[uu,vv],uu]/2},
            {D[ee[uu,vv],vv]/2,ee[uu,vv],ff[uu,vv]},
            {D[gg[uu,vv],uu]/2,ff[uu,vv],gg[uu,vv]}}])/
        (ee[uu,vv]*gg[uu,vv] - ff[uu,vv]^2)^2   /. {uu->u,vv->v}

(* Christoffel Symbols *)

     christoffel[1,1,1][x_][u_,v_]:=
        Simplify[(D[ee[x][uu,vv],vv]*ff[x][uu,vv] -
        2*ff[x][uu,vv]*D[ff[x][uu,vv],uu] +
        D[ee[x][uu,vv],uu]*gg[x][uu,vv])/
       (2*(ee[x][uu,vv]*gg[x][uu,vv] -
        ff[x][uu,vv]^2))] /. {uu->u,vv->v}

     christoffel[2,1,1][x_][u_,v_]:=
        Simplify[-(D[ee[x][uu,vv],vv]*ee[x][uu,vv] -
        2*ee[x][uu,vv]*D[ff[x][uu,vv],uu] +
        D[ee[x][uu,vv],uu]*ff[x][uu,vv])/
       (2*(ee[x][uu,vv]*gg[x][uu,vv] -
        ff[x][uu,vv]^2))] /. {uu->u,vv->v}
                        
    christoffel[1,1,2][x_][u_,v_]:=
        Simplify[(D[ee[x][uu,vv],vv]*gg[x][uu,vv] -
        ff[x][uu,vv]*D[gg[x][uu,vv],uu])/
        (2*(ee[x][uu,vv]*gg[x][uu,vv] -
        ff[x][uu,vv]^2))] /. {uu->u,vv->v}

    christoffel[1,2,1][x_][u_,v_]:=
        christoffel[1,1,2][x][u,v]
                        
    christoffel[2,1,2][x_][u_,v_]:=
        Simplify[-(D[ee[x][uu,vv],vv]*ff[x][uu,vv] -
         ee[x][uu,vv]*D[gg[x][uu,vv],uu])/
        (2*(ee[x][uu,vv]*gg[x][uu,vv] -
        ff[x][uu,vv]^2))] /. {uu->u,vv->v}

    christoffel[2,2,1][x_][u_,v_]:=
         christoffel[2,1,2][x][u,v]
                        
    christoffel[1,2,2][x_][u_,v_]:=
        Simplify[(-D[gg[x][uu,vv],vv]*ff[x][uu,vv] +
        2*gg[x][uu,vv]*D[ff[x][uu,vv],vv] -
        D[gg[x][uu,vv],uu]*gg[x][uu,vv])/
        (2*(ee[x][uu,vv]*gg[x][uu,vv] -
        ff[x][uu,vv]^2))] /. {uu->u,vv->v}

    christoffel[2,2,2][x_][u_,v_]:=
        Simplify[(D[gg[x][uu,vv],vv]*ee[x][uu,vv] -
        2*ff[x][uu,vv]*D[ff[x][uu,vv],vv] +
        D[gg[x][uu,vv],uu]*ff[x][uu,vv])/
        (2*(ee[x][uu,vv]*gg[x][uu,vv] -
         ff[x][uu,vv]^2))] /. {uu->u,vv->v}

    christoffel[1,1,1][ee_,ff_,gg_][u_,v_]:=
         Simplify[(D[ee[u,v],v]*ff[u,v] -
         2*ff[u,v]*D[ff[u,v],u] +
         D[ee[u,v],u]*gg[u,v])/
        (2*(ee[u,v]*gg[u,v] - ff[u,v]^2))]

    christoffel[2,1,1][ee_,ff_,gg_][u_,v_]:=
         Simplify[-(D[ee[u,v],v]*ee[u,v] -
         2ee[u,v]*D[ff[u,v],u] +
         D[ee[u,v],u]*ff[u,v])/
        (2*(ee[u,v]*gg[u,v] - ff[u,v]^2))]

    christoffel[1,1,2][ee_,ff_,gg_][u_,v_]:=
         Simplify[(D[ee[u,v],v]*gg[u,v] -
         ff[u,v]*D[gg[u,v],u])/
         (2*(ee[u,v]*gg[u,v] - ff[u,v]^2))]

    christoffel[1,2,1][ee_,ff_,gg_][u_,v_]:=
         christoffel[1,1,2][ee,ff,gg][u,v]

    christoffel[2,1,2][ee_,ff_,gg_][u_,v_]:=
         Simplify[-(D[ee[u,v],v]*ff[u,v] -
          ee[u,v]*D[gg[u,v],u])/
         (2*(ee[u,v]*gg[u,v] - ff[u,v]^2))]

    christoffel[2,2,1][ee_,ff_,gg_][u_,v_]:=
          christoffel[2,1,2][ee,ff,gg][u,v]

    christoffel[1,2,2][ee_,ff_,gg_][u_,v_]:=
         Simplify[(-D[gg[u,v],v]*ff[u,v] +
         2*gg[u,v]*D[ff[u,v],v] -
         D[gg[u,v],u]*gg[u,v])/
         (2*(ee[u,v]*gg[u,v] - ff[u,v]^2))]

     christoffel[2,2,2][ee_,ff_,gg_][u_,v_]:=
         Simplify[(D[gg[u,v],v]*ee[u,v] -
         2*ff[u,v]*D[ff[u,v],v] +
         D[gg[u,v],u]*ff[u,v])/
         (2*(ee[u,v]*gg[u,v] - ff[u,v]^2))]

(* Complex to Real to Complex\label{ComplextoRealtoComplex *)

    ctor[z_]:= ComplexExpand[{Re[z],Im[z]},{z}]

    rtoc[{x_,y_}]:= x + I*y

    cdot[z_,w_]:= ComplexExpand[Re[z.Conjugate[w]],{z,w}]

(* Complex Structure *)

    J[{p1_,p2_}]:= {-p2,p1}

(* Cross Product *)

    Cross[u_List,v_List]:=
        RotateLeft[u*RotateLeft[v] -
        RotateLeft[u]*v] /;
        Length[u] == Length[v] == 3

(* Construction of Plane Curves from
Plane Curves  *)

    circle3point[{p1_,q1_},{p2_,q2_},{p3_,q3_}]:=
      Circle[{(p2^2*q1 - p3^2*q1 - p1^2*q2 + p3^2*q2 -
      q1^2*q2 + q1*q2^2 + p1^2*q3 - p2^2*q3 + q1^2*q3 -
      q2^2*q3 - q1*q3^2 + q2*q3^2)/
      (2*(p2*q1 - p3*q1 - p1*q2 + p3*q2 + p1*q3 - p2*q3)), 
      (p1^2*p2 - p1*p2^2 - p1^2*p3 + p2^2*p3 + p1*p3^2 - 
      p2*p3^2 + p2*q1^2 - p3*q1^2 - p1*q2^2 + p3*q2^2 + 
      p1*q3^2 - p2*q3^2)/
      (2*(p2*q1 - p3*q1 - p1*q2 + p3*q2 + p1*q3 - p2*q3))},
      Abs[((p1^2 - 2*p1*p2 + p2^2 + q1^2 -
      2*q1*q2 + q2^2)^(1/2)*(p1^2 - 2*p1*p3 + p3^2 +
      q1^2 - 2*q1*q3 + q3^2)^(1/2)*(p2^2 - 2*p2*p3 +
      p3^2 + q2^2 - 2*q2*q3 + q3^2)^(1/2))/(2*(-(p2*q1) +
      p3*q1 + p1*q2 - p3*q2 - p1*q3 + p2*q3))]]
                
    evolute[alpha_][t_]:= alpha[tt] +
        D[alpha[tt],tt].D[alpha[tt],tt]/
       (D[alpha[tt],{tt,2}].
        J[D[alpha[tt],tt]])*J[D[alpha[tt],tt]] /. {tt->t}

    evolute[1,alpha_][t_]:= Simplify[alpha[tt] +
        D[alpha[tt],tt].D[alpha[tt],tt]/
        (D[alpha[tt],{tt,2}].J[D[alpha[tt],tt]])*
         J[D[alpha[tt],tt]] /. {tt->t}]

    evolute[n_,alpha_][t_]:= evolute[n,alpha][t] =
        evolute[1,evolute[n - 1,alpha]][t]//Simplify

    involute[alpha_,c_][t_]:= alpha[tt] +
                        (arclength[alpha][c] -
                         arclength[alpha][tt])*D[alpha[tt],tt]/
                         Sqrt[Factor[D[alpha[tt],tt].
                         D[alpha[tt],tt]]] /. tt->t      

    involute[alpha_][t_]:= involute[alpha,0][t]

    ndsinvolute[a_,b_][alpha_][t_]:=
        Module[{aaa,bbb,s,tt},
            aaa=NDSolve[{s'[tt]==
                arclengthprime[alpha][tt],
                s[a]==0},s,{tt,a,b}];
            bbb=alpha[tt] +
                (s[a] - s[tt])*D[alpha[tt],tt]/
                 Sqrt[D[alpha[tt],tt].
                 D[alpha[tt],tt]]  /.aaa;
            bbb /. tt->t]
                
    normalline[alpha_][s_,t_]:=
        Line[{alpha[tt],alpha[tt] + s*J[D[alpha[tt],tt]]/
        (D[alpha[tt],tt].D[alpha[tt],tt])^(1/2)}] /. tt->t
                
    osculatingcircle[alpha_][t_]:=
        Circle[alpha[tt] + alpha'[tt].D[alpha[tt],tt]/
        (D[alpha[tt],{tt,2}].
        J[D[alpha[tt],tt]])*J[D[alpha[tt],tt]],
        (D[alpha[tt],tt].D[alpha[tt],tt])^(3/2)/
        (D[alpha[tt],{tt,2}].J[D[alpha[tt],
                tt]])] /. tt->t

    tangentline[alpha_][s_,t_]:= Line[{alpha[tt],alpha[tt] +
          s*D[alpha[tt],tt]/Sqrt[Factor[D[alpha[tt],tt].
                  D[alpha[tt],tt]]]}] /. tt->t

    plottl[curve_,t0_,{tmin_,tmax_},opts___]:=
          ParametricPlot[curve[t0 + t]//Evaluate,{t,tmin,tmax},
          Epilog->{tangentline[curve][t0][tmax]//Evaluate},
          AspectRatio->Automatic,opts];
                
    parcurve[alpha_][s_][t_]:= alpha[tt] +
          s*J[D[alpha[tt],tt]]/Sqrt[Simplify[
              D[alpha[tt],tt].D[alpha[tt],tt]]] /. tt->t
                
    pedal[{a_,b_},alpha_][t_]:= {a,b} +
        ((alpha[tt] - {a,b}).J[D[alpha[tt],tt]]/
        (D[alpha[tt],tt].D[alpha[tt],tt]))*J[
         D[alpha[tt],tt]] /. {tt->t}

    pedal[alpha_][t_]:= pedal[{0,0},alpha][t]

(* Construction of Space Curves from
Plane Curves  *)

     helical[alpha_,a_][t_]:=
              {alpha[t][[1]],alpha[t][[2]],a*t}

     helical[alpha_][f_][t_]:=
              {alpha[t][[1]],alpha[t][[2]],f[t]}

     inscurv[n_,a_,alpha_][t_]:= Insert[alpha[t],a,n]
              
     kappa23d[alpha_][t_]:= Append[alpha[t],kappa2[alpha][t]]

     kappa23d[mag_,alpha_][t_]:=
           Append[alpha[t],mag*kappa2[alpha][t]]

     writhe[alpha_][t_]:= {alpha[t][[1]],
            Cos[t]*alpha[t][[2]],Sin[t]*alpha[t][[2]]}

(* Construction of Surfaces from
Plane Curves  *)

    cone[{p1_,p2_,p3_},gamma_][u_,v_]:=
        {p1,p2,p3} + v*gamma[u]

    cylinder[{q1_,q2_,q3_},gamma_][u_,v_]:=
        gamma[u] + v*{q1,q2,q3}

    genhel[alpha_][sl_][u_,v_]:= {Cos[u]*alpha[v][[1]],
                                  Sin[u]*alpha[v][[1]],
                                  sl*u + alpha[v][[2]]}
     
    gensurfrev[alpha_,gamma_][u_,v_]:= {alpha[u][[1]]*gamma[v][[1]],
                   alpha[u][[2]]*gamma[v][[1]],gamma[v][[2]]}

    ruled[alpha_,gamma_][u_,v_]:= alpha[u] + v*gamma[u]

    surfast[alpha_][u_,v_]:= {Cos[u]^3*alpha[v][[1]],
                              Sin[u]^3*alpha[v][[1]],
                                       alpha[v][[2]]}

    surfast[n_,alpha_][u_,v_]:= {Cos[u]^n*alpha[v][[1]],
                                 Sin[u]^n*alpha[v][[1]],
                                          alpha[v][[2]]}
 
    surfrev[alpha_][u_,v_]:= {Cos[u]*alpha[v][[1]],
        Sin[u]*alpha[v][[1]],alpha[v][[2]]}
 
    twist[a_,b_,gamma_][u_,v_]:=
        (a + Cos[b*u]*gamma[v][[1]] - Sin[b*u]*gamma[v][[2]])*
        {Cos[u],Sin[u],0} +
        (Sin[b*u]*gamma[v][[1]] + Cos[b*u]*gamma[v][[2]])*
        {0,0,1}

(* Construction of Surfaces from
Space Curves  *)

    binormaldev[alpha_][u_,v_]:=
        alpha[u] + v*binormal[alpha][u]
         
    normaldev[alpha_][u_,v_]:=
        alpha[u] + v*normal[alpha][u]

    perpdev[phi_,alpha_][u_, v_]:= 
        alpha[u] + v*(Cos[phi]*normal[alpha][u] +
                Sin[phi]*binormal[alpha][u])

    seashell[gamma_][r_][t_,theta]:= gamma[t] + 
        r*t*(Cos[theta]*normal[gamma][t] +
          Sin[theta]*binormal[gamma][t])

    tandev[alpha_][u_,v_]:= alpha[uu] + v*D[alpha[uu],uu]/
          Sqrt[D[alpha[uu],uu].D[alpha[uu],uu]] /. uu->u
                
    tubecurve[gamma_][r_][t_,theta_]:= Simplify[gamma[t] + 
        r*(Cos[theta]*normal[gamma][t] +
          Sin[theta]*binormal[gamma][t])]

    vartube[gamma_][r_][t_,theta_]:= Simplify[gamma[t] + 
        r[t]*(Cos[theta]*normal[gamma][t] +
             Sin[theta]*binormal[gamma][t])]

(* Construction of Surfaces from
Surfaces  *)
                
    parsurf[x_][r_][u_,v_]:= Module[{xu,xv,n1,n2,n3,n4},
        xu = D[x[uu,vv],uu];
        xv = D[x[uu,vv],vv];
        n1 = Cross[xu,xv]//Simplify;
        n2 = Simplify[Factor[n1.n1]];
        n3 = PowerExpand[n2^(1/2)]//Simplify;
        n4 = Simplify[n1/n3];
        Simplify[Together[x[uu,vv] +
                r*n4]]] /. {uu->u,vv->v}
                
    pedsurf[{a_,b_,c_},x_][u_,v_]:=
         Module[{xu,xv,n1,s2,s3,n4},
         xu = D[x[uu,vv],uu];
        xv = D[x[uu,vv],vv];
        n1 = Cross[xu, xv]//Simplify;
        s2 = Simplify[Factor[n1.n1]];
        s3 = PowerExpand[s2^(1/2)]//Simplify;
        n4 = Simplify[n1/s3];
        {a,b,c} + Simplify[Together[((x[uu,vv] -
        {a,b,c}).n4)*n4]]] /. {uu->u,vv->v}
           
     pedsurf[x_][u_,v_]:= Module[{xu,xv,n1,s2,s3,n4},
        xu = D[x[uu,vv],uu];
        xv = D[x[uu,vv],vv];
        n1 = Cross[xu, xv]//Simplify;
        s2 = Simplify[Factor[n1.n1]];
        s3 = PowerExpand[s2^(1/2)]//Simplify;
        n4 = Simplify[n1/s3];
        Simplify[Together[(x[uu,vv].n4)*n4]]]  /.
             {uu->u,vv->v}

    stretch[a_,b_,c_,surf_][u_,v_]:=
           {a*surf[u,v][[1]],b*surf[u,v][[2]],c*surf[u,v][[3]]}
    
    tangentplane[x_][u0_,v0_][u_,v_]:= x[u0,v0] +
          D[x[uu0,vv0],uu0]*u + D[x[uu0,vv0],vv0]*v /.
          {uu0->u0,vv0->v0}

    plottp[surf_,{u0_,v0_},{umin_,umax_},
          {vmin_,vmax_},opts___]:=
          ParametricPlot3D[{surf[u0 + u,v0 + v],
             tangentplane[surf][u0,v0][u,v]}//Evaluate,
             {u,umin,umax},{v,vmin,vmax},
             Boxed->False,Axes->None,opts];

 (* Curvature of Surfaces (Fast) *)

     gcurvature[x_][u_,v_]:= Simplify[
       (Det[{D[x[uu,vv],uu,uu],D[x[uu,vv],uu],
             D[x[uu,vv],vv]}]*
        Det[{D[x[uu,vv],vv,vv],D[x[uu,vv],uu],
             D[x[uu,vv],vv]}] -
        Det[{D[x[uu,vv],uu,vv],D[x[uu,vv],uu],
             D[x[uu,vv],vv]}]^2)/
       (D[x[uu,vv],uu].D[x[uu,vv],uu]*
             D[x[uu,vv],vv].D[x[uu,vv],vv] -
        D[x[uu,vv],uu].
             D[x[uu,vv],vv]^2)^2] /. {uu->u,vv->v}
 
     mcurvature[x_][u_,v_]:= Simplify[
       (Det[{D[x[uu,vv],uu,uu],D[x[uu,vv],uu],
            D[x[uu,vv],vv]}]*
            D[x[uu,vv],vv].D[x[uu,vv],vv] -
       2*Det[{D[x[uu,vv],uu,vv],D[x[uu,vv],uu],
            D[x[uu,vv],vv]}]*
            D[x[uu,vv],uu].D[x[uu,vv],vv] +
        Det[{D[x[uu,vv],vv,vv],D[x[uu,vv],uu],
            D[x[uu,vv],vv]}]*
            D[x[uu,vv],uu].D[x[uu,vv],uu])/
      (2*(D[x[uu,vv],uu].D[x[uu,vv],uu]*
            D[x[uu,vv],vv].D[x[uu,vv],vv] -
            D[x[uu,vv],uu].
            D[x[uu,vv],vv]^2)^(3/2))] /. {uu->u,vv->v}

(* Curvature of Surfaces (Slow) *)

     gaussiancurvature[x_][u_,v_]:= Simplify[
        (eee[x][uu,vv]*ggg[x][uu,vv] - fff[x][uu,vv]^2)/
        (ee[x][uu,vv]*gg[x][uu,vv] -
                ff[x][uu,vv]^2)] /. {uu->u,vv->v}

     meancurvature[x_][u_,v_]:= Simplify[
        (eee[x][uu,vv]*gg[x][uu,vv] -
                2*fff[x][uu,vv]*ff[x][uu,vv] +
         ggg[x][uu,vv]*ee[x][uu,vv])/
        (2*(ee[x][uu,vv]*gg[x][uu,vv] -
                ff[x][uu,vv]^2))] /. {uu->u,vv->v}

(* Curvature and Torsion of Curves *)
         
    kappa[alpha_][t_]:=
        Simplify[Factor[Cross[D[alpha[tt],tt],
        D[alpha[tt],{tt,2}]].
        Cross[D[alpha[tt],tt],D[alpha[tt],{tt,2}]]]]^(1/2)/
        Simplify[Factor[D[alpha[tt],tt].
        D[alpha[tt],tt]]]^(3/2) /. tt->t

    kappa2[alpha_][t_]:=
        D[alpha[tt],{tt,2}].J[D[alpha[tt],tt]]/
         Simplify[D[alpha[tt],tt].
              D[alpha[tt],tt]]^(3/2) /. tt->t

    kappa2c[alpha_][t_]:= Simplify[
        ComplexExpand[Re[D[alpha[t],{t,2}]*
        Conjugate[I*D[alpha[t],t]]]]/
        ComplexExpand[Abs[D[alpha[t],t]],
        TargetFunctions->{Re,Im}]^3]

    kappa2implicit[g_][x_,y_]:=
        (D[g[xx,yy],{xx,2}]*D[g[xx,yy],yy]^2 -
         2D[g[xx,yy],xx,yy]*D[g[xx,yy],xx]*D[g[xx,yy],yy] +
         D[g[xx,yy],{yy,2}]*D[g[xx,yy],xx]^2)/
        (D[g[xx,yy],xx]^2 + D[g[xx,yy],yy]^2)^(3/2) /. 
                         {xx->x,yy->y}

    kappan[alpha_][t_]:=
        Simplify[Sqrt[Factor[
        (Sqrt[Factor[D[alpha[tt],tt].
                     D[alpha[tt],tt]]]*D[alpha[tt],{tt,2}] -
        D[Sqrt[Factor[D[alpha[tt],tt].
                     D[alpha[tt],tt]]],tt]*D[alpha[tt],tt]).
        (Sqrt[Factor[D[alpha[tt],tt].
                     D[alpha[tt],tt]]]*D[alpha[tt],{tt,2}] -
        D[Sqrt[Factor[D[alpha[tt],tt].
                     D[alpha[tt],tt]]],tt]*
        D[alpha[tt],tt])]]]/
        Simplify[Factor[D[alpha[tt],tt].
                     D[alpha[tt],tt]]]^(3/2) /. tt->t

    kappa2polar[rho_][t_]:= (-D[rho[tt],{tt,2}]*rho[tt]
          + 2*D[rho[tt],tt]^2 + rho[tt]^2)/
         (D[rho[tt],tt]^2 + rho[tt]^2)^(3/2) /. tt->t

    tau[alpha_][t_]:= Simplify[Det[
        {D[alpha[tt],tt],D[alpha[tt],{tt,2}],
        D[alpha[tt],{tt,3}]}]]/
        Simplify[Factor[Cross[D[alpha[tt],tt],
        D[alpha[tt],{tt,2}]].Cross[D[alpha[tt],tt],
        D[alpha[tt],{tt,2}]]]] /. tt->t
              
    totalacsalrcurvature[a_,b_,alpha_]:= Integrate[(kappa2[alpha][t]*
        Sqrt[Simplify[D[alpha[t],u].
        D[alpha[t],u]]])//PowerExpand,{u,a,b}]

    turningangle[alpha_][t_]:= ArcTan[D[alpha[tt],tt][[1]],
               D[alpha[tt],tt][[2]]] /. tt->t
                
    turningnumber[alpha_,{a_,b_},opts___]:=
        Module[{tt,tmp1,tmp2},
            tmp1=kappa2[alpha][tt]*Sqrt[alpha'[tt].
                alpha'[tt]];
            tmp2=NIntegrate[tmp1,{tt,a,b},
                AccuracyGoal->2,opts];
            Round[tmp2/N[2Pi]]]

(* Euclidean Motions Acting on
Plane Curves *)

    tl[p_,alpha_][t_]:= p + alpha[t]

    rt[theta_,alpha_][t_]:=
       {Cos[theta]*alpha[t][[1]] - Sin[theta]*alpha[t][[2]],
        Sin[theta]*alpha[t][[1]] + Cos[theta]*alpha[t][[2]]}

    rf[alpha_][t_]:= {alpha[t][[1]],-alpha[t][[2]]}

(* Focal Set Programs *)

    focalkmu[surf_][u_,v_]:= (surf[uu,vv] +
        (1/kmu[surf][uu,vv])*
        unitnormal[surf][uu,vv] /. {uu->u,vv->v})

    focalkpi[surf_][u_,v_]:= (surf[uu,vv] +
        (1/kpi[surf][uu,vv])*
        unitnormal[surf][uu,vv] /. {uu->u,vv->v}) 

    focalk1[surf_][u_,v_]:= surf[uu,vv] +
        (1/k1[surf][uu,vv])*
        unitnormal[surf][uu,vv] /. {uu->u,vv->v}

    focalk2[surf_][u_,v_]:= surf[uu,vv] +
        (1/k2[surf][uu,vv])*
        unitnormal[surf][uu,vv] /. {uu->u,vv->v}

    tangentkmu[x_][u_,v_]:= (D[x[u,vv],vv]/
                Sqrt[gg[x][u,vv]] /. vv->v) /;
                   ff[x][u,v]==fff[x][u,v]==0

    tangentkpi[x_][u_,v_]:= (D[x[uu,v],uu]/
                Sqrt[gg[x][uu,v]] /. uu->u) /;
                   ff[x][u,v]==fff[x][u,v]==0
                
(* Frenet Frame Programs *)

    binormal[alpha_][t_]:= Simplify[Cross[
        D[alpha[tt],tt],D[alpha[tt],{tt,2}]]]/
        Simplify[Factor[Cross[D[alpha[tt],tt],
        D[alpha[tt],{tt,2}]].
        Cross[D[alpha[tt],tt],D[alpha[tt],
        {tt,2}]]]]^(1/2)  /. tt->t
 
    darbouxvector[gamma_][t_]:=
         Sqrt[D[gamma[tt],tt].D[gamma[tt],tt]]*(
                tau[gamma][tt]*tangent[gamma][tt] +
                kappa[gamma][tt]*binormal[gamma][tt]) /. tt->t

    normal[alpha_][t_]:= Cross[binormal[alpha][t],
                           tangent[alpha][t]]

    tangent[alpha_][t_]:= D[alpha[tt],tt]/
        Simplify[Factor[D[alpha[tt],tt].
          D[alpha[tt],tt]]]^(1/2)  /. tt->t

(* Geodesics *)

    geo[x_][p_,q_,u_,v_]:=
        {p,q,-christoffel[1,1,1][x][u,v]*p^2 -
        2*christoffel[1,1,2][x][u,v]*p*q -
         christoffel[1,2,2][x][u,v]*q^2,
        -christoffel[2,1,1][x][u,v]*p^2 -
        2*christoffel[2,1,2][x][u,v]*p*q -
         christoffel[2,2,2][x][u,v]*q^2}
   
     geo[ee_,ff_,gg_][p_,q_,u_,v_]:=
         {p,q,-christoffel[1,1,1][ee,ff,gg][u,v]*p^2 -
         2*christoffel[1,1,2][ee,ff,gg][u,v]*p*q -
          christoffel[1,2,2][ee,ff,gg][u,v]*q^2,
         -christoffel[2,1,1][ee,ff,gg][u,v]*p^2 -
         2*christoffel[2,1,2][ee,ff,gg][u,v]*p*q -
          christoffel[2,2,2][ee,ff,gg][u,v]*q^2}

    solvegeoeqs[x_,a_:0,{u0_:0,v0_:0},tfin_,div_,optsnd___]:=
        Flatten[Map[NDSolve[#,{u,v},{ss,0,tfin}]&,
        Table[Join[MapThread[Equal,
        {{u''[ss],v''[ss]},
        Take[geo[x][u'[ss],v'[ss],u[ss],v[ss]],{3,4}]}],
        {u[a]==u0,v[a]==v0,
        u'[a]==Cos[theta],v'[a]==Sin[theta]}],
        {theta,2Pi/div,2Pi,2Pi/div}]],1]

(* Intrinsic Curvature of Curves *)

     intrinsic[fun_,a_:0,{c_:0,d_:0,theta0_:0},optsnd___,
        {smin_:-10,smax_:10}][t_]:=
        Flatten[Module[{x,y,theta},
        {x[t],y[t]} /.
        NDSolve[{x'[ss]==Cos[theta[ss]],
             y'[ss]==Sin[theta[ss]],
             theta'[ss]==fun[ss],
             x[a]==c,y[a]==d,theta[a]==theta0},
             {x,y,theta},{ss,smin,smax},optsnd]]]

     plotintrinsic[fun_,a_:0,{c_:0,d_:0,theta0_:0},optsnd___,
                  {smin_:-10,smax_:10},optspp___]:=
          ParametricPlot[Module[{x,y,theta},
          {x[t],y[t]} /.
          NDSolve[{x'[ss]==Cos[theta[ss]],
             y'[ss]==Sin[theta[ss]],
             theta'[ss]==fun[ss],
             x[a]==c,y[a]==d,theta[a]==theta0},
             {x,y,theta},{ss,smin,smax},optsnd]]//Evaluate,
          {t,smin,smax},AspectRatio->Automatic,optspp];

    plotintrinsic[fun_,optsnd___,
       {smin_:-10,smax_:10},optspp___]:=
       plotintrinsic[fun,0,{0,0,0},{smin,smax},optspp]

    plotintrinsic3d[{kk_,tt_},{a_:0,{p1_:0,p2_:0,p3_:0},
           {q1_:1,q2_:0,q3_:0},{r1_:0,r2_:1,r3_:0}},
                  {smin_:-10,smax_:10},opts___]:=
        ParametricPlot3D[Module[
           {x1,x2,x3,t1,t2,t3,n1,n2,n3,b1,b2,b3},
          {x1[s],x2[s],x3[s]} /.
          NDSolve[{x1'[ss]==t1[ss],
                 x2'[ss]==t2[ss],
                 x3'[ss]==t3[ss],
             t1'[ss]==kk[ss]*n1[ss],
             t2'[ss]==kk[ss]*n2[ss],
             t3'[ss]==kk[ss]*n3[ss],
                 n1'[ss]==-kk[ss]*t1[ss] + tt[ss]*b1[ss],
                 n2'[ss]==-kk[ss]*t2[ss] + tt[ss]*b2[ss],
                 n3'[ss]==-kk[ss]*t3[ss] + tt[ss]*b3[ss],
                     b1'[ss]==-tt[ss]*n1[ss],
                     b2'[ss]==-tt[ss]*n2[ss],
                     b3'[ss]==-tt[ss]*n2[ss],
             x1[a]==p1,x2[a]==p2,x3[a]==p3,
             t1[a]==q1,t2[a]==q2,t3[a]==q3,
             n1[a]==r1,n2[a]==r2,n3[a]==r3,
                 b1[a]==q2*r3 - q3*r2,
                 b2[a]==q3*r1 - q1*r3,
                 b3[a]==q1*r2 - q2*r1},
           {x1,x2,x3,t1,t2,t3,n1,n2,n3,b1,b2,b3},
             {ss,smin,smax}]]//Evaluate,
             {s,smin,smax},opts];

(* Inversion Programs *)
 
    inversion[p_,rho_,x_]:= p + (rho^2/((x - p).(x - p)))*(x - p)

    inversion[rho_,x_]:= (rho^2/x.x)*x

    inversioncurve[p_List,rho_,curve_][t_]:=
                               inversion[p,rho,curve[t]]

    inversioncurve[rho_,curve_][t_]:=
                        inversion[rho,curve[t]]

    inversionsurf[p_List,rho_,surf_][u_,v_]:=
                        inversion[p,rho,surf[u,v]]

    inversionsurf[rho_,surf_][u_,v_]:=
                        inversion[rho,surf[u,v]]

(* Jacobians *)


    surfacejacobianmatrix[x_][u_,v_]:= {D[x[uu,vv],uu],
          D[x[uu,vv],vv]}  /. {uu->u,vv->v}

    surfacejacobian[x_][u_,v_]:= Det[
        surfacejacobianmatrix[x][u,v].
        Transpose[surfacejacobianmatrix[x][u,v]]]

(* Metric (First Fundamendal Form) *)

    ee[x_][u_,v_]:= Simplify[D[x[uu,vv],uu].
                        D[x[uu,vv],uu]] /. {uu->u,vv->v}

    ff[x_][u_,v_]:= Simplify[D[x[uu,vv],uu].
                        D[x[uu,vv],vv]] /. {uu->u,vv->v}

    gg[x_][u_,v_]:= Simplify[D[x[uu,vv],vv].
                        D[x[uu,vv],vv]] /. {uu->u,vv->v}

    infarea[x_][u_,v_]:= Sqrt[Simplify[ee[x][u,v]*gg[x][u,v] -
                        ff[x][u,v]^2]]

(* Minimal Surface Programs *)

   analytictominimal[Psi_][t_][u_,v_]:=
        ComplexExpand[Re[E^(-I*t)*Psi[zz]]
        /. zz->u + I*v]

     geominsurf[alpha_][z_]:= Append[alpha[z],
            I*Simplify[PowerExpand[arclength[alpha][z]]]]

     harmonictoanalytic[x_][z_]:=
        Simplify[2x[z/2,z/(2*I)] - x[0,0]]

     harmonictoanalytic[x_][x0_,y0_][z_]:=
        Simplify[2x[(z + x0 - I*y0)/2,
        (z - x0 + I*y0)/(2*I)] - x[x0,y0]]
        
    wei[f_,g_][z_]:= Module[{zz,tmp},
           tmp=Integrate[{(1 - g[zz]^2)*f[zz]/2,
              I*(1 + g[zz]^2)*f[zz]/2,f[zz]g[zz]},zz];
           Simplify[tmp /. zz->z]]

    weicurvature[f_,g_][z_]:=
        -4*Abs[D[g[zz],zz]]/(Abs[f[z]]*(1 + g[zz]^2)^4)

(* Miscellaneous Programs *)
              
    asymeq1[x_][u_,v_]:=
        Module[{uuu,a,xxx,eqa1,eqa2},
        xxx=Solve[eeee[x][uuu,v]*a^2 +
        2*ffff[x][uuu,v]*a +
        gggg[x][uuu,v]==0,a];
        eqa1=u'[v]== a /. xxx[[1]];
        eqa2=u'[v]== a /. xxx[[2]];
        {eqa1,eqa2} /. {a->u'[v],uuu->u[v]}]

    asymeq2[x_][u_,v_]:=
        Module[{vvv,a,xxx,eqa1,eqa2},
        xxx=Solve[eeee[x][u,vvv] +
        2*ffff[x][u,vvv]*b +
        gggg[x][u,vvv]*b^2==0,b];
        eqa1=v'[u]== b /. xxx[[1]];
        eqa2=v'[u]== b /. xxx[[2]];
        {eqa1,eqa2} /. {b->v'[u],vvv->v[u]}]

     monge[h_][u_,v_]:= {u,v,h[u,v]}
        
     nateq[alpha_][k2_,s_]:= Eliminate[
        {s==arclength[alpha][t]//PowerExpand,
         k2==kappa2[alpha][t]//PowerExpand//Simplify},t]

     nateq[alpha_][k2_,s_]:= ToRules[Eliminate[
        {s==arclength[alpha][t]//PowerExpand,
         k2==kappa2[alpha][t]//PowerExpand//Simplify},t]]

     paramtoimp[alpha_,{x_,y_}]:=
            Eliminate[{x,y}==alpha[t],t]
           
     polarize[func_][r_,theta_]:=
        func[r*Cos[theta],r*Sin[theta]]

     umeqs[k_][x_][u_,v_]:=
       {ee[x][u,v],ff[x][u,v],gg[x][u,v]} -
       k*{eee[x][u,v],fff[x][u,v],ggg[x][u,v]}


(* Principal Curvatures *)

     kmu[x_][u_,v_]:= ggg[x][u,v]/gg[x][u,v] /;
                   ff[x][u,v]==fff[x][u,v]==0

     kpi[x_][u_,v_]:= eee[x][u,v]/ee[x][u,v] /;
                   ff[x][u,v]==fff[x][u,v]==0

    k1[x_][u_,v_]:= mcurvature[x][u,v] +
         Sqrt[Simplify[mcurvature[x][u,v]^2 -
          gcurvature[x][u,v]]]

    k2[x_][u_,v_]:= mcurvature[x][u,v] -
          Sqrt[Simplify[mcurvature[x][u,v]^2 -
           gcurvature[x][u,v]]]

(* Second Fundamental Form and
Gauss Map Programs *)

    eee[x_][u_,v_]:= Simplify[Det[{
        D[x[uu,vv],uu,uu],D[x[uu,vv],uu],D[x[uu,vv],vv]}]/
        Sqrt[Simplify[D[x[uu,vv],uu].D[x[uu,vv],uu]*
             D[x[uu,vv],vv].D[x[uu,vv],vv] -
             D[x[uu,vv],uu].
             D[x[uu,vv],vv]^2]]] /. {uu->u,vv->v}


    fff[x_][u_,v_]:= Simplify[Det[{
        D[x[uu,vv],uu,vv],D[x[uu,vv],uu],D[x[uu,vv],vv]}]/
        Sqrt[Simplify[D[x[uu,vv],uu].D[x[uu,vv],uu]*
             D[x[uu,vv],vv].D[x[uu,vv],vv] -
             D[x[uu,vv],uu].
             D[x[uu,vv],vv]^2]]] /. {uu->u,vv->v}

    ggg[x_][u_,v_]:= Simplify[Det[{
        D[x[uu,vv],vv,vv],D[x[uu,vv],uu],D[x[uu,vv],vv]}]/
        Sqrt[Simplify[D[x[uu,vv],uu].D[x[uu,vv],uu]*
             D[x[uu,vv],vv].D[x[uu,vv],vv] -
             D[x[uu,vv],uu].
             D[x[uu,vv],vv]^2]]] /. {uu->u,vv->v}


    eeee[x_][u_,v_]:= Simplify[Det[{
        D[x[uu,vv],uu,uu],
        D[x[uu,vv],uu],
        D[x[uu,vv],vv]}]] /. {uu->u,vv->v}

    ffff[x_][u_,v_]:= Simplify[Det[{
        D[x[uu,vv],uu,vv],
        D[x[uu,vv],uu],
        D[x[uu,vv],vv]}]] /. {uu->u,vv->v}

    gggg[x_][u_,v_]:= Simplify[Det[{
        D[x[uu,vv],vv,vv],
        D[x[uu,vv],uu],
        D[x[uu,vv],vv]}]] /. {uu->u,vv->v}
 
          
    unitnormal[x_][u_,v_]:= Module[{xu,xv,N1},
        xu = D[x[uu,vv],uu];
        xv = D[x[uu,vv],vv];
        N1 = Cross[xu,xv];
        Simplify[N1.N1]^(-1/2)*N1] /. {uu->u,vv->v}

    jacobian2[x_][u_,v_]:= Module[{xu,xv,N1},
        xu = D[x[uu,vv],uu];
        xv = D[x[uu,vv],vv];
        N1 = Cross[xu,xv];
        Simplify[N1.N1]] /. {uu->u,vv->v}
                  
    weingarten[x_][u_,v_]:=
       {{eee[x][u,v]*gg[x][u,v] - fff[x][u,v]*ff[x][u,v],
         fff[x][u,v]*ee[x][u,v] - eee[x][u,v]*ff[x][u,v]},
        {fff[x][u,v]*gg[x][u,v] - ggg[x][u,v]*ff[x][u,v],
         ggg[x][u,v]*ee[x][u,v] - fff[x][u,v]*ff[x][u,v]}}/
        (ee[x][u,v]*gg[x][u,v] - ff[x][u,v]^2)

(* Vector Analysis Programs *)
 
    der[v_List][ww_][x_,y_]:= {v[[1]]*D[ww[xx,yy],xx],
              v[[2]]*D[ww[xx,yy],yy]} /. {xx->x,yy->y}

    der[v_List][ww_][x_,y_,z_]:= {v[[1]]*D[ww[xx,yy,zz],xx],
                                  v[[2]]*D[ww[xx,yy,zz],yy],
                  v[[3]]*D[ww[xx,yy,zz],zz]} /. {xx->x,yy->y,zz->z}

    grad[g_][x_,y_]:= {D[g[xx,yy],xx],
        D[g[xx,yy],yy]} /. {xx->x,yy->y}

    grad[g_][x_,y_,z_]:=
        {D[g[xx,yy,zz],xx],D[g[xx,yy,zz],yy],
         D[g[xx,yy,zz],zz]} /. {xx->x,yy->y,zz->z}

    gmap[g_][x_,y_]:= grad[g][x,y]/
        (grad[g][x,y].grad[g][x,y])^(1/2)

    gmap[g_][x_,y_,z_]:= grad[g][x,y,z]/
        (grad[g][x,y,z].grad[g][x,y,z])^(1/2)

    divergence[ww_][x_,y_]:= D[ww[xx,yy],xx][[1]] +
            D[ww[xx,yy],yy][[2]] /. {xx->x,yy->y}

    divergence[ww_][x_,y_,z_]:= D[ww[xx,yy,zz],xx][[1]] +
                                D[ww[xx,yy,zz],yy][[2]] +
            D[ww[xx,yy,zz],zz][[3]] /. {xx->x,yy->y,zz->z}

    laplace[g_][x_,y_]:= D[g[xx,yy],{xx,2}] +
        D[g[xx,yy],{yy,2}] /. {xx->x,yy->y}

    laplace[g_][x_,y_,z_]:= D[g[xx,yy,zz],{xx,2}] +
        D[g[xx,yy,zz],{yy,2}] + D[g[xx,yy,zz],{zz,2}] /.
        {xx->x,yy->y,zz->z}

End[]

EndPackage[]
