Surfaces

Exploring differential geometry of surfaces with Mathematica

Professor Mohammad Ghomi
Math 426, PSU, Spring 2004

Introduction

This notebook  contains several miniprograms  for studying the differential geometry of surfaces in  Euclidean space. It also loads up the programs written by Alfred Gray, including the definitions for dozens of surfaces. To study a specific surace, you may either define it directly, or find its parametrization from the list of preprogramed curves.

How to define a parametrized surface?

To use the programs in this notebook, you first need  to define the curve you wish to study. The general format is

s[u_, v_] := { x[u, v]   ,   y[u, v],  z[u, v] }

The Following cell for, instance,

s[u_, v_] := {u, v, (u^3) - 3 u (v^2)}

defines the "monkey saddle". You may also use any of the predefined parametrizations contained in this notebook.

How to find a parametrization?

If you do not know the coordinate functions for the curve which you wish to study, but  know the name of the curve, then the chances are that the formula for it is already included in the Gray's package CURVES which is loaded below. To check that this is the case, just type in the name followed by "?". For instance:

?monkeysaddle

So the monkey saddle is already in the package. To see the formula for it, just type:

monkeysaddle[u, v]

To define our surface, we may then simple set:

s[u_, v_] := {u, v, (u^3) - 3 u (v^2)}

The complete list of preprogramed surfaces is available in a following section. Also, to see the list you may simply type "curves" at any time.

surfs

Examples of  parametrizations

Here are some other examples of standard curves:

s[u_, v_] := torus[3, 1][u, v];
(Subscript[u, i]) = 0; (Subscript[u, f]) = 2  Pi;
(Subscript[v, i]) = 0; (Subscript[v, f]) = 2  Pi;

s[u_, v_] := {u, v, u ^2 + v^2};
(Subscript[u, i]) = -1; (Subscript[u, f]) = 1;
(Subscript[v, i]) = -1; (Subscript[v, f]) = 1;

s[u_, v_] := monkeysaddle[u, v];
(Subscript[u, i]) = -1; (Subscript[u, f]) = 1;
(Subscript[v, i]) = -1; (Subscript[v, f]) = 1;

s[u_, v_] := 
        monkeysaddle[u Cos[v],  u Sin[v]];
(Subscript[u, i]) = 0; (Subscript[u, f]) = 1;
(Subscript[v, i]) = 0; (Subscript[v, f]) = 2  Pi;

s[u_, v_] := {u, v, Sin[u] Sin[v]};
(Subscript[u, i]) := -Pi; (Subscript[u, f]) := Pi;
(Subscript[v, i]) := -Pi; (Subscript[v, f]) := Pi;

s[u_, v_] := ellipsoid[1, 1, 3][u, v];
(Subscript[u, i]) := -Pi; (Subscript[u, f]) := Pi;
(Subscript[v, i]) := -Pi/2; (Subscript[v, f]) := Pi/2;

Initialization

General Initialization

Off[General::spell1]

Off[ParametricPlot::ppcom]

Off[$MaxExtraPrecision::meprec]

Off[Graphics::gptn]

<< Graphics`Colors`

<< Graphics`Shapes`

<< Graphics`Arrow`

Gray's Packages

"PLTPROGS" (Plotting Enhancements for Mathematica)

(*:Sat May 23 21:47:33 MET DST 1998*)

(*:Title: PLTPROGS.m *)

(*:Author: Alfred Gray *)

(* Copyright 1994-1998 by Alfred Gray *)

(*:Package Version: 1.0 *)

(*:Mathematica Versions: 2.2, 3.0 *)

(*:Summary: 
This package consists of several plotting enhancements 
useful for plotting curves and surfaces.
*)

(*:Keywords:
borderedpolygon, centerlines, checkerboard, dip, elrot, 
filledcurve, graphics2dtosurfrev, lathe, linetolines, 
nomesh, orth4d, paramcurveplot1, paramcurveplot2, 
plotpatch2d, polygonholes, primitiveplot3d, selectcolgraphics3d,
 selectgraphics3d, ShowInNotebook, surfacecurveplot, stereo4d, 
 wireframe, WriteAcroLines, WriteAcroPoints, Write2COLMESH, 
 Write3DCMESH, Write3DCNMESH, Write3DMESH
*)

BeginPackage["PLTPROGS`"]

     pltprogs:=Print[Select[Names["PLTPROGS`*"],#!="pltprogs"&]]

(* Miscellaneous Plotting Programs *)

     barycenter::usage="barycenter[Polygon[{pt1,...,ptn}]] subdivides a polygon using its barycenter."

     borderedpolygon::usage="borderedpolygon[intsty,borsty,versty,pts] creates from the list pts a polygonal region whose vertices are the elements of pts. The directives constituting intsty are applied to the interior, the directives constituting borsty are applied to the boundary, and the directives constituting versty are applied to the vertices."

     centerlines::usage="centerlines[...] has the same syntax as ParametricPlot3D[...]. It plots a parametric patch in R^3, but adds lines connecting the center of each quadrilateral to its vertices. centerlines[graphic] adds lines connecting the center of each quadrilateral to its vertices."

     checkerboard::usage="checkerboard[...] has the same syntax as ParametricPlot3D[...]. It plots a parametric patch in R^3 in a checker board fashion.  checkerboard[graphic] gives a checkerboard pattern to a graphic."

    dip::usage="DisplayFunction->dip[ins] inserts a Graphics3D primitive at the beginning of a Graphics3D list when the list is displayed. For example, DisplayFunction->dip[EdgeForm[]] displays polygons with no edges."

    elrot::usage="elrot[p,q][theta][aa_List] is the elementary rotation by an angle theta of pth and qth elements of the list aa."

    filledcurve::usage="filledcurve parametrically plots a plane curve and fills it in. It has the same syntax as ParametricPlot.  filledcurve[graphic] fills in a graphic."

    graphics2dtosurfrev::usage="graphics2dtosurfrev[r,n][list] constructs a surface of revolution of radius r from a Graphics list of lines. Each circle of latitude is approximated by a regular polygon with n vertices."

    lathe::usage="lathe[dots] creates a surface of revolution with no edges from a list of points in the plane."

    linetolines::usage="linetolines[Line[{pt1,...,ptn}]] converts a broken line segment Line[{pt1,...,ptn}] into a list of line segments connecting consecutive pairs of the points pt1,...,ptn."

    nomesh::usage="nomesh[graphic] converts a graphic to a graphic with no edges."

    orth4d::usage="orth4d[surf,{u,umin,umax},{v,vmin,vmax},opts] uses GraphicsArray to plot a surface in R^4 via the four orthogonal projections from R^4 to R^3."

    paramcurveplot1::usage="paramcurveplot1 has the same syntax as ParametricPlot3D, except that a PlotPoints specification is obligatory. It plots the parameter lines u->x[u,v0] of a parametric patch surface:U->R^3."

    paramcurveplot2::usage="paramcurveplot2 has the same syntax as
ParametricPlot3D, except that a PlotPoints specification is obligatory. It plots the parameter lines v->x[u0,v] of a parametric patch surface:U->R^3."

    plotpatch2d::usage="plotpatch2d has the same syntax as ParametricPlot. It plots the images of horizontal and vertical lines under the function f:U->R^2. plotpatch2d has the same functionality as Maeder's CartesianMap and PolarMap, but draws curved instead of straight lines."

    polygonholes::usage="polygonholes[fac][...] has the same syntax as ParametricPlot3D[...]. It plots a parametric patch in R^3 with holes in each polygon, where fac is the ratio of the size of each hole to that of the original polygon containing it.  polygonholes[fac][graphic]  puts a hole in each polygon, where fac is the ratio of the size of each hole to that of the original polygon containing it."

    selectgraphics3d::usage="selectgraphics3d[graphics3dobj,bound,opts] selects from a list of polygons those whose vertices have their coordinates bounded by bound, and then displays the results."

    selectcolgraphics3d::usage="selectcolgraphics3d[graphics3dobj,bound,opts] selects from a list of colored polygons those whose vertices have their coordinates bounded by bound, and then displays the results."

    stereo4d::usage="stereo4d[{a,b,c,d},surf,{u,umin,umax},{v,vmin,vmax},opts] uses GraphicsArray to plot a surface in R^4 via the four stereographic projections from R^4 to R^3 with respect to the point {a,b,c,d} in R^4."

    ShowInNotebook::usage="ShowInNotebook[graphics,nbopts] displays a graphic in a separate notebook."

    surfacecurveplot::usage="surfacecurveplot[surface,{u,u0,u1},{v,v0,v1}, PlotLines->{ulines,vlines},PlotPoints->{upoints,vpoints},opts] plots both the u- and v-parameter curves of a patch x:U->R^3. The number of curves is specified with PlotLines, and the resolution of the lines is specified with PlotPoints."

    wireframe::usage="wireframe[...] has the same syntax as ParametricPlot3D[...]. It plots a wireframe of a parametric patch in R^3.  wirframe[graphic] converts a solid graphic into a wireframe graphic."

(* Mathematica to Acrospin Programs *)

primitiveplot3d::usage="primitiveplot3d[{fx,fy,fz},{t,tmin,tmax}]
produces a list of points of a space curve parametrized by a variable
t which runs from tmin to tmax. primitiveplot3d[{fx,fy,fz},{t,tmin,tmax},{u,umin,umax}] produces a list of points by polygons of a surface parametrized by
t and u."

WriteAcroLines::usage="WriteAcroLines[{color},file,{fx,fy,fz},
{t,tmin,tmax}] outputs to a file in acrospin line format the
parametric plot of a space curve. WriteAcroLines[{color1,color2},file,{fx,fy,fz},{t,tmin,tmax},{u,umin,umax}] outputs to a file in acrospin line format the parametric plot of a surface."

WriteAcroPoints::usage="WriteAcroPoints[{color},file,{fx,fy,fz},
{t,tmin,tmax}] outputs to a file in acrospin point format the
parametric plot of a space curve.
WriteAcroPoints[{color},file,{fx,fy,fz},{t,tmin,tmax},{u,umin,umax}]
outputs to a file in acrospin point format the parametric plot of a surface."

(* Mathematica to Geomview Programs *)

Write3DMESH::usage="Write3DMESH[file,expression,{u,umin,umax},{v,vmin,vmax},
PlotPoints->{udiv,vdiv},opts] outputs to a file in mesh format the
parametric plot of a surface."

Write3DCMESH::usage="Write3DCMESH[file,{expression,color},
{u,umin,umax},{v,vmin,vmax},PlotPoints->{udiv,vdiv},opts] outputs to a file
in mesh format the parametric plot of expression colored by color."

Write3DCNMESH::usage="Write3DCNMESH[file,{expression,unitn,color},
{u,umin,umax},{v,vmin,vmax},
PlotPoints->{udiv,vdiv},opts] outputs to a file in mesh format the
parametric plot of expression colored by color with normal information unitn."

Write2COLMESH::usage="Write2COLMESH[{color1,color2},file,expression,
{u,umin,umax},{v,vmin,vmax},
PlotPoints->{udiv,vdiv},opts] outputs to a file in mesh format the parametric plot of a surface with one side color1 and the other color2."

Begin["`Private`"]

(* Miscellaneous Plotting Programs *)

  barycenter[Polygon[pts_]]:=
        Module[{ptstmp},
        ptstm=Append[pts,First[pts]];
	Table[Polygon[{ptstm[[i]],ptstm[[i + 1]],
        Plus@@pts/Length[ptstm]}],{i,Length[ptstm] - 1}]]

  borderedpolygon[intsty_,borsty_,versty_,pts_]:=
        {Flatten[{intsty,Polygon[pts]}],
         Flatten[{borsty,Line[Append[pts,First[pts]]]}],
         Flatten[{versty,Point/@pts}]}

    centerlines[surface_,{u_,u0_,u1_},{v_,v0_,v1_},opts___]:=
        Module[{plottmp,grtmp},
                plottmp = ParametricPlot3D[
                surface//Evaluate,{u,u0,u1},{v,v0,v1},
                DisplayFunction->Identity,opts];
                grtmp = plottmp /.
     (Polygon[pts_]:>
     {Polygon[Append[pts[[{1,2}]],Plus@@pts/4]],
      Polygon[Append[pts[[{2,3}]],Plus@@pts/4]],
      Polygon[Append[pts[[{3,4}]],Plus@@pts/4]],
      Polygon[Append[pts[[{4,1}]],Plus@@pts/4]]});
        Show[grtmp,DisplayFunction->$DisplayFunction]]

    centerlines[graphic_]:=Show[graphic /. (Polygon[pts_]:>
     {Polygon[Append[pts[[{1,2}]],Plus@@pts/4]],
      Polygon[Append[pts[[{2,3}]],Plus@@pts/4]],
      Polygon[Append[pts[[{3,4}]],Plus@@pts/4]],
      Polygon[Append[pts[[{4,1}]],Plus@@pts/4]]})]

    checkerboard[surface_,{u_,u0_,u1_},{v_,v0_,v1_},opts___]:=
        Module[{plottmp,grtmp},
                plottmp = ParametricPlot3D[
                surface,{u,u0,u1},{v,v0,v1},
                DisplayFunction->Identity,opts];
                grtmp = plottmp /.(Polygon[pts_]:> 
                Polygon[(pts + RotateRight[pts])/2]);
        Show[grtmp,DisplayFunction->$DisplayFunction]]

    checkerboard[graphic_]:=Show[graphic /. (Polygon[pts_]:>
        Polygon[(pts + RotateRight[pts])/2])]

    dip[ins_][g_]:= $DisplayFunction[Insert[g,ins,{1,1}]]

    elrot[p_,q_][theta_][aa_List]:=ReplacePart[ReplacePart[
        aa,aa[[p]]*Cos[theta] + aa[[q]]*Sin[theta],p],
        -aa[[p]]*Sin[theta] + aa[[q]]*Cos[theta],q]

    filledcurve[curve_,{u_,u0_,u1_},opts___]:=
        Module[{plottmp,grtmp},
            plottmp = ParametricPlot[
            curve//Evaluate,{u,u0,u1},AspectRatio->Automatic,
            DisplayFunction->Identity,opts];
            grtmp = plottmp /. (Line[pts_]:> Polygon[pts]);
        Show[grtmp,DisplayFunction->$DisplayFunction]]
    
    filledcurve[graphic_]:=Show[graphic /. (Line[pts_]:> Polygon[pts])]

    graphics2dtosurfrev[r_,n_][Line[{{a_,b_},{c_,d_}}]]:=
        Table[Polygon[{
        {r*a*Cos[2*Pi*(k - 1)/n],r*a*Sin[2*Pi*(k - 1)/n],b},
        {r*a*Cos[2*Pi*k/n],r*a*Sin[2*Pi*k/n],b},
        {r*c*Cos[2*Pi*k/n],r*c*Sin[2*Pi*k/n],d},
        {r*c*Cos[2*Pi*(k - 1)/n],
                r*c*Sin[2*Pi*(k - 1)/n],d}}],{k,n}]

   lathe[dots_]:=Show[Graphics3D[
      Map[graphics2dtosurfrev[0.5,60],
      linetolines[Line[dots]]]],Boxed->False,
      DisplayFunction->dip[EdgeForm[]]]

    linetolines[Line[pts_]]:= Table[Line[{pts[[i]],
             pts[[i + 1]]}],{i,Length[pts] - 1}]

   nomesh[graphics_]:=
      Show[graphics,DisplayFunction->dip[EdgeForm[]]]

    orth4d[surf_,{u_,umin_,umax_},{v_,vmin_,vmax_},opts___]:=
        Module[{g1,g2,g3,g4},
        g1=ParametricPlot3D[Evaluate[
                Drop[surf,{1}]],{u,umin,umax},{v,vmin,vmax},
                DisplayFunction->Identity,opts];
        g2=ParametricPlot3D[Evaluate[
                Drop[surf,{2}]],{u,umin,umax},{v,vmin,vmax},
                DisplayFunction->Identity,opts];
        g3=ParametricPlot3D[Evaluate[
                Drop[surf,{3}]],{u,umin,umax},{v,vmin,vmax},
                DisplayFunction->Identity,opts];
        g4=ParametricPlot3D[Evaluate[
                Drop[surf,{4}]],{u,umin,umax},{v,vmin,vmax},
                DisplayFunction->Identity,opts];
        Show[GraphicsArray[{{g1,g2},{g3,g4}}],
                GraphicsSpacing->{0.2,0},
                DisplayFunction->$DisplayFunction]]

    paramcurveplot1[surface_,{u_,u0_,u1_},{v_,v0_,v1_},
        PlotPoints->{upoints_,vpoints_},opts___]:=
        ParametricPlot3D[Table[
        surface,{v,v0,v1,(v1-v0)/(vpoints - 1)}]//Evaluate,
        {u,u0,u1},PlotPoints->upoints,opts]

    paramcurveplot2[surface_,{u_,u0_,u1_},{v_,v0_,v1_},
        PlotPoints->{upoints_,vpoints_},opts___]:=
        ParametricPlot3D[Table[
        surface,{u,u0,u1,(u1-u0)/(upoints - 1)}]//Evaluate,
        {v,v0,v1},PlotPoints->vpoints,opts]

    plotpatch2d[expr_,{x_,xmin_,xmax_},{y_,ymin_,ymax_},
        PlotPoints->{xpoints_,ypoints_},opts___]:=
        Show[ParametricPlot[Evaluate[
        Table[ComplexExpand[
        {Re[expr],Im[expr]}],
        {y,ymin,ymax,(ymax - ymin)/(ypoints - 1)}]],
        {x,xmin,xmax},PlotPoints->xpoints,
        DisplayFunction->Identity],
        ParametricPlot[Evaluate[
        Table[ComplexExpand[
        {Re[expr],Im[expr]}],
        {x,xmin,xmax,(xmax - xmin)/(xpoints -1)}]],
        {y,ymin,ymax},PlotPoints->ypoints,
        DisplayFunction->Identity],
        AspectRatio->Automatic,
        DisplayFunction->$DisplayFunction,opts]

    polygonholes[fac_][surface_,{u_,u0_,u1_},{v_,v0_,v1_},opts___]:=
        Module[{plottmp,grtmp},
                plottmp = ParametricPlot3D[
                surface,{u,u0,u1},{v,v0,v1},
                DisplayFunction->Identity,opts];
                grtmp = plottmp /.(Polygon[pts_]:> 
                (cntr=(Plus@@pts)/Length[pts];
                sp=(fac*(# - cntr) + cntr)&/@pts;
                Polygon/@Transpose[{pts,RotateLeft[pts],
                RotateLeft[sp],sp}]));
        Show[grtmp,DisplayFunction->$DisplayFunction]]

    polygonholes[fac_][graphic_]:=Show[graphic /. (Polygon[pts_]:>
        (cntr=(Plus@@pts)/Length[pts];
                sp=(fac*(# - cntr) + cntr)&/@pts;
                Polygon/@Transpose[{pts,RotateLeft[pts],
                RotateLeft[sp],sp}]))]

    selectgraphics3d[graphics3dobj_,bound_,opts___]:=
        Show[Graphics3D[Select[graphics3dobj,
        (Abs[#[[1,1,1]]]<bound && Abs[#[[1,1,2]]]<bound &&
         Abs[#[[1,1,3]]]<bound && Abs[#[[1,2,1]]]<bound &&
         Abs[#[[1,2,2]]]<bound && Abs[#[[1,2,3]]]<bound &&
         Abs[#[[1,3,1]]]<bound && Abs[#[[1,3,2]]]<bound &&
         Abs[#[[1,3,3]]]<bound && Abs[#[[1,4,1]]]<bound &&
         Abs[#[[1,4,2]]]<bound && Abs[#[[1,4,2]]]<bound
         )&]],opts]

    selectcolgraphics3d[graphics3dobj_,bound_,opts___]:=
        Show[Graphics3D[Select[graphics3dobj,
        (Abs[#[[2,1,1,1]]]<bound && Abs[#[[2,1,1,2]]]<bound &&
        Abs[#[[2,1,1,3]]]<bound && Abs[#[[2,1,2,1]]]<bound &&
        Abs[#[[2,1,2,2]]]<bound && Abs[#[[2,1,2,3]]]<bound &&
        Abs[#[[2,1,3,1]]]<bound && Abs[#[[2,1,3,2]]]<bound &&
        Abs[#[[2,1,3,3]]]<bound && Abs[#[[2,1,4,1]]]<bound &&
        Abs[#[[2,1,4,2]]]<bound && Abs[#[[2,1,4,2]]]<bound
         )&]],opts]

    stereo4d[{a_,b_,c_,d_},surf_,
        {u_,umin_,umax_},{v_,vmin_,vmax_},opts___]:=
        Module[{g1,g2,g3,g4},
        g1=ParametricPlot3D[Evaluate[
                Drop[surf,{1}]/(a-surf[[1]])],
                {u,umin,umax},{v,vmin,vmax},
                DisplayFunction->Identity,opts];
        g2=ParametricPlot3D[Evaluate[
                Drop[surf,{2}]/(b-surf[[2]])],
                {u,umin,umax},{v,vmin,vmax},
                DisplayFunction->Identity,opts];
        g3=ParametricPlot3D[Evaluate[
                Drop[surf,{3}]/(c-surf[[3]])],
                {u,umin,umax},{v,vmin,vmax},
                DisplayFunction->Identity,opts];
        g4=ParametricPlot3D[Evaluate[
                Drop[surf,{4}]/(d-surf[[4]])],
                {u,umin,umax},{v,vmin,vmax},
                DisplayFunction->Identity,opts];
        Show[GraphicsArray[{{g1,g2},{g3,g4}}],
                GraphicsSpacing->{0.2,0},
                DisplayFunction->$DisplayFunction]]

    ShowInNotebook[gr_,nbopts___]:=
    Module[{nb},
        nb = NotebookPut[Notebook[{
            Cell[GraphicsData["PostScript",
                              DisplayString[gr]],
                 "Graphics",
                 ImageSize->{350,350},
                 ImageMargins -> {{0,0},{0,0}},
                 CellMargins -> {{0,0},{0,0}}]},
            nbopts,WindowSize->Fit,Visible -> False,
            ShowCellBracket -> False,WindowElements->{},
            WindowFrameElements->"CloseBox"]];
        SetOptions[nb,WindowSize->Fit];
        SetOptions[nb,Visible->True];
        nb]

    surfacecurveplot[surface_,{u_,u0_,u1_},{v_,v0_,v1_},
        PlotLines->{ulines_,vlines_},
        PlotPoints->{upoints_,vpoints_},opts___]:=
            Show[paramcurveplot1[surface,{u,u0,u1},{v,v0,v1},
                PlotPoints->{upoints,ulines},
                DisplayFunction->Identity,opts],
                paramcurveplot2[surface,{u,u0,u1},{v,v0,v1},
                PlotPoints->{vlines,vpoints},
                DisplayFunction->Identity,opts],
                DisplayFunction->$DisplayFunction]

    wireframe[surface_,{u_,u0_,u1_},{v_,v0_,v1_},opts___]:=
        Module[{plottmp,grtmp},
                plottmp = ParametricPlot3D[
                surface,{u,u0,u1},{v,v0,v1},
                DisplayFunction->Identity,opts];
                grtmp = plottmp /.(Polygon[pts_]:> 
                Line[Append[pts,First[pts]]]);
        Show[grtmp,DisplayFunction->$DisplayFunction]]

    wireframe[graphic_]:=Show[graphic /. (Polygon[pts_]:> 
                Line[Append[pts,First[pts]]])]

(* Mathematica to Acrospin Programs *)

blue             = 1
green            = 2
cyan             = 3
red              = 4
magenta          = 5
brown            = 6
white            = 7
darkgray         = 8
lightblue        = 9
lightgreen       = 10
lightcyan        = 11
lightred         = 12
lightmagenta     = 13
yellow           = 14
intensifiedwhite = 15

primitiveplot3d[exp_,{t_,tmin_,tmax_},opts___]:=
First[ParametricPlot3D[exp,
{t,tmin,tmax},opts,
DisplayFunction->Identity]] /.
(Line[pts_] :> pts)

primitiveplot3d[exp_,{u_,umin_,umax_},{v_,vmin_,vmax_},opts___]:=
First[ParametricPlot3D[exp,
{u,umin,umax},{v,vmin,vmax},opts,
DisplayFunction->Identity]] /.
(Polygon[pts_] :> pts)

WriteAcroLines[{color_},file_String,curve_,
{t_,tmin_,tmax_},opts___]:=
Module[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7,tmp8},
OpenWrite[file];
WriteString[file,"Set Color"];
WriteString[file,color];
WriteString[file,"\n"];
WriteString[file,"EndPointList name X Y Z\n"];
tmp1=primitiveplot3d[curve//Evaluate,
{t,tmin,tmax},opts];
tmp2=Flatten[tmp1,1];
tmp3=Chop[tmp2,10^-3];
tmp4=Map[PaddedForm[#,{4,4}]&,tmp3,{2}];
tmp5=Table[
{StringForm["VERTEX``",j],tmp4[[j]]},
{j,Length[tmp3]}];
tmp6=Map[Flatten[#,1]&,tmp5,{1}];
Scan[
(WriteString[file,First[#]];
Scan[WriteString[file,"   ",#]&,Rest[#]];
WriteString[file,"\n"])&,
tmp6];
WriteString[file,"\nLineList from to\n"];
tmp7=Table[
StringForm["VERTEX`` VERTEX``\n",j,j+1],
{j,Length[tmp3]-1}];
Scan[WriteString[file,#,"\n"]&,tmp7];
Close[file]]

WriteAcroLines[{color1_:12,color2_:14},
file_String,surface_,
{u_,umin_,umax_},{v_,vmin_,vmax_},opts___]:=
Module[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,
tmp7,tmp8,tmp9,tmp10,tmp11,tmp12},
OpenWrite[file];
WriteString[file,"EndPointList name X Y Z\n\n"];
tmp1=primitiveplot3d[surface//Evaluate,
{u,umin,umax},{v,vmin,vmax},opts];
tmp2=Chop[tmp1,10^-3];
tmp3=Map[PaddedForm[#,{4,4}]&,tmp2,{3}];
tmp4=Table[
{StringForm["POLY``VERTEX``",j,k],tmp3[[j,k]]},
{j,Length[tmp3]},{k,4}];
tmp5=Map[Flatten[#,1]&,tmp4,{2}];
tmp6=Flatten[tmp5,1];
Scan[WriteString[file,"   ",#]&,Rest[#]];
Scan[
(WriteString[file,First[#]];
Scan[WriteString[file,"   ",#]&,Rest[#]];
WriteString[file,"\n"])&,
tmp6];
WriteString[file,"\nLineList from to Color Layer\n\n"];
tmp7=Table[
StringForm["POLY``VERTEX`` POLY``VERTEX`` `` ``\n",
j,1,j,2,color1,1],{j,Length[tmp3]}];
tmp8=Table[
StringForm["POLY``VERTEX`` POLY``VERTEX`` `` ``\n",
j,4,j,3,color1,1],{j,Length[tmp3]}];
tmp9=Join[tmp7,tmp8];
Scan[WriteString[file,#,"\n"]&,tmp9];
tmp10=Table[
StringForm["POLY``VERTEX`` POLY``VERTEX`` `` ``\n",
j,1,j,4,color2,2],{j,Length[tmp3]}];
tmp11=Table[
StringForm["POLY``VERTEX`` POLY``VERTEX`` `` ``\n",
j,3,j,2,color2,2],{j,Length[tmp3]}];
tmp12=Join[tmp10,tmp11];
Scan[WriteString[file,#,"\n"]&,tmp12];
Close[file]]

WriteAcroPoints[{color_},file_String,curve_,
{t_,tmin_,tmax_},opts___]:=
Module[{tmp1,tmp2,tmp3,tmp4},
OpenWrite[file];
WriteString[file,"Set Color"];
WriteString[file,color];
WriteString[file,"\n"];
WriteString[file,"PointList X Y Z\n"];
tmp1=primitiveplot3d[curve//Evaluate,
{t,tmin,tmax},opts];
tmp2=Union[Flatten[tmp1,1]];
tmp3=Chop[tmp2,10^-3];
tmp4=Map[PaddedForm[#,{4,4}]&,tmp3,{2}];
Scan[
(WriteString[file,First[#]];
Scan[WriteString[file,"   ",#]&,Rest[#]];
WriteString[file,"\n"])&,
tmp4];
Close[file]]

WriteAcroPoints[{color_},file_String,surface_,
{u_,umin_,umax_},{v_,vmin_,vmax_},opts___]:=
Module[{tmp1,tmp2,tmp3,tmp4},
OpenWrite[file];
WriteString[file,"Set Color"];
WriteString[file,color];
WriteString[file,"\n"];
WriteString[file,"PointList X Y Z\n"];
tmp1=primitiveplot3d[surface//Evaluate,
{u,umin,umax},{v,vmin,vmax},opts];
tmp2=Union[Flatten[tmp1,1]];
tmp3=Chop[tmp2,10^-3];
tmp4=Map[PaddedForm[#,{4,4}]&,tmp3,{2}];
Scan[
(WriteString[file,First[#]];
Scan[WriteString[file,"   ",#]&,Rest[#]];
WriteString[file,"\n"])&,
tmp4];
Close[file]]

(* Mathematica to Geomview Programs *)

Write3DMESH[file_String,expression_,
{u_,umin_,umax_},{v_,vmin_,vmax_},
PlotPoints->{udiv_,vdiv_},otmp___]:=
Module[{tmp1,tmp2,tmp3,tmp4},
OpenWrite[file];
WriteString[file,"MESH\n"];
WriteString[file,ToString[vdiv+1]];
WriteString[file," "];
WriteString[file,ToString[udiv+1]];
WriteString[file,"\n"];
tmp1=N[Table[expression,
    {u,umin,umax,(umax-umin)/udiv},
    {v,vmin,vmax,(vmax-vmin)/vdiv}]];
tmp2=Flatten[tmp1,1];
tmp3=Chop[tmp2,10^-4];
tmp4=Map[PaddedForm[#,{4,4}]&,tmp3,{2}];
Scan[
   (WriteString[file,First[#]];
    Scan[WriteString[file,"   ",#]&,Rest[#]];
    WriteString[file,"\n"])&,tmp4];
Close[file]]

Write3DCMESH[file_String,{expression_,color_},
{u_,umin_,umax_},{v_,vmin_,vmax_},
PlotPoints->{udiv_,vdiv_},otmp___]:=
Module[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7},
OpenWrite[file];
WriteString[file,"CMESH\n"];
WriteString[file,ToString[vdiv+1]];
WriteString[file," "];
WriteString[file,ToString[udiv+1]];
WriteString[file,"\n"];
tmp1=N[Table[{expression,HSBColor[color,1,1]},
    {u,umin,umax,(umax-umin)/udiv},
    {v,vmin,vmax,(vmax-vmin)/vdiv}]];
tmp2=Apply[List,tmp1,{3}];
tmp3=Map[Flatten[#]&,tmp2,{2}];
tmp4=Flatten[tmp3,1];
tmp5=Map[Append[#,1]&,tmp4];
tmp6=Chop[tmp5,10^-4];
tmp7=Map[PaddedForm[#,{4,4}]&,tmp6,{2}];
Scan[
   (WriteString[file,First[#]];
    Scan[WriteString[file,"   ",#]&,Rest[#]];
    WriteString[file,"\n"])&,tmp7];
Close[file]]



Write3DCNMESH[file_String,{expression_,unitn_,color_},
{u_,umin_,umax_},{v_,vmin_,vmax_},
PlotPoints->{udiv_,vdiv_},otmp___]:=
Module[{tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,tmp7},
OpenWrite[file];
WriteString[file,"CNMESH\n"];
WriteString[file,ToString[vdiv+1]];
WriteString[file," "];
WriteString[file,ToString[udiv+1]];
WriteString[file,"\n"];
tmp1=N[Table[{expression,unitn,HSBColor[color,1,1]},
    {u,umin,umax,(umax-umin)/udiv},
    {v,vmin,vmax,(vmax-vmin)/vdiv}]];
tmp2=Apply[List,tmp1,{3}];
tmp3=Map[Flatten[#]&,tmp2,{2}];
tmp4=Flatten[tmp3,1];
tmp5=Map[Append[#,1]&,tmp4];
tmp6=Chop[tmp5,10^-4];
tmp7=Map[PaddedForm[#,{4,4}]&,tmp6,{2}];
Scan[
   (WriteString[file,First[#]];
    Scan[WriteString[file,"   ",#]&,Rest[#]];
    WriteString[file,"\n"])&,tmp7];
Close[file]]

Write2COLMESH[file_String,expression_,
{u_,umin_,umax_},{v_,vmin_,vmax_},
PlotPoints->{udiv_,vdiv_},otmp___]:=
Module[{tmp1,tmp2,tmp3,tmp4},
OpenWrite[file];
WriteString[file,"appearance {\n"];
WriteString[file,"material {\n"];
WriteString[file,"emission 1 1 0\n"];
WriteString[file,"diffuse  1 0 0}}\n"];
WriteString[file,"MESH\n"];
WriteString[file,ToString[vdiv+1]];
WriteString[file," "];
WriteString[file,ToString[udiv+1]];
WriteString[file,"\n"];
tmp1=N[Table[expression,
    {u,umin,umax,(umax-umin)/udiv},
    {v,vmin,vmax,(vmax-vmin)/vdiv}]];
tmp2=Flatten[tmp1,1];
tmp3=Chop[tmp2,10^-4];
tmp4=Map[PaddedForm[#,{4,4}]&,tmp3,{2}];
Scan[
   (WriteString[file,First[#]];
    Scan[WriteString[file,"   ",#]&,Rest[#]];
    WriteString[file,"\n"])&,tmp4];
Close[file]]

Write2COLMESH[{color1_,color2_},file_String,expression_,
{u_,umin_,umax_},{v_,vmin_,vmax_},
PlotPoints->{udiv_,vdiv_},otmp___]:=
Module[{cc1,cc2,cc3,dd1,dd2,dd3,df,em,tmp1,tmp2,tmp3,tmp4},
OpenWrite[file];
WriteString[file,"appearance {\n"];
WriteString[file,"material {\n"];
{cc1,cc2,cc3}=List@@color1;
em=StringForm["emission `` `` ``\n",cc1,cc2,cc3];
WriteString[file,em];
{dd1,dd2,dd3}=List@@color2;
df=StringForm["diffuse `` `` ``}}\n",dd1,dd2,dd3];
WriteString[file,df];
WriteString[file,"MESH\n"];
WriteString[file,ToString[vdiv+1]];
WriteString[file," "];
WriteString[file,ToString[udiv+1]];
WriteString[file,"\n"];
tmp1=N[Table[expression,
    {u,umin,umax,(umax-umin)/udiv},
    {v,vmin,vmax,(vmax-vmin)/vdiv}]];
tmp2=Flatten[tmp1,1];
tmp3=Chop[tmp2,10^-4];
tmp4=Map[PaddedForm[#,{4,4}]&,tmp3,{2}];
Scan[
   (WriteString[file,First[#]];
    Scan[WriteString[file,"   ",#]&,Rest[#]];
    WriteString[file,"\n"])&,tmp4];
Close[file]]

End[]

EndPackage[]

"CSPROGS" (Programs for Manipulating Curves and Surfaces)

(*:Title: CSPROGS *)

(*:Author: Alfred Gray *)

(* Copyright 1994  by Alfred Gray *)

(*:Package Version: 1.0 *)

(*:Mathematica Version: 2.2 *)

(*:Summary: 
	This package consists of miniprgrams to compute curvature, torsion and related geometric quantities associated with curves and surfaces.
*)

(*:Keywords:
analytictominimal, arclength, arclengthprime, asymeq1, asymeq2, binormal, binormalsurf, brioschicurvature, christoffel, circle3point, cone, Cross, ctor, curveonsurf, cylinder, darbouxvector, der, divergence, ee, eee, eeee, evolute, ff, fff, ffff, focalkmu, focalkpi, focalk1, focalk2, frenetframe, gaussiancurvature, gcurvature, genhel, gensurfrev, geo, geodesiccurvature, geodesictorsion, geomincurve, gg, ggg, gggg, gmap, grad, harmonictoanalytic, helical, infarea, inscurve, intrinsic, inversion, inversioncurve, inversionsurf, involute, involutends, J, jacobian2, kappa, kappan, kappa2, kappa2c, kappa2implicit, kappa2polar, kappa23d, kmu, kpi, k1, k2, laplace, length, lengthn, lengthnds, mcurvature, meancurvature, meancurvaturevf, monge, mongepolar, nateq, normal, normalcurvature, normalline, normalsurf, osculatingcircle, paramtoimp, parcurve, parsurf, pedal, pedsurf, perpsurf, plotintrinsic, plotintrinsic3d, plottl, plottp, polarize, rf, rollcurve, rt, rtoc, ruled, seashell, solvegeoeqs, stereoinj, stereoinjcomplex, stereoproj, stretch, surfacejacobian, surfacejacobianmatrix, surfacenormal, surfast, surfrev, tandev, tangent, tangentkmu, tangentkpi, tangentline, tangentplane, tau, tl, totalscalarcurvature, transsurf, tubecurve, turningangle, turningnumber, twist, umeqs, unitnormal, vartube, wei, weicurvature, weimetric, weingarten, weiprime, writhe
*)

Print["For a list of programs type csprogs."]

BeginPackage["CSPROGS`"]

    csprogs:=Print[Select[Names["CSPROGS`*"],#!="csprogs"&]]

(* Arc Length Functions *)

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

    arclength::usage="t->arclength[alpha][t] is the arc length 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."

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

    lengthnds::usage="t->lengthnds[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]."

(* 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."

(* 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[v,w] is the cross product of the vectors
v={v1,v2,v3} and w={w1,w2,w3}."

(* Construction of Plane Curves from
Plane Curves  *)

    evolute::usage="t->evolute[alpha][t] is the
evolute of the plane curve alpha.  t->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]."

     involutends::usage="t->involutends[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[alpha,t0,{tmin,tmax},opts] plots
the tangent line emanating from alpha at alpha[t0]."

    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[{a,b},alpha][t] is the pedal curve of the
plane curve alpha with respect to the point {a,b}.  t->pedal[alpha][t]
is the pedal curve of the plane curve alpha with respect to the origin."

    rollcurve::usage="t->rollcurve[alpha,b,h,c][t] is the parametrized
curve traced out by a point P attached to a circle of radius b
rolling along a plane curve alpha.  The distance of P to the center of the
rolling circle is h, and alpha[c] is the point at which the circle starts
rolling."

(* Construction of Space Curves from
Plane Curves  *)

    curveonsurf::usage="t->curveonsurf[surf,alpha][t] is the composition
of the patch surf:U->R^3 with the curve alpha:(a,b)->U."

    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."

    inscurve::usage="t->inscurve[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[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  *)

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

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

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

    transsurf::usage="{u,v}->transsurf[alpha,gamma][u,v] is the
translation surface {u,v}->alpha[u,v] + gamma[u,v]
generated by the curve alpha,gamma:(a,b)->R^n."

    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."

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

(* Curvature and Torsion of Curves *)

    geodesiccurvature::usage="geodesiccurvature[surf,alpha][t] 
calculates the geodesic curvature of the curve t->surf[alpha[t]] 
contained in the parametrized surface surf:U->R^3, where alpha
is a plane curve whose trace is contained in the subset U of R^2."

    geodesictorsion::usage="geodesictorsion[surf,alpha][t] 
calculates the geodesic torsion of the curve t->surf[alpha[t]] 
contained in the parametrized surface surf:U->R^3, where alpha
is a plane curve whose trace is contained in the subset U of R^2."

    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."

    normalcurvature::usage="normalcurvature[surf,alpha][t] 
calculates the normal curvature of the curve t->surf[alpha[t]] 
contained in the parametrized surface surf:U->R^3, where alpha
is a plane curve whose image is contained in the subset U of R^2."

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

    totalscalarcurvature::usage="totalscalarcurvature[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 *)

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

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

    frenetframe::usage="t->frenetframe[alpha][t] is the
Frenet frame 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,{u,v},{u0:0,v0:0},tfin,div,optsnd___]
solves the geodesic equations of the patch 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 in R^n of radius rho centered at p.
inversion[rho][x] is the inversion of a point x with respect to a
sphere in R^n 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 in R^n of radius rho
centered at p.  t->inversioncurve[rho][alpha][t] is the inversion of a
curve alpha with respect to a sphere in R^n 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 in R^n of
radius rho centered at p.  {u,v}->inversionsurf[rho][surf][u,v] is the
inversion of a patch surf with respect to a sphere in R^n of radius rho
centered at the origin."

(* Jacobians and Normals *)

    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."

    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 product of the Jacobian matrix of the
patch x:U->R^n with its transpose."

    surfacenormal::usage="{u,v}->surfacenormal[x][u,v] is the
cross product of the tangent vectors xu and xv of the patch x:U->R^3."

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

(* Metric (First Fundamental 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]."

     geomincurve::usage="z->geomincurve[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."

    weimetric::usage="z->weimetric[f,g][z] is the scaling function of the
metric of the Weierstrass patch determined by meromorphic functions f and g."

    weiprime::usage="z->weiprime[f,g][z] is integrand of the Weierstrass
representation."

(* 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.  The Gauss symbol is D."

    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.  The Gauss symbol is D'."

    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.  The Gauss symbol is D''."

    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."

    gmap::usage="gmap[g][x,y] is the Gauss map for a level curve of
a function g:R^2->R.  gmap[g][x,y,z] is the Gauss map for a level surface of
a function g:R^3->R."

    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."

(* 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."

    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}."

    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,
where U is an open subset of R^2."

    mongepolar::usage="{r,theta}->mongepolar[g][r,theta] is the parametrization
{r,theta}->{r*Cos[theta],r*Sin[theta],g[r,theta]} of the graph of g:U->R,
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="{r,theta}->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 {u,v}->func[u,v] is an n-tuple of functions."

    stereoinj::usage="p->stereoinj[a][p] is the stereographic
injection from R^n to a sphere of radius a centered at the origin
in R^(n+1)."

    stereoinjcomplex::usage="z->stereoinjcomplex[g][z] is the composition
stereographic of stereoinj:C->S^2(1) with a function g:C->C."

    stereoproj::usage="q->stereoproj[a][q] is the stereographic
projection from the sphere of radius a centered at the origin in
R^(n+1) to R^n."

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

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}]

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

    lengthnds[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]]

(* 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

(* Complex Structure *)

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

(* Cross Product *)

    Cross[v_List,w_List]:=
        RotateLeft[v*RotateLeft[w] -
        RotateLeft[v]*w] /;
        Length[v] == Length[w] == 3

(* Construction of Plane Curves from
Plane Curves  *)
                
    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]

    involutends[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] + D[alpha[tt],tt].D[alpha[tt],tt]/
        (D[alpha[tt],{tt,2}].
        J[D[alpha[tt],tt]])*J[D[alpha[tt],tt]],
        Abs[(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]

    rollcurve[alpha_,b_,h_,c_][t_]:= (alpha[tt] - 
        h*Sin[(arclength[alpha][tt]-arclength[alpha][c])/b]*
        D[alpha[tt],tt]/Sqrt[Factor[D[alpha[tt],tt].
        D[alpha[tt],tt]]] +
        (b - h*Cos[(arclength[alpha][tt]-arclength[alpha][c])/b])*
        J[D[alpha[tt],tt]]/Sqrt[Factor[D[alpha[tt],tt].
        D[alpha[tt],tt]]]) /. tt->t

(* Construction of Space Curves from
Plane Curves  *)

     curveonsurf[surf_,alpha_][t_]:=
             surf[alpha[t][[1]],alpha[t][[2]]]

     helical[alpha_,c_][t_]:=
        {alpha[t][[1]],alpha[t][[2]],c t}

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

     inscurve[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,v*q2,v*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[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  *)

    binormalsurf[alpha_][u_,v_]:=
        alpha[u] + v*binormal[alpha][u]
         
    normalsurf[alpha_][u_,v_]:=
        alpha[u] + v*normal[alpha][u]

    perpsurf[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

    transsurf[alpha_,gamma_][u_,v_]:= alpha[u] + gamma[v]
               
    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}

    meancurvaturevf[x_][u_,v_]:= 
        Module[{h,h1,h2,h3,h4},
        h=gg[x][uu,vv]*D[x[uu,vv],uu,uu] -
           ff[x][uu,vv]*D[x[uu,vv],uu,vv] +
           ee[x][uu,vv]*D[x[uu,vv],vv,vv];
        h1 = (h.D[x[uu,vv],uu])*D[x[uu,vv],uu];
        h2 = (h.D[x[uu,vv],vv])*D[x[uu,vv],vv];
        h3 = h - h1 - h2;
        h4 = 2*(ee[x][uu,vv]*gg[x][uu,vv] -
                ff[x][uu,vv]^2);
         Simplify[h3/h4] /. {uu->u,vv->v}]

(* Curvature and Torsion of Curves *)

    geodesiccurvature[surf_,alpha_][t_]:=
        Module[{unc,tmp,crv1,crv2},
        unc=unitnormal[surf][alpha[tt][[1]],alpha[tt][[2]]];
        tmp=surf[alpha[tt][[1]],alpha[tt][[2]]];
        crv1=D[tmp,tt];crv2=D[tmp,{tt,2}];
        Simplify[PowerExpand[Det[
        {crv2,unc,crv1}]/(crv1.crv1)^(3/2) /. tt->t]]]

    geodesictorsion[surf_,alpha_][t_]:=
        Module[{unc,tmp,crv1,unc1},
        unc=unitnormal[surf][alpha[tt][[1]],alpha[tt][[2]]];
        tmp=surf[alpha[tt][[1]],alpha[tt][[2]]];
        crv1=D[tmp,tt];unc1=D[unc,tt];
        Simplify[PowerExpand[Det[
        {unc,unc1,crv1}]/(crv1.crv1)^(3/2) /. tt->t]]]

    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

    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
    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}

    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

    normalcurvature[surf_,alpha_][t_]:=
        Module[{unc,tmp,crv1,crv2},
        unc=unitnormal[surf][alpha[tt][[1]],alpha[tt][[2]]];
        tmp=surf[alpha[tt][[1]],alpha[tt][[2]]];
        crv1=D[tmp,tt];crv2=D[tmp,{tt,2}];
        Simplify[PowerExpand[(crv2.unc)/(crv1.crv1) /. 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
              
    totalscalarcurvature[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[Evaluate[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[alpha_][t_]:=
         Sqrt[D[alpha[tt],tt].D[alpha[tt],tt]]*(
                tau[alpha][tt]*tangent[alpha][tt] +
                kappa[alpha][tt]*binormal[alpha][tt]) /. tt->t

    frenetframe[alpha_][t_]:=
        {Line[{alpha[t],alpha[t] + tangent[alpha][t]}],
         Line[{alpha[t],alpha[t] + normal[alpha][t]}],
         Line[{alpha[t],alpha[t] + binormal[alpha][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,{u_,v_},{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]*n3[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 and Normals*)

   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}

    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]]]

   surfacenormal[x_][u_,v_]:= Module[{xu,xv},
        xu = D[x[uu,vv],uu];
        xv = D[x[uu,vv],vv];
        Simplify[Cross[xu,xv]] /. {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}

(* Metric (First Fundamental 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]

     geomincurve[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_]:= -16*Abs[D[g[zz],zz]]^2/
            (Abs[f[zz]]^2*(1 + Abs[g[zz]]^2)^4) /. zz->z

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

    weiprime[f_,g_][z_]:= {(1 - g[z]^2)*f[z]/2,
              I*(1 + g[z]^2)*f[z]/2,f[z]*g[z]}

(* 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}

        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}

(* 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]}]

    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))]]

     monge[h_][u_,v_]:= {u,v,h[u,v]}

     mongepolar[g_][r_,theta_]:=
        {r*Cos[theta],r*Sin[theta],g[r,theta]}
     
     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]]
     
     stereoinj[a_][p_]:=Append[2*a^2*p,
                          a*(p.p - a^2)]/(p.p + a^2)

     stereoproj[a_][q_]:=a*Drop[q,-1]/(a - Last[q])

     stereoinjcomplex[g_][z_]:=
            {(2*Re[g[z]])/(Im[g[z]]^2 + Re[g[z]]^2 + 1), 
             (2*Im[g[z]])/(Im[g[z]]^2 + Re[g[z]]^2 + 1), 
             (Im[g[z]]^2 + Re[g[z]]^2 - 1)/
             (Im[g[z]]^2 + Re[g[z]]^2 + 1)}

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

End[]

EndPackage[]

"SURFS" (List of parametrized surfaces)

(*Fri May 29 18:42:40 BST 1998*)

(*:Title: SURFS.m *)

(*:Author: Alfred Gray *)

(* Copyright 1994-1998 by Alfred Gray *)

(*:Package Version: 1.0 *)

(*:Mathematica Versions: 2.2, 3.0 *)

(*:Summary: 
This package consists of definitions of parametrizations of 
surfaces.
*)

(*:Keywords:
agnesirev, apple, astell, astroidmincurve, bohdom, bour, bourmincurve, bourpolar, boy, cardioidmincurve, catalan, catalandef, catalanmincurve, catenoid, catenoidasym, catenoidiso, catenoidmincurve, cayleysexticmincurve, circleinvolutemin, circleinvolutemincurve, circularcone, circularcylinder, circularmembrane, circularmembranemovie, circularmembraneplot, cissoidmincurve, clothoidmin, clothoidmincurve, cnccosurfrev, cnccotcheb, cnchyprofilemincurve, cnchysurfrev, cnchytcheb, cossurface, cossurfaceimplicit, costa, costamincurve, cpcprofilemincurve, cpcsurfrev, crosscap, crosscapimplicit, deltoidinvolutemin, deltoidinvolutemincurve, deltoidmin, deltoidmincurve, dini, dinitcheb, earth, eighteight, eightsurface, ellhypcyclide, ellipsemin, ellipsemincurve, ellipsoid, ellipsoidimplicit, ellipticalhyperboloidminus, ellipticalhyperboloidplus, ellipticcoor, ellipticparaboloid, ellipticparaboloidimplicit, enneper, ennepercatenoidmincurve, ennepercatenoidpolar, ennepermincurve, enneperpolar, epicycloidmincurve, equihom1implicit, equihom2implicit, equihom3implicit, equihom4implicit, equihom5implicit, equihom6implicit, exptwist, exptwistasym, flattorus, funnel, funnelasym, gencube, genhypparab, genoctahedron, goursatimplicit, handkerchief, heart, helicoid, helicoidprin, heltocat, henneberg, hennebergmincurve, hyperbolamincurve, hyperbolicparaboloid, hyperbolicparaboloidimplicit, hyperboloid, hyperboloidimplicit, hypocycloidinvolutemin, hypocycloidinvolutemincurve, hypocycloidmincurve, hy2sheet, hy2sheetimplicit, invertedcircularcone, invertedcircularcylinder, invertedtorus, jorgemeeksmincurve, jorgemeeksmincurveprime, jorgemeekspolar, kidney, kleinbottle, kleinbottlebis, kleinbottlewri, klein4, kuen, kuentcheb, kazoolaimplicit, kummerimplicit, lemniscatemin, lemniscatemincurve, logmincurve, logminpolar, logspiralmincurve, menn, mercatorellipsoid, mercatorsphere, mincurvs, moebcirc, moebiusstrip, monkeysaddle, monkeysaddlepolar, nephroidmincurve, nielsenspiralmincurve, nodoidnds, parabolamin, parabolamincurve, paraboliccoor, paraboliccyclide, paraboloid, pear, pearbis, perturbedms, perturbedmspolar, pillow, plane, planepolar, planetocat, plucker, pluckerpolar, pseudocrosscap, pseudosphere, pseudospherebis, pseudospheretcheb, rconoid, richmond, richmondmincurve, richmondpolar, roman, rosemin, rosemincurve, scherk, scherkimplicit, scherktowermincurveprime, scherk5, scherk5implicit, scmincurve, shoe, shoeasym, sievert, sinovaloid, sinsurface, sinsurfaceimplicit, snail, sphere, sphereimplicit, squaremembrane, squaremembranemovie, squaremembraneplot, stereographicellipsoid, stereographicsphere, stereographicspherepolar, swallowtail, tetrahedral, thomsen, torus, torusasym, torusbis, torusimplicit, tractrixmincurve, twiflat, twisphere, twocuspsimplicit, veronese, wallis, whitneyumbrella, wrinkles
*)

BeginPackage["SURFS`","CURVES`","CSPROGS`","NumericalMath`BesselZeros`"]

    surfs:=Print[Select[Names["SURFS`*"],#!="surfs"&]]

(* Parametrically Defined Surfaces *)

    agnesirev::usage="{u,v}->agnesirev[a][u,v] is a parametrization of the surface of revolution formed from a rotated witch of Agnesi.  To plot try agnesirev[1][0,2Pi,-1,1]."

    apple::usage="{u,v}->apple[u,v] is a parametrization of the surface of revolution that resembles an apple.  To plot try apple[0,3Pi/2,-Pi,Pi]."

    astell::usage="{u,v}->astell[a,b,c][u,v] is a parametrization of an astroidial ellipsoid with axes of lengths a, b and c.  To plot try astell[1,1,1][0,2Pi,-Pi/2,Pi/2,PlotPoints->{40,40},PlotRange->All]."

    bohdom::usage="{u,v}->bohdom[a,b,c][u,v] is a parametrization of a Bohemian dome. It is formed by moving an ellipse along a circle in a perpendicular plane so that the ellipse remains parallel to a plane.  To plot try bohdom[4,4,4][0,2Pi,-Pi,Pi]."

    bour::usage="{u,v}->bour[u,v] is a parametrization of Bour's minimal surface.  To plot try bour[-1,1,-1,1]."

    bourpolar::usage="{r,theta}->bourpolar[n][t][r,theta] is the polar parametrization of a 1-parameter family of minimal surfaces such that {r,theta}->bourpolar[n][0][r,theta] is Bour's minimal surface of degree n. {r,theta}->bourpolar[m,n][t][r,theta] is the polar parametrization of a 1-parameter family of minimal surfaces such that {r,theta}->bourpolar[m,n][0][r,theta] is Bour's minimal surface of degree {m,n}.  To plot try bourpolar[2,1][0][0,1.7,-Pi,Pi,PlotPoints->{20,40},PlotRange->All]."

    boy::usage="{u,v}->boy[a,b,c][u,v] is a parametrization of Boy's surface.  To plot try boy[4,4,1][-Pi/2,Pi/2,-Pi/2,Pi/2]."

    catalan::usage="{u,v}->catalan[a][u,v] is a parametrization of Catalan's minimal surface.  To plot try catalan[1][-2Pi,2Pi,-Pi/2,Pi/2]."

    catalandef::usage="{u,v}->catalandef[a][t][u,v] is a 1-parameter family of minimal surfaces connecting Catalan's minimal surface to its conjugate.  To plot try catalandef[1][Pi/4][-2Pi,2Pi,-Pi/2,Pi/2]."

    catenoid::usage="{u,v}->catenoid[a][u,v] is the minimal surface of revolution generated by t->catenary[a][t].  To plot try catenoid[1][0,3Pi/2,-1.5,1.5]."

    catenoidasym::usage="{p,q}->catenoidasym[a][p,q] is a parametrization of a catenoid by asymptotic curves.  To plot try catenoidasym[1][-Pi,Pi,-Pi,Pi,PlotRange->All]."

    catenoidiso::usage="{u,v}->catenoidiso[a][u,v] is an isothermal parametrization of a catenoid.  To plot try catenoidiso[1][-Pi,Pi,-Pi,Pi,PlotRange->All]."

    circleinvolutemin::usage="{u,v}->circleinvolutemin[1,a][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->circleinvolutemin[1,a][0][u,v] is a minimal surface containing the first involute of a circle as a geodesic.  To plot try circleinvolutemin[1,1][1][0,Pi,-1,1]."

    circularcone::usage="{u,v}->circularcone[a,b][u,v] is a parametrization of a circular cone of radius a and slope b/a.  To plot try circularcone[1,0.5][-Pi,Pi,-1,1]."

    circularcylinder::usage="{u,v}->circularcylinder[a][u,v] is a parametrization of a circular cylinder of radius a.  To plot try circularcylinder[1][-Pi,Pi,-4,4]."

    clothoidmin::usage="{u,v}->clothoidmin[a][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->clothoidmin[a][0][u,v] is a minimal surface containing a clothoid as a geodesic.  To plot try clothoidmin[1][0][-5,5,-0.3,0.3,PlotPoints->{60,10}]."

    cnccosurfrev::usage="{u,v}->cnccosurfrev[a,b][u,v] is a parametrization of a surface of revolution of constant negative curvature -1/a^2 of conic type.  To plot try cnccosurfrev[1,0.83][-Pi,Pi,-0.6,0.6]."

    cnccotcheb::usage="{u,v}->cnccotcheb[a,b][u,v] is a Tchebyshef parametrization of a constant negative curvature surface of revolution of conic type.  To plot try cnccotcheb[1,1.2][-3.5,3.5,-Pi/1.2,Pi/1.2]."

    cnchysurfrev::usage="{u,v}->cnchysurfrev[a,b][u,v] is a parametrization of a surface of revolution of constant negative curvature -1/a^2 of hyperboloid type.  To plot try cnchysurfrev[1,0.83][-Pi,Pi,-0.6,0.6]."

    cnchytcheb::usage="{u,v}->cnchytcheb[a,b][u,v] is a Tchebyshef parametrization of a constant negative curvature surface of revolution of hyperboloid type.  To plot try cnchytcheb[1,0.8][-3.5,3.5,-Pi/0.8,Pi/0.8]."

    cossurface::usage="{u,v}->cossurface[a][u,v] is the surface {u,v}->a*{Cos[u],Cos[v],Cos[u + v]}. The implicit equation is (z - a*x*y) == a^2*(1 - x^2)*(1 - y^2).  To plot try cossurface[1][-Pi,Pi,-Pi,Pi]."

    costa::usage="{u,v}->costa[u,v] is a parametrization of Costa's minimal surface."

    cpcsurfrev::usage="{u,v}->cpcsurfrev[a,b][u,v] is a parametrization of a surface of revolution of constant positive curvature 1/a^2. There are 3 cases: football type (a>b), sphere (a=b), barrel type (a<b).  To plot try cpcsurfrev[1,0.5][-Pi,Pi,-1,1]."

    crosscap::usage="{u,v}->crosscap[a][u,v] is a parametrization of a cross cap in R^3.  To plot try crosscap[1][-Pi/2,Pi/2,-Pi/2,Pi/2]."

    deltoidinvolutemin::usage="{u,v}->deltoidinvolutemin[a][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->deltoidinvolutemin[a][0][u,v] is a minimal surface containing the involute of a deltoid as a geodesic.  To plot try deltoidinvolutemin[1][0][-2Pi,2Pi,-0.3,0.3,
PlotRange->All,PlotPoints->{40,10}."

    deltoidmin::usage="{u,v}->deltoidmin[a][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->deltoidmin[a][0][u,v] is a minimal surface containing a deltoid as a geodesic.  To plot try deltoidmin[1][0][-Pi,Pi,-0.3,0.3,PlotRange->All,PlotPoints->{40,10}]."

    dini::usage="{u,v}->dini[a,b][u,v] is the standard parametrization of Dini's surface of constant negative curvature -1/a^2. It is the generalized helicoid of slant b generated by a tractrix. The case b=0 is the standard parametrization of a pseudosphere.  To plot try dini[1,0.2][0,4Pi,0.05,1,
PlotPoints->{55,10},DisplayFunction->dip[EdgeForm[]]]."

    dinihelicoid::usage="{u,v}->dinihelicoid[pm,a,b,c][u,v] is a generalized helicoid of constant negative curvature -1/a^2 and slope c. The parameter pm is plus or minus 1."

    dinitcheb::usage="{u,v}->dinitcheb[a,b][u,v] is a Tchebyshef principal patch and {p,q}->dinitcheb[a,b][p + q,p - q] is a Tchebyshef asymptotic patch for Dini's surface of constant negative curvature -1/a^2.  To plot try dinitcheb[1,0.2][-3,3,0,6Pi,PlotPoints->{10,55},DisplayFunction->dip[EdgeForm[]]]."

    earth::usage="{u,v}->earth[u,v] is a parametrization of the earth (with equatorial diameter 12756.4 kilometers and polar diameter 12713.2 kilometers).  To plot try earth[-Pi,Pi,-Pi/2,Pi/2]."

    eighteight::usage="{u,v}->eighteight[u,v] is a parametrization of the generalized surface of revolution formed by moving a figure eight about a figure eight.  To plot try eighteight[-Pi,Pi,-Pi/2,Pi/2]."

    eightsurface::usage="{u,v}->eightsurface[u,v] is a surface of revolution generated by t->eight[t].  To plot try eightsurface[-Pi,Pi,-Pi/2,Pi/2]."

    ellhypcyclide::usage="{u,v}->ellhypcyclide[a,c,k][u,v] is a parametrization of the cyclide of Dupin of radius k whose focal sets are the ellipse u->{a*Cos[u],Sqrt[a^2 - c^2]*Sin[u],0} and the hyperbola v->{c*Sec[v],0,Sqrt[a^2 - c^2]*Tan[v]}.  To plot try ellhypcyclide[8,3,5][-Pi,Pi,-Pi,Pi]."

    ellipsemin::usage="{u,v}->ellipsemin[a,b][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->ellipsemin[a,b][0][u,v] is an elliptical catenoid and {u,v}->ellipsemin[a,b][Pi/2][u,v] is an elliptical helicoid. The patch {u,v}->ellipsemin[a,b][0][u,v] contains the ellipse t->ellipse[a,b][t] as a geodesic.  To plot try ellipsemin[2,1][0][-Pi,Pi,-1,1]."

    ellipsoid::usage="{u,v}->ellipsoid[a,b,c][u,v] is a parametrization of an ellipsoid with axes of lengths a, b and c. {u,v}->ellipsoid[a,c][u,v] is a parametrization of an ellipsoid with axes of lengths a, a and c.  To plot try ellipsoid[1,2,3][-Pi,Pi,-Pi/2,Pi/2]."

    ellipticalhyperboloidminus::usage= "{u,v}->ellipticalhyperboloidminus[a,b,c][u,v] is a parametrization of a hyperboloid of one sheet as a ruled surface.  To plot try ellipticalhyperboloidminus[1,2,3][-Pi,Pi,-Pi/2,Pi/2]."

    ellipticalhyperboloidplus::usage= "{u,v}->ellipticalhyperboloidplus[a,b,c][u,v] is a parametrization of a hyperboloid of one sheet as a ruled surface.  To plot try ellipticalhyperboloidplus[1,2,3][-Pi,Pi,-Pi/2,Pi/2]."

    ellipticcoor::usage="{u,v,w}->ellipticcoor[d,e,f][a,b,c][w][u,v] is an elliptic coordinate patch on R^3. The parameters d, e and f are plus 1 or minus 1.  To plot try ellipticcoor[1,1,1][12,5,1][-1][1,5,5,12]."

    ellipticparaboloid::usage="{u,v}->ellipticparaboloid[a,b][u,v] is the surface {u,v}->{u,v,u^2/a^2 + v^2/b^2}.  To plot try ellipticparaboloid[1,2][-1,1,-1,1]."

    enneper::usage="{u,v}->enneper[u,v] is a parametrization of Enneper's minimal surface.  To plot try enneper[-1.5,1.5,-1.5,1.5]."

    ennepercatenoidpolar::usage="{r,theta}->ennepercatenoidpolar[n,a,b][t][r,theta] is the polar parametrization of a 1-parameter family of minimal surfaces constructed from the minimal curve whose Weierstrass representation is given by f[z_]:=2 and g[z_]:=a/z + b*z^n. For t=0 and a=b=1 the surface has a catenoid-like bottom and an enneper-like top with n + 1 lobes.  To plot try ennepercatenoidpolar[1,1,1][0][0.40,1.08,-Pi,Pi]."

    enneperpolar::usage="{r,theta}->enneperpolar[n][t][r,theta] is the polar parametrization of a 1-parameter family of minimal surfaces such that {r,theta}->enneperpolar[n][0][r,theta] is Enneper's minimal surface of degree n.  To plot try enneperpolar[1][0][0,2,-Pi,Pi,PlotPoints->{15,69}]."

    exptwist::usage="{u,v}->exptwist[a,c][u,v] is a parametrization of a helicoid-like surface whose twisting varies exponentially.  To plot try exptwist[1,0.2][0,4Pi,-1,1,PlotPoints->{60,10}]."

    exptwistasym::usage="{p,q}->exptwistasym[a,c][p,q] is the parametrization by asymptotic curves of {u,v}->exptwist[a,c][u,v].  To plot try exptwistasym[1,0.3][-Pi,2Pi,-Pi,2Pi,PlotPoints->{48,15}]."

    flattorus::usage="{u,v}->flattorus[a,b][u,v] is a parametrization of a torus into R^4 as the product of two circles. {u,v}->flattorus[a,b,c,d][u,v]} is a parametrization of an elliptical torus into R^4 as the product of two ellipses. The Gaussian curvature is zero."

    funnel::usage="{u,v}->funnel[1,1,c][u,v] is a parametrization of the surface of revolution formed by rotating the graph of t->c*Log[t].  To plot try funnel[5,1,2][-Pi,Pi,0.01,1]."

    funnelasym::usage="{p,q}->funnelasym[c][p,q] is the parametrization by asymptotic curves of the surface of revolution formed by rotating the graph of t->c*Log[t].  To plot try funnelasym[1][-Pi/2,Pi/2,0.0,Pi,PlotRange->All]."

    gencube::usage="{u,v}->gencube[n,a,b,c][u,v] is a generalized cube. {u,v}->gencube[1,a/Sqrt[2],a/Sqrt[2],a][u,v] is a sphere. {u,v}->gencube[2,a,a,a][u,v] is a parametrization of an ordinary cube. {u,v}->genoctahedron[3,a,a,a][u,v] is an astroidial cube.  To plot try gencube[2,1,1,1][-Pi/2,Pi/2,-Pi,Pi,PlotRange->All,PlotPoints->{25,25}]."

    genhypparab::usage="{u,v}->genhypparab[a,b][u,v] is a parametrization of a hyperbolic paraboloid as a ruled surface.  To plot try genhypparab[1,1][-1,1,-1,1]."

    genoctahedron::usage="{u,v}->genoctahedron[n,a,b,c][u,v] is a generalized octahedron. {u,v}->genoctahedron[1,a,b,c][u,v] is an ellipsoid. {u,v}->genoctahedron[2,a,a,a][u,v] is a parametrization of an ordinary octahedron. {u,v}->genoctahedron[3,a,b,c][u,v] is an astroidial ellipsoid.  To plot try genoctahedron[2,1,1,1][-Pi/2,Pi/2,-Pi,Pi,PlotRange->All,PlotPoints->{25,25}]."

    handkerchief::usage="{u,v}->handkerchief[a][u,v] is a parametrization of a handkerchief shaped surface.  To plot try handkerchief[1][-1,1,-1,1,PlotRange->All]."

    heart::usage="{u,v}->heart[a][u,v] is a surface of revolution generated by t->cardioid[a][t] with the axis of revolution passing through the cusp.  To plot try heart[1][-Pi/2,Pi/2,-Pi,Pi,PlotPoints->{15,45}]."

    helicoid::usage="{u,v}->helicoid[a,c][u,v] is the standard parametrization of a circular helicoid of radius a and slant c. {u,v}->helicoid[a,b,c][u,v] is the standard parametrization of an elliptical helicoid of slant c.  To plot try helicoid[1,1,0.3][0,4Pi,-1,1,PlotPoints->{75,15}]."

    helicoidprin::usage="{u,v}->helicoidprin[b][u,v] is a parametrization of a circular helicoid of radius b and slant b by principal curves. The parametrization is also isothermal.  To plot try helicoidprin[1][-Pi/2,Pi/2,-Pi/2,Pi/2,PlotRange->All]."

    heltocat::usage="{u,v}->heltocat[t][u,v] is a 1-parameter family of minimal surfaces connecting a helicoid to a catenoid.  To plot try heltocat[Pi/4][0,4Pi,-1.5,1.5,PlotPoints->{75,15}]."

    henneberg::usage="{u,v}->henneberg[u,v] is a parametrization of Henneberg's minimal surface.  To plot try henneberg[0.3,1,-Pi,Pi,PlotPoints->{15,45}]."

    hyperbolicparaboloid::usage="{u,v}->hyperbolicparaboloid[u,v] is a parametrization of a hyperbolic paraboloid.  To plot try hyperbolicparaboloid[-1,1,-1,1]."

    hyperboloid::usage="{u,v}->hyperboloid[a,c][u,v] is the standard parametrization of a circular hyperboloid of one sheet. {u,v}->hyperboloid[a,b,c][u,v] is the standard parametrization of an elliptical hyperboloid of one sheet.  To plot try hyperboloid[1,2,3][0,3Pi/2,-2,2]."

    hypocycloidinvolutemin::usage="{u,v}->hypocycloidinvolutemin[a,b][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->hypocycloidinvolutemin[a,b][0][u,v] is a minimal surface containing the involute of a hypocycloid as a geodesic.  To plot try hypocycloidinvolutemin[4,1][0][-Pi,Pi,-1,1,PlotPoints->{99,15}]."

    hy2sheet::usage="{u,v}->hy2sheet[a,c][u,v] is the standard parametrization of a circular hyperboloid of two sheets. {u,v}->hy2sheet[a,b,c][u,v] is the standard parametrization an elliptical hyperboloid of two sheets.  To plot try hy2sheet[1,1,1][-1,1,-1,1]."

    invertedcircularcone::usage="{u,v}->invertedcircularcone[{p1,p2,p3},rho,a,b][u,v] is a parametrization of {u,v}->circularcone[a,b][u,v] inverted with respect to a sphere of radius rho centered at {p1,p2,p3}. {u,v}->invertedcircularcone[rho,a,b][u,v] is a parametrization of {u,v}->{u,v}->circularcone[a,b][u,v] inverted with respect to a sphere of radius rho centered at the origin. An inverted circular cone is a cyclide of Dupin.  To plot try invertedcircularcone[40,1,1][-Pi,Pi,-40,40,PlotRange->All]."

    invertedcircularcylinder::usage="{u,v}->invertedcircularcylinder[{p1,p2,p3},rho,a][u,v] is a parametrization of {u,v}->circularcylinder[a][u,v] inverted with respect to a sphere of radius rho centered at {p1,p2,p3}. {u,v}->invertedcircularcylinder[rho,a][u,v] is a parametrization of {u,v}->circularcylinder[a][u,v] inverted with respect to a sphere of radius rho centered at the origin. An inverted circular cylinder is a cyclide of Dupin.  To plot try invertedcircularcylinder[40,1][-Pi,Pi,-40,40,PlotRange->All]."

    invertedtorus::usage="{u,v}->invertedtorus[{p1,p2,p3},rho,a,b][u,v] is a parametrization of {u,v}->torus[a,b][u,v] inverted with respect to a sphere of radius rho centered at {p1,p2,p3}. {u,v}->invertedtorus[rho,a,b][u,v] is a parametrization of {u,v}->torus[a,b][u,v] inverted with respect to a sphere of radius rho centered at the origin. An inverted torus is a cyclide of Dupin.  To plot try invertedtorus[40,8,3][-Pi,Pi,-Pi,Pi,PlotRange->All]."

    jorgemeekspolar::usage="jorgemeekspolar[n,a][r,theta] is a polar parametrization of a Jorge-Meeks n-noid.  To plot try jorgemeekspolar[3,1][0.001,1.999,-Pi,Pi]."

    kidney::usage="{u,v}->kidney[a][u,v] is a surface of revolution generated by t->nephroid[a][t].  To plot try kidney[1][0,3Pi/2,-Pi/2,Pi/2]."

    kleinbottle::usage="{u,v}->kleinbottle[a][u,v] is a parametrization in R^3 of the surface formed by moving and twisting a figure eight along a circle of radius a. The surface is nonorientable and a neighborhood of the self-intersection curve is nonorientable.  To plot try kleinbottle[2][-Pi/4,3Pi/2,-Pi,Pi]."

    kleinbottlebis::usage="{u,v}->kleinbottlebis[e,f][u,v] is a Klein bottle in R^3. A neighborhood of the self-intersection curve is orientable.  To plot try kleinbottlebis[1,2][-Pi,29.9517,-Pi,Pi]."

    kleinbottlewri::usage="{u,v}->kleinbottlewri[u,v] is WRI's version of a Klein bottle in R^3. A neighborhood of the self-intersection curve is orientable.  To plot try kleinbottlewri[0,2Pi,0,2Pi]."

    klein4::usage="{u,v}->klein4[a,b][u,v] is a parametrization of a Klein bottle as a surface in R^4."

    kuen::usage="{u,v}->kuen[a][u,v] is a parametrization of Kuen's surface. It has constant negative curvature -1/a^2.  To plot try kuen[1][-4,4,0.01,Pi - 0.01]."

    kuentcheb::usage="{u,v}->kuentcheb[a][u,v] is a Tchebyshef principal patch and {p,q}->kuentcheb[a][p + q,p - q] is a Tchebyshef asymptotic patch for Kuen's surface of constant negative curvature -1/a^2.  To plot try kuentcheb[1][-4,4,-4,4]."

    lemniscatemin::usage="{u,v}->lemniscatemin[a][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->lemniscatemin[a][0][u,v] is a minimal surface containing a lemniscate as a geodesic.  To plot try lemniscatemin[1][0][-Pi,Pi,-0.2Pi,0.2Pi,PlotPoints->{100,20}]."

    logminpolar::usage="{r,theta}->logminpolar[n][t][r_,theta_] is a 1-parameter family of minimal surfaces whose height function is the identity and whose Gauss map is an iterated logarithm.  To plot try logminpolar[0][0][0.005,2.005,-3Pi,3Pi,PlotPoints->{10,60},ViewPoint->{-2.5,-1.5,1.6}]."

    menn::usage="{u,v}->menn[a][u,v] is a parametrization of Menn's surface.  To plot try menn[1][-1.5,1.5,-1.5,1.5]."

    mercatorellipsoid::usage="{u,v}->mercatorellipsoid[a,b,c][u,v] is the Mercator injection of an ellipsoid with axes of lengths a, b and c into R^3.  To plot try mercatorellipsoid[1,2,3][0,3Pi/2,-0.49Pi,0.49Pi]."

    mercatorsphere::usage="{u,v}->mercatorsphere[a][u,v] is the Mercator injection of a sphere of radius a into R^3.  To plot try mercatorsphere[1][0,3Pi/2,-0.49Pi,0.49Pi]."

    moebcirc::usage="{u,v}->moebcirc[a][u,v] is a parametrization of a Moebius strip that has a circle as boundary.  To plot try moebcirc[1][0,Pi,0,Pi]."

    moebiusstrip::usage="{u,v}->moebiusstrip[a][u,v] is the standard parametrization of a Moebius strip.  To plot try moebiusstrip[1][0,2Pi,-0.3,0.3,PlotPoints->{40,5}]."

    monkeysaddle::usage="{u,v}->monkeysaddle[u,v] is a parametrization of a monkey saddle. {u,v}->monkeysaddle[n][u,v] is a parametrization of a monkey saddle for a monkey with n-2 tails.  To plot try monkeysaddle[-1.1,1.1,-1.1,1.1]."

    monkeysaddlepolar::usage="{r,theta}->monkeysaddlepolar[n][r,theta] is a polar parametrization of a monkey saddle with n - 2 tails.  To plot try monkeysaddlepolar[3][0,1,0,2Pi,PlotRange->All,PlotPoints->{10,33}]."

    nodoid::usage="See nodoidnds.  To plot try nodoid[-9,9][8,3,1,0.6][0,3Pi/2,-9,9,PlotPoints->{21,51}]."

    nodoidnds::usage="nodoidnds[vmin,vmax][a,b,c,e][u,v] is an approximation over the interval vmin<v<vmax found by NDSolve to a nodoid. It is a surface of revolution of a curve of Delaunay of eccentricity e and has constant mean curvature."

    parabolamin::usage="{u,v}->parabolamin[a][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->parabolamin[a][0][u,v] is a minimal surface containing a parabola as a geodesic.  To plot try parabolamin[1][0][-2,2,-0.3,0.3,PlotPoints->{40,10}]."

    paraboliccoor::usage="{u,v,w}->paraboliccoor[d,e,f][a,b,c][w][u,v] is a parabolic coordinate patch on R^3. The parameters d, e and f are plus 1 or minus 1.  To plot try paraboliccoor[1,1,1][2,1,3][4][-9,1,1,2]."

    paraboliccyclide::usage="{u,v}->paraboliccyclide[a,k][u,v] is a parametrization of the cyclide of Dupin of radius k whose focal sets are the parabolas u->{u,0,-u^2/(8*a) + a} and v->{0,v,v^2/(8*a) - a}.  To plot try paraboliccyclide[2,1][-40,40,-40,40]."

    paraboloid::usage="{u,v}->paraboloid[a][u,v] is a circular paraboloid of radius a. {u,v}->paraboloid[a,b][u,v] and {u,v}->paraboloid[a,b,c][u,v] is either an elliptic or a hyperbolic paraboloid.  To plot try paraboloid[1][-1,1,-1,1]."

    paraboloidpolar::usage="{r,theta}->paraboloidpolar[n,a,b][r,theta] is a polar parametrization of a generalized paraboloid.  To plot try paraboloidpolar[2,2,1][0,1,0,3Pi/2]."

    pear::usage="{u,v}->pear[a,b][u,v] is a surface of revolution generated by t->piriform[a,b][t] with the axis of revolution passing through the cusp.  To plot try pear[1,1][0,3Pi/2,-Pi/2,Pi/2]."

    pearbis::usage="{u,v}->pearbis[u,v] is a parametrization of the surface of revolution that resembles an pear.  To plot try pearbis[0,3Pi/2,-Pi,Pi]."

    perturbedms::usage="{u,v}->perturbedms[a][u,v] is a parametrization of a monkey saddle perturbed by a circular paraboloid.  To plot try perturbedms[-1.0][-1,1,-1,1]."

    perturbedmspolar::usage="{r,theta}->perturbedmspolar[n,a][r,theta] is a polar parametrization of a monkey saddle of order n perturbed by a circular paraboloid.  To plot try perturbedmspolar[3,-1][0,1,0,2Pi,PlotPoints->{10,40}]."
    
    pillow::usage="{u,v}->pillow[a][u,v] is a pillow-shaped surface.  To plot try pillow[1][-Pi,Pi,-Pi,Pi]."

    plane::usage="{u,v}->plane[a1,a2,b1,b2,c1,c2][u,v] is a parametrization of a plane in R^3."

    planepolar::usage="{r,theta}->planepolar[a][r,theta] is a polar parametrization of the plane."

    planetocat::usage="{u,v}->planetocat[t][u,v] is a 1-parameter family of minimal surfaces containing a plane and a catenoid. The principal curves of each surface are planar.  To plot try planetocat[0.1][0,3Pi/2,-2,2]."

    plucker::usage="{u,v}->plucker[u,v] is the standard parametrization of Plucker's surface. {u,v}->plucker[n][u,v] is a parametrization of Plucker's surface with n folds. {u,v}->plucker[m,n,a][u,v] is a generalization of the monkey saddle and Plucker's surface.  To plot try plucker[-1,1,-1,1]."

    pluckerpolar::usage="{r,theta}->pluckerpolar[n][r,theta] is the polar parametrization of Plucker's surface with n folds. {r,theta}->pluckerpolar[m,n,a][r,theta] is a generalization of the polar parametrizations of the monkey saddle and Plucker's surface.  To plot try pluckerpolar[0,5,1][0,1,-Pi,Pi,PlotPoints->{10,75}]."

    pseudocrosscap::usage="{u,v}->pseudocrosscap[u,v] is a parametrization of a pseudocrosscap in R^3.  To plot try pseudocrosscap[-1,0.5,0,2Pi,PlotPoints->{24,48}]."

    pseudosphere::usage="{u,v}->pseudosphere[a][u,v] is the standard parametrization of a pseudosphere of constant negative curvature -1/a^2.  To plot try pseudosphere[1][0,3Pi/2,0.001,Pi-0.001]."

    pseudospherebis::usage="{u,v}->pseudospherebis[a][u,v] is a parametrization of a pseudosphere of constant negative curvature -1/a^2.  To plot try pseudospherebis[1][0,3Pi/2,0.01,2,PlotRange->All]."

    pseudospheretcheb::usage="{u,v}->pseudospheretcheb[a][u,v] is a Tchebyshef principal patch and {p,q}->pseudospheretcheb[a][p + q,p - q] is a Tchebyshef asymptotic patch for a pseudosphere of constant negative curvature -1/a^2.  To plot try pseudospheretcheb[1][-2.8,2.8,0,3Pi/2]."

    rconoid::usage="{u,v}->rconoid[u,v] is a right conoid. It coincides with {u,v}->2*pluckerpolar[1][v/2,u].  To plot try rconoid[-Pi,Pi,-2,2]."

    richmond::usage="{u,v}->richmond[u,v] is a minimal surface with one planar end.  To plot try richmond[-1,1,-1,1,PlotRange->{{-1,1},{-1,1},Automatic}]."

    richmondpolar::usage="{r,theta}->richmondpolar[n][t][r,theta] is the polar parametrization of a 1-parameter family of minimal surfaces such that {r,theta}->richmondpolar[n][0][r,theta] is a minimal surface with one planar end of degree n.  To plot try richmondpolar[1][0][0.3,1.3,0,2Pi]."

    roman::usage="{u,v}->roman[a][u,v] is a parametrization of Steiner's Roman surface of radius a in R^3.  To plot try roman[1][0,Pi,-Pi/2,Pi/2]."

    rosemin::usage="{u,v}->rosemin[n,a][t][u,v] is a 1-parameter family of minimal surfaces such that {u,v}->rosemin[n,a][0][u,v] contains the rose t->{a*Cos[t]*Cos[n*t],a*Cos[n*t]*Sin[t]} as a geodesic.  To plot try rosemin[3,1][0][0,Pi,-0.3,0.3,PlotPoints->{45,15}]."

    scherk::usage="{u,v}->scherk[a][u,v] is a parametrization of Scherk's minimal surface. {u,v}->scherk[a,theta][u,v] is a parametrization of twisted Scherk's minimal surface of angle theta.  To plot try scherk[1][-Pi/2 + 0.01,3Pi/2 + 0.01,-Pi/2 + .01,3Pi/2 + 0.01]."

    scherk5::usage="{u,v}->scherk5[a,b,c][u,v] is a parametrization of Scherk's fifth minimal surface, where a, b, c are plus 1 or minus 1.  {u,v}->scherk5[u,v] is the same as {u,v}->scherk5[a,b,c][u,v]. To plot try scherk5[-1,1,-1,1]."

    shoe::usage="{u,v}->shoe[a,-a][u,v] resembles the instep of a shoe.  To plot try shoe[1,-1][-1.5,1.5,-1.5,1.5]."

    shoeasym::usage="{u,v}->shoeasym[a,b][u,v] is an asymptotic parametrization of a shoe surface.  To plot try shoeasym[1,-1][-1,1.2,-1,1.2]."

    sievert::usage="{u,v}->sievert[a][u,v] is a parametrization of Sievert's surface of constant positive curvature a^2.  To plot try sievert[1][-Pi/2,Pi/2,0.001,Pi - 0.001]."

    sinovaloid::usage="{u,v}->sinovaloid[n,a][u,v] is the surface of revolution generated by t->sinoval[n,a][t].  To plot try sinovaloid[3,1][0,3Pi/2,-Pi/2,Pi/2]."

    sinsurface::usage="{u,v}->sinsurface[a][u,v] is the {u,v}->a*{Sin[u],Sin[v],Sin[u + v]}. The implicit equation is 4*x^2*y^2*(a^2 - z^2) == a^2*(x^2 + y^2 - z^2)^2.  To plot try sinsurface[1][-Pi,Pi,-Pi,Pi,PlotPoints->{40,40}]."

    snail::usage="{u,v}->snail[u,v] is a surface resembling a snail.  To plot try snail[-Pi,Pi,-Pi,Pi,PlotRange->All]."

    sphere::usage="{u,v}->sphere[a][u,v] is the standard parametrization of a sphere of radius a centered at the origin.  To plot try sphere[1][0,3Pi/2,-Pi/2,Pi/2]."

    stereographicellipsoid::usage="{u,v}->stereographicellipsoid[a,b,c][u,v] is a stereographic parametrization of an ellipsoid with axes of lengths a, b and c.  To plot try stereographicellipsoid[5,3,1][-2.8,2.8,-2.8,2.8]."

    stereographicsphere::usage="{u,v}->stereographicsphere[a][u,v] is a stereographic parametrization of a sphere of radius a in R^3.  To plot try stereographicsphere[1][-2.8,2.8,-2.8,2.8,PlotPoints->{40,40}]."

    stereographicspherepolar::usage="{r,theta}->stereographicspherepolar[a][r,theta] is a polar stereographic parametrization of a sphere of radius a in R^3.  To plot try stereographicspherepolar[1][0,3,0,2Pi,PlotPoints->{10,40}]."

    swallowtail::usage="{u,v}->swallowtail[u,v] is a parametrization of a swallow tail shaped surface.  To plot try swallowtail[-3,2,-0.8,0.8,PlotRange->All,PlotPoints->{10,15},ViewPoint->{4,1,-4}]."

    tetrahedral::usage="{u,v}->tetrahedral[m,n][A,B,C][a,b,c][u,v] is a parametrization of a tetrahedral surface."

    thomsen::usage="{u,v}->thomsen[a,b][u,v] is a parametrization of Thomsen's minimal surface.  To plot try thomsen[1,-2][-Pi/2,Pi/2,0,2Pi]."

    torus::usage="{u,v}->torus[a,b][u,v] is a parametrization of the torus formed by revolving a circle of radius b in the xz-plane about the z-axis along a circle of radius a in the xy-plane. {u,v}->torus[a,b,c][u,v] is a parametrization of the torus formed by revolving an ellipse with axes of lengths b and c in the xz-plane about the z-axis along a circle of radius a in the xy-plane. {u,v}->torus[a,b,b,d,d][u,v] is a surface formed by moving a circle of radius b in the xz-plane along an ellipse in the xy-plane. {u,v}->torus[a,b,c,d,e][u,v] is a surface formed by moving an ellipse in the xz-plane along an ellipse in the xy-plane.  To plot try torus[8,3][0,3Pi/2,-Pi,Pi]."

    torusasym::usage="{p,q}->torusasym[pm,a][p,q] is an asymptotic parametrization of {u,v}->torus[a,a][u,v]. The parameter pm is plus or minus 1.  To plot try torusasym[1,1][-Pi,Pi,-Pi,Pi]."

    torusbis::usage="{u,v}->torusbis[pm,a,b][u,v] is a parametrization of a torus. The parameter pm is plus or minus 1.  To plot try torusbis[1,8,3][5,11,0,2Pi,PlotPoints->{10,33}]."

    twiflat::usage="{u,v}->twiflat[pm,a,c][sl][u,v] is a parametrization of a generalized helicoid of zero Gaussian curvature of slant sl. The parameter pm is plus or minus 1.  To plot try twiflat[1,1,1][1][0,7Pi/2,1,5,PlotPoints->{40,15}]." 

    twisphere::usage="{u,v}->twisphere[a,b][u,v] is a parametrization of a generalized helicoid formed by rotating and twisting a semicircle. It is also called a corkscrew surface.  To plot try twisphere[1,1][0,2Pi,-Pi,Pi]." 

    veronese::usage="{u,v}->veronese[a][u,v] is a parametrization of a Veronese surface in R^5."

    wallis::usage="{u,v}->wallis[a,b,c][u,v] is a parametrization of a conical edge of Wallis.  To plot try wallis[1,1,3][Pi/4,3Pi/4,0,1]."

    whitneyumbrella::usage="{u,v}->whitneyumbrella[u,v] is a parametrization of a Whitney umbrella.  To plot try whitneyumbrella[-1,1,-1,1]."

    wrinkles::usage="{u,v}->wrinkles[a][u,v] is the Monge patch {u,v}->{u,v,a*(u/v + v/u)}.  To plot try wrinkles[0.05][-1,1,-1,1]."

(* Implicitly Defined Surfaces *)

    cossurfaceimplicit::usage="cossurfaceimplicit[a][x,y,z]==0 is the nonparametric form of the surface {u,v}->{a*Cos[u],a*Cos[v],a*Cos[u + v]}."

    crosscapimplicit::usage="crosscapimplicit[a,b][x,y,z]==0 is the nonparametric form of a crosscap."

    ellipsoidimplicit::usage="ellipsoidimplicit[a,b,c][x,y,z]==0 is the nonparametric form of an ellipsoid."

    ellipticparaboloidimplicit::usage="ellipticparaboloidimplicit[a,b,c][x,y,z]==0 is the nonparametric form of an elliptic paraboloid."

    equihom1implicit::usage="equihom1implicit[x,y,z]==0 is the nonparametric form of the first equiaffinely homogeneous surface."

    equihom2implicit::usage="equihom2implicit[x,y,z]==0 is the nonparametric form of the second equiaffinely homogeneous surface."

    equihom3implicit::usage="equihom3implicit[x,y,z]==0 is the nonparametric form of the third equiaffinely homogeneous surface."

    equihom4implicit::usage="equihom4implicit[x,y,z]==0 is the nonparametric form of the fourth equiaffinely homogeneous surface."

    equihom5implicit::usage="equihom5implicit[x,y,z]==0 is the nonparametric form of the fifth equiaffinely homogeneous surface."

    equihom6implicit::usage="equihom6implicit[x,y,z]==0 is the nonparametric form of the sixth equiaffinely homogeneous surface."

    goursatimplicit::usage="goursatimplicit[a,b,c][x,y,z]==0 is the nonparametric form of Goursat's surface."

    hyperbolicparaboloidimplicit::usage="hyperbolicparaboloidimplicit[a,b,c][x,y,z]==0 is the nonparametric form of a hyperbolic paraboloid."

    hyperboloidimplicit::usage="hyperboloidimplicit[a,b,c][x,y,z]==0 is the nonparametric form of a hyperboloid of one sheet."

    hy2sheetimplicit::usage="hy2sheetimplicit[a,b,c][x,y,z]==0 is the nonparametric form of a hyperboloid of two sheets."

    kazoolaimplicit::usage="kazoolaimplicit[a,b,c,d][x,y,z]==0 is the nonparametric form of a kazoola."

    kummerimplicit::usage="kummerimplicit[x,y,z]==0 is the nonparametric form of Kummer's surface."

    scherkimplicit::usage="scherkimplicit[a][x,y,z]==0 is the nonparametric form of Scherk's minimal surface. scherkimplicit[a,phi][x,y,z]==0 is the nonparametric form of a twisted Scherk's minimal surface."

    scherk5implicit::usage="scherk5implicit[x,y,z]==0 is the nonparametric form of a Scherk's fifth minimal surface."

    sinsurfaceimplicit::usage="sinsurfaceimplicit[a][x,y,z]==0 is the nonparametric form of the surface {u,v}->{a*Sin[u],a*Sin[v],a*Sin[u + v]}."

    sphereimplicit::usage="sphereimplicit[a][x,y,z]==0 is the nonparametric form of a sphere. sphereimplicit[n,a][x,y,z]==0 is the nonparametric form of the surface x^(2*n) + y^(2*n) +z^(2*n) = a^(2*n)."

    torusimplicit::usage="torusimplicit[a,b][x,y,z]==0 is the nonparametric form of a torus."

    twocuspsimplicit::usage="twocuspsimplicit[x,y,z]==0 is the nonparametric form of a surface with two cusps."

(* Drum Plots *)

    circularmembrane::usage="{r,theta}->circularmembrane[m,n,t][r,theta] is a parametrization of a circular membrane of degree {m,n} at time t."

    circularmembraneplot::usage="circularmembraneplot[m,n,t] plots a circular membrane of degree {m,n} at time t."

    circularmembranemovie::usage=" circularmembranemovie[m,n,{t,tmin,tmax,numf}] animates a circular membrane of degree {m,n} from tmin to tmax with numf frames."

    squaremembrane::usage="{u,v}->squaremembrane[m,n,t][u,v] is a parametrization of a square membrane of degree {m,n} at time t."

    squaremembraneplot::usage="squaremembraneplot[m,n,t] plots a square membrane of degree {m,n} at time t."

    squaremembranemovie::usage="squaremembranemovie[m,n,{t,tmin,tmax,numf}] animates a square membrane of degree {m,n} from tmin to tmax with numf frames."

(* Minimal Curves *)

    astroidmincurve::usage="z->astroidmincurve[a][z] is the minimal curve generated by the astroid t->astroid[a][t]."

    bourmincurve::usage="z->bourmincurve[m][z] is the minimal curve corresponding to Bour's minimal surface of order m. z->bourmincurve[m,n][z] is the minimal curve corresponding to Bour's minimal surface of order {m,n}. Note that bourmincurve[m]=bourmincurve[m,1], enneper[n]=bourmincurve[2,n] and richmond[n][z]=bourmincurve[0,n + 1][z]/2."

    cardioidmincurve::usage="z->cardioidmincurve[a][z] is the minimal curve generated by the cardioid t->cardioid[a][t]."

    catalanmincurve::usage="z->catalanmincurve[a][z] is the minimal curve corresponding to Catalan's minimal surface. It contains a cycloid as a geodesic."

    catenoidmincurve::usage="z->catenoidmincurve[a][z] is the minimal curve corresponding to a catenoid. It contains a circle as a geodesic."

    cayleysexticmincurve::usage="z->cayleysexticmincurve[a][z] is the minimal curve generated by the Cayley's sextic t->cayleysextic[a][t]."

    circleinvolutemincurve::usage="z->circleinvolutemincurve[n,a][z] is the minimal curve generated by the nth involute starting at {a,0} of the circle t->{a*Cos[t],a*Sin[t]}."

    cissoidmincurve::usage="z->cissoidmincurve[a][z] is the minimal curve generated by the cissoid t->cissoid[a][t]."

    clothoidmincurve::usage="z->clothoidmincurve[a][z] is the minimal curve generated by the clothoid t->clothoid[1,a][z]."

    cnchyprofilemincurve::usage="z->cnchyprofilemincurve[pm,a,b][z] is the minimal curve generated by the curve t->cnchyprofile[pm,a,b][t]. The parameter pm is plus 1 or -1."

    costamincurve::usage="z->costamincurve[z] is the minimal curve that determines Costa's minimal surface."

    cpcprofilemincurve::usage="z->cpcprofilemincurve[pm,a,b][z] is the minimal curve generated by the curve t->cpcprofile[pm,a,b][t]. The parameter pm is plus 1 or -1."

    deltoidinvolutemincurve::usage="z->deltoidinvolutemincurve[a][z] is the minimal curve generated by the involute starting at t=0 of t->deltoid[a][t]."

    deltoidmincurve::usage="z->deltoidmincurve[a][z] is the minimal curve generated by the deltoid t->deltoid[a][t]."

    ellipsemincurve::usage="z->ellipsemincurve[a,b][z] is the minimal curve generated by the ellipse t->ellipse[a,b][t]."

    ennepercatenoidmincurve::usage="z->ennepercatenoidmincurve[n,a,b][z] is the minimal curve whose Weierstrass representation is given by f[z_]:=2 and g[z_]:=a/z + b*z^n."

    ennepermincurve::usage="z->ennepermincurve[n][z] is the minimal curve corresponding to Enneper's minimal surface of order n."

    epicycloidmincurve::usage="z->epicycloidmincurve[a,b][z] is the minimal curve generated by the epicycloid t->epicycloid[a,b][t]."

    hennebergmincurve::usage="z->hennebergmincurve[z] is the minimal curve corresponding to Henneberg's minimal surface."

    hyperbolamincurve::usage="z->hyperbolamincurve[a,b][z] is the minimal curve generated by the hyperbola t->hyperbola[a,b][t]."

    hypocycloidinvolutemincurve::usage="z->hypocycloidinvolutemincurve[a,b][z] is the minimal curve generated by the involute starting at t=0 of t->hypocycloid[a,b][t]."

    hypocycloidmincurve::usage="z->hypocycloidmincurve[a,b][z] is the minimal curve generated by the hypocycloid t->hypocycloid[a,b][t]."

    jorgemeeksmincurve::usage="z->jorgemeeksmincurve[n,a][z] is the minimal curve used to generate the Jorge-Meeks n-noid."

    jorgemeeksmincurveprime::usage="z->jorgemeeksmincurveprime[n,a][z] is the derivative of z->jorgemeeksmincurve[n,a][z]."

    lemniscatemincurve::usage="z->lemniscatemincurve[a][z] is the minimal curve generated by the lemniscate t->lemniscate[a][t]."

    logmincurve::usage="z->logmincurve[n][z] is a minimal curve whose height function is the identity and whose Gauss map is an iterated logarithm."

    logspiralmincurve::usage="z->logspiralmincurve[a,b][z] is the minimal curve generated by the log spiral t->logspiral[a,b][t]."

    nephroidmincurve::usage="z->nephroidmincurve[a][z] is the minimal curve generated by the nephroid t->nephroid[a][t]."

    nielsenspiralmincurve::usage="z->nielsenspiralmincurve[a][z] is the minimal curve generated by Nielsen's spiral t->nielsenspiral[a][t]."

    parabolamincurve::usage="z->parabolamincurve[a][z] is the minimal curve generated by the parabola t->parabola[a][t]."

    richmondmincurve::usage="z->richmondmincurve[n][z] is the minimal curve corresponding to Richmond's minimal surface of order n."

    rosemincurve::usage="z->rosemincurve[n,a][z] is the minimal curve generated by the rose t->rose[n,a][t]."

    scherktowermincurveprime::usage="z->scherktowermincurveprime[n][z] is the integrand of the Scherk tower minimal curve."

    scmincurve::usage="z->scmincurve[z] is the minimal curve generated by the semicubical parabola t->sc[t]."

    tractrixmincurve::usage="z->tractrixmincurve[a][z] is the minimal curve generated by the tractrix t->tractrix[a][t]."

(* Surface Metrics *)

    geopolar::usage="geopolar[f] is the geodesic polar coordinate metric relative to the function f."

    liouville::usage="liouville[f,g] is the metric of a Liouville patch corresponding to the functions f and g."

    poincare::usage="poincare[lambda] is the metric on the Poincare upper half-plane with constant curvature -1/lambda^2. poincare[m,n,lambda] is a generalization of the metric poincare[lambda]."

    poindisk::usage="poindisk[lambda] is the metric on the Poincare disk with constant curvature -1/lambda^2."

    schwarzschildplane::usage="schwarzschildplane[m] is the 2-dimensional Schwarzschild metric. schwarzschildplane[m,f] is the 2-dimensional Schwarzschild metric relative to the function f."

Begin["`Private`"]

(* Parametrically Defined Surfaces *)

    agnesirev[a_][u_,v_]:=
    {-2*a*Cos[u]*Cos[v]^2,-2*a*Cos[v]^2*Sin[u],
     2*a*Tan[v]}

    agnesirev[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	agnesirev[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    apple[u_,v_]:= {Cos[u]*(4 + 3.8*Cos[v]),
                    Sin[u]*(4 + 3.8*Cos[v]),
        (Cos[v] + Sin[v]-1)*(1 + Sin[v])*
         Log[1 - Pi*v/10] + 7.5*Sin[v]}

    apple[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	apple[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    astell[a_,b_,c_][u_,v_]:= {a*Cos[u]*Cos[v],
                               b*Sin[u]*Cos[v],c*Sin[v]}^3

    astell[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	astell[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    bohdom[a_,b_,c_][u_,v_]:= {a*Cos[u],
                               a*Sin[u] + b*Cos[v],c*Sin[v]}

    bohdom[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	bohdom[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    bour[u_,v_]:= {(2*u^2 - u^4 - 2*v^2 + 6*u^2*v^2 - v^4)/4,
                    u*v*(-1 - u^2 + v^2),(2/3)*(u^3 - 3*u*v^2)}

    bour[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 bour[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    bour[m_][t_][u_,v_]:= Simplify[ComplexExpand[
        analytictominimal[bourmincurve[m]][t][u,v],
        TargetFunctions->{Re,Im}]]

    bourpolar[m_][t_][r_,theta_]:=
      {r^(m - 1)*Cos[t - (m - 1)*theta]/(m - 1) -
       r^(m + 1)*Cos[t - (m + 1)*theta]/(m + 1),
       r^(m - 1)*Sin[t - (m - 1)*theta]/(m - 1) +
       r^(m + 1)*Sin[t - (m + 1)*theta]/(m + 1),
       2*r^m*Cos[t - m*theta]/m}

    bourpolar[m_,n_][t_][r_,theta_]:=
      {r^(m - 1)*Cos[t - (m - 1)*theta]/(m - 1) -
       r^(m + 2*n - 1)*Cos[t - (m + 2*n - 1)*theta]/(m + 2*n - 1),
       r^(m - 1)*Sin[t - (m - 1)*theta]/(m - 1) +
       r^(m + 2*n - 1)*Sin[t - (m + 2*n - 1)*theta]/(m + 2*n - 1),
       2*r^(m + n - 1)*Cos[t - (m + n - 1)*theta]/(m + n - 1)}

    bourpolar[m_,n_][t_:0][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	bourpolar[m,n][t][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    boy[a_,b_,c_][u_,v_]:=
      {(a/2)*(-Cos[v]^2 + 2*Cos[u]^2*Sin[v]^2 + 
        Cos[v]*Cos[u]*Sin[v]*(-Cos[v]^2 + Cos[u]^2*Sin[v]^2) - 
        Sin[v]^2*Sin[u]^2 + 2*Cos[v]*Sin[v]*Sin[u]*
        (-Cos[v]^2 + Sin[v]^2*Sin[u]^2) - Sin[v]^4*Sin[4*u]/4),
       (Sqrt[3]/2)*b*(-Cos[v]^2 + 
        Cos[v]*Cos[u]*Sin[v]*(Cos[v]^2 - Cos[u]^2*Sin[v]^2) + 
        Sin[v]^2*Sin[u]^2 - (Sin[v]^4*Sin[4*u])/4),
       c*(Cos[v] + Cos[u]*Sin[v] + Sin[u]*Sin[v])*
         (-4*Sin[v]*(Cos[u] - Sin[u])*(-Cos[v] + Cos[u]*Sin[v])*
         (Cos[v] - Sin[u]*Sin[v]) + 
         (Cos[v] + Cos[u]*Sin[v] + Sin[u]*Sin[v])^3)}

    boy[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 boy[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    catalan[a_][u_,v_]:= {a*(u - Sin[u]*Cosh[v]),
                          a*(1 - Cos[u]*Cosh[v]),
                          -4*a*Sin[u/2]*Sinh[v/2]}

    catalan[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 catalan[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    catalandef[a_][t_][u_,v_]:=
        {a*Cos[t]*(u - Cosh[v]*Sin[u]) + 
         a*Sin[t]*(v - Cos[u]*Sinh[v]),
         a*Cos[t]*(1 - Cos[u]*Cosh[v]) + 
         a*Sin[t]*Sin[u]*Sinh[v],
         4*a*(1 - Cos[u/2]*Cosh[v/2])*Sin[t] - 
         4*a*Cos[t]*Sin[u/2]*Sinh[v/2]}

    catalandef[a_][t_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 catalandef[a][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    catenoid[a_][u_,v_]:= {a*Cos[u]*Cosh[v/a],
                           a*Sin[u]*Cosh[v/a],v}

    catenoid[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	catenoid[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    catenoid[a_,b_,c_][u_,v_]:= {a*Cos[u]*Cosh[v],
                                 b*Sin[u]*Cosh[v],c*v}

    catenoidasym[a_][p_,q_]:=
        {a*Cos[(p + q)/2]*Cosh[(q - p)/2],
         a*Sin[(p + q)/2]*Cosh[(q - p)/2],
         a*(q - p)/2}

    catenoidasym[a_][pmin_,pmax_,qmin_,qmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 catenoidasym[a][p,q]],{p,pmin,pmax},{q,qmin,qmax},
	opt,Axes->None];

    catenoidiso[a_][u_,v_]:= {a*Cos[u/a]*Cosh[v/a],
                              a*Sin[u/a]*Cosh[v/a],v}

    catenoidiso[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	catenoidiso[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    circleinvolutemin[1,a_][t_][u_,v_]:=
       {a*Cos[t]*(Cos[u]*Cosh[v] +
           u*Cosh[v]*Sin[u] - v*Cos[u]*Sinh[v]) + 
        a*Sin[t]*(v*Cosh[v]*Sin[u] +
           u*Cos[u]*Sinh[v] - Sin[u]*Sinh[v]),
        a*Sin[t]*(-v*Cos[u]*Cosh[v] + Cos[u]*Sinh[v] + 
           u*Sin[u]*Sinh[v]) + a*Cos[t]*(-u*Cos[u]*Cosh[v] +
           Cosh[v]*Sin[u] - v*Sin[u]*Sinh[v]),
        -a*u*v*Cos[t] + a*(u^2 - v^2)*Sin[t]/2}

    circleinvolutemin[1,a_][t_:0][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	circleinvolutemin[1,a][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    circularcone[a_,b_][u_,v_]:= {a*v*Cos[u],a*v*Sin[u],b*v}

    circularcone[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	circularcone[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    circularcylinder[a_][u_,v_]:= {a*Cos[u],a*Sin[u],v}

    circularcylinder[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	circularcylinder[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    clothoidmin[a_][t_][u_,v_]:=
      {a*Sqrt[Pi]*Cos[t]*Re[FresnelS[(u + I*v)/Sqrt[Pi]]] + 
       a*Sqrt[Pi]*Im[FresnelS[(u + I*v)/Sqrt[Pi]]]*Sin[t],
       a*Sqrt[Pi]*Cos[t]*Re[FresnelC[(u + I*v)/Sqrt[Pi]]] + 
       a*Sqrt[Pi]*Im[FresnelC[(u + I*v)/Sqrt[Pi]]]*Sin[t],
       -(a*v*Cos[t]) + a*u*Sin[t]}

    clothoidmin[a_][t_:0][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	clothoidmin[a][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    cnccosurfrev[a_,b_][u_,v_]:= {b*Cos[u]*Sinh[v/a],
                                  b*Sin[u]*Sinh[v/a],
                      -I*Sqrt[a^2 - b^2]*EllipticE[I*v/a,
                      -b^2/(a^2 - b^2)]}

    cnccosurfrev[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 cnccosurfrev[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    cnccotcheb[a_,b_][u_,v_]:=
        {(a/b)*JacobiCN[b*u,1/b^2]*Cos[b*v],
         (a/b)*JacobiCN[b*u,1/b^2]*Sin[b*v],
         a*(b*u - EllipticE[JacobiAmplitude[b*u,1/b^2],1/b^2])}
         
    cnccotcheb[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 cnccotcheb[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    cnccotcheb[a_,b_,c_][u_,v_]:=
        {(a/b)*JacobiCN[b*u,1/b^2]*Cos[b*v],
         (a/b)*JacobiCN[b*u,1/b^2]*Sin[b*v],
         c*(b*u - EllipticE[JacobiAmplitude[b*u,1/b^2],1/b^2])}

    cnchysurfrev[a_,b_][u_,v_]:= {b*Cos[u]*Cosh[v/a],
                                  b*Sin[u]*Cosh[v/a],
                           - I*a*EllipticE[I*v/a,-b^2/a^2]}

    cnchysurfrev[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 cnchysurfrev[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    cnchytcheb[a_,b_][u_,v_]:=
        {(a/b)*JacobiDN[u,b^2]*Cos[b*v],
         (a/b)*JacobiDN[u,b^2]*Sin[b*v],
         (a/b)*(u - EllipticE[JacobiAmplitude[u,b^2],b^2])}
         
    cnchytcheb[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	cnchysurfrev[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    cnchytcheb[a_,b_,c_][u_,v_]:=
        {(a/b)*JacobiDN[u,b^2]*Cos[b*v],
         (a/b)*JacobiDN[u,b^2]*Sin[b*v],
         (c/b)*(u - EllipticE[JacobiAmplitude[u,b^2],b^2])}

    cossurface[a_][u_,v_]:= {a*Cos[u],a*Cos[v],a*Cos[u + v]}

    cossurface[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	cossurface[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    costa[u_,v_]:=
        {(1/2)*Re[-WeierstrassZeta[u + I*v,{189.07272,0}] +
        3.14159*u + 0.35889 +
        0.228473*(WeierstrassZeta[u + I*v - 1/2,{189.07272,0}] -
             WeierstrassZeta[u + I*v - I/2,{189.07272,0}])],
        (1/2)*Re[-I*WeierstrassZeta[u + I*v,{189.07272,0}] +
        3.14159*v + 0.35889 -
        0.228473*(I*WeierstrassZeta[u + I*v - 1/2,{189.07272,0}] -
             I*WeierstrassZeta[u + I*v - I/2,{189.07272,0}])],
        0.626657*Log[Abs[(WeierstrassP[u + I*v,{189.07272,0}] -
             6.87519)/
             (WeierstrassP[u + I*v,{189.07272,0}] + 6.87519)]]}

    costa[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	costa[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    cpcsurfrev[a_,b_][u_,v_]:= {b*Cos[u]*Cos[v/a],
                                b*Sin[u]*Cos[v/a],
                                a*EllipticE[v/a,b^2/a^2]}

    cpcsurfrev[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	cpcsurfrev[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    crosscap[a_][u_,v_]:= {a^2*Sin[u]*Sin[2*v]/2,
                 a^2*Sin[2*u]*Cos[v]^2,a^2*Cos[2*u]*Cos[v]^2}

    crosscap[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	crosscap[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    deltoidinvolutemin[a_][t_][u_,v_]:=
          {(a/3)*Cos[t]*(8*Cos[u/2]*Cosh[v/2] + 2*Cos[u]*Cosh[v] - 
            Cos[2*u]*Cosh[2*v]) + 
          (a/3)*Sin[t]*(-8*Sin[u/2]*Sinh[v/2] - 2*Sin[u]*Sinh[v] + 
            Sin[2*u]*Sinh[2*v]),
           (Cos[t]*(-8*a*Cosh[v/2]*Sin[u/2] + 2*a*Cosh[v]*Sin[u] + 
            a*Cosh[2*v]*Sin[2*u]))/3 + 
           (Sin[t]*(-8*a*Cos[u/2]*Sinh[v/2] + 2*a*Cos[u]*Sinh[v] + 
            a*Cos[2*u]*Sinh[2*v]))/3,
           (4*a*Sin[t]*(3*u - 2*Cosh[3*v/2]*Sin[3*u/2]))/9 - 
           (4*a*Cos[t]*(3*v - 2*Cos[3*u/2]*Sinh[3*v/2]))/9}

    deltoidinvolutemin[a_][t_:0][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	deltoidinvolutemin[a][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    deltoidmin[a_][t_][u_,v_]:=
          {a*Cos[t]*(2*Cos[u]*Cosh[v] + Cos[2*u]*Cosh[2*v]) + 
           2*a*Sin[t]*Sin[u]*(-Sinh[v] - Cos[u]*Sinh[2*v]),
           -2*a*(Cos[u]*Sinh[v]*((-1 + Cos[u]*Cosh[v])*Sin[t] + 
            2*a*Cos[t]*Sin[u]*Sinh[v])) + 
            Cosh[v]*Sin[u]*(2*a*Cos[t]*(1 - Cos[u]*Cosh[v]) + 
            2*a*Sin[t]*Sin[u]*Sinh[v]),
            -(8/3)*a*Cos[3*u/2]*Cosh[3*v/2]*Sin[t] - 
           (8/3)*a*Cos[t]*Sin[3*u/2]*Sinh[3*v/2]}

    deltoidmin[a_][t_:0][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	deltoidmin[a][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    dini[a_,b_][u_,v_]:= {a*Cos[u]*Sin[v],
                          a*Sin[u]*Sin[v],
                          a*(Cos[v] + Log[Tan[v/2]]) + b*u}

    dini[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	dini[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    dinihelicoid[pm_,a_,b_,c_][u_,v_]:=
          {v*Cos[u],v*Sin[u],c*u +
           pm*(Integrate[Sqrt[(b*c^2 + a^2*vv^2 + b*vv^2 - c^2*vv^2 - vv^4)/
                          (vv^4 -b*vv^2)],vv])} /. vv->v

    dinitcheb[a_,b_][u_,v_]:=
          {a*Cos[b]*Cos[v]/Cosh[(u - v*Sin[b])/Cos[b]],
           a*Cos[b]*Sin[v]/Cosh[(u - v*Sin[b])/Cos[b]],
           a*(u - Cos[b]*Tanh[(u - v*Sin[b])/Cos[b]])}

    dinitcheb[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	dinitcheb[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    dinitcheb[a_,b_,c_][u_,v_]:=
          {a*Cos[c]*Cos[v]/Cosh[(u - v*Sin[c])/Cos[c]],
           a*Cos[c]*Sin[v]/Cosh[(u - v*Sin[c])/Cos[c]],
           b*(u - Cos[c]*Tanh[(u - v*Sin[c])/Cos[c]])}

    earth[u_,v_]:= {6378.2*Cos[u]*Cos[v],6378.2*Sin[u]*Cos[v],
                    6356.6*Sin[v]}

    earth[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	earth[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    eighteight[u_,v_]:= {Sin[u]*Sin[v],
                         Cos[u]*Sin[u]*Sin[v],Cos[v]*Sin[v]}

    eighteight[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	eighteight[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    eightsurface[u_,v_]:=
        {Cos[u]*Sin[2*v],Sin[u]*Sin[2*v],Sin[v]}

    eightsurface[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	eightsurface[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    ellhypcyclide[a_,c_,k_][u_,v_]:=
        {(c*(k - c*Cos[u]) + a*Cos[u]*(a - k*Cos[v]))/
                 (a - c*Cos[u]*Cos[v]),
         (Sqrt[a^2 - c^2]*(a - k*Cos[v])*Sin[u])/
                 (a - c*Cos[u]*Cos[v]),
         (Sqrt[a^2 - c^2]*(k - c*Cos[u])*Sin[v])/
                 (a - c*Cos[u]*Cos[v])}

    ellhypcyclide[a_,c_,k_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	ellhypcyclide[a,c,k][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    ellipsemin[a_,b_][t_][u_,v_]:=
        {a*Cos[t]*Cos[u]*Cosh[v] - a*Sin[t]*Sin[u]*Sinh[v],
         b*Cos[t]*Sin[u]*Cosh[v] + b*Sin[t]*Cos[u]*Sinh[v],
        -b*(Cos[t]*Im[EllipticE[u + I*v,1 - a^2/b^2]] - 
            Sin[t]*Re[EllipticE[u + I*v,1 - a^2/b^2]])}

    ellipsemin[a_,b_][t:0][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	ellipsemin[a,b][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    ellipsoid[a_,c_][u_,v_]:= {a*Cos[v]*Cos[u],
                               a*Cos[v]*Sin[u],c*Sin[v]}
                              
    ellipsoid[a_,b_,c_][u_,v_]:= {a*Cos[v]*Cos[u],
                                  b*Cos[v]*Sin[u],c*Sin[v]}

    ellipsoid[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	ellipsoid[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    ellipticalhyperboloidminus[a_,b_,c_][u_,v_]:=
        {a*Cos[u] + a*v*Sin[u],-b*v*Cos[u] + b*Sin[u],c*v}

    ellipticalhyperboloidminus[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	ellipticalhyperboloidminus[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    ellipticalhyperboloidplus[a_,b_,c_][u_,v_]:=
        {a*Cos[u] - a*v*Sin[u],b*v*Cos[u] + b*Sin[u],c*v}

    ellipticalhyperboloidplus[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	ellipticalhyperboloidplus[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    ellipticcoor[d_,e_,f_][a_,b_,c_][w_][u_,v_]:=
         {d*Sqrt[(a - u)*(a - v)*(a - w)/((a - b)*(a - c))],
          e*Sqrt[(b - u)*(b - v)*(b - w)/((b - a)*(b - c))],
          f*Sqrt[(c - u)*(c - v)*(c - w)/((c - a)*(c - b))]}

    ellipticcoor[d_,e_,f_][a_,b_,c_][w_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    ellipticcoor[d,e,f][a,b,c][w][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    ellipticparaboloid[a_,b_][u_,v_]:=
                       {u,v,u^2/a^2 + v^2/b^2}

    ellipticparaboloid[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	ellipticparaboloid[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    enneper[u_,v_]:= {u - u^3/3 + u*v^2,
                     -v + v^3/3 - v*u^2,u^2 - v^2}

    enneper[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	enneper[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
	
    enneper[n_][t_][u_,v_]:= Simplify[ComplexExpand[
        analytictominimal[ennepermincurve[n]][t][u,v],
        TargetFunctions->{Re,Im}]]

    ennepercatenoidpolar[n_][t_][r_,theta_]:=
        {(n*(2*n + 1)*r^2*Cos[t - theta] +
          n*(2*n + 1)*Cos[t + theta] - 
          r^(n + 1)*(n*r^(n + 1)*Cos[t - theta - 2*n*theta] + 
          2*Cos[t - n*theta] + 4*n*Cos[t - n*theta]))/
         (n*(2*n + 1)*r),
         (n*(2*n + 1)*r^2*Sin[t - theta] -
          n*(2*n + 1)*Sin[t + theta] +
          r^(n + 1)*(2*Sin[t - n*theta] + 4*n*Sin[t - n*theta] + 
          n*r^(n + 1)*Sin[t + theta - 2*(n + 1)*theta]))/
         (n*(2*n + 1)*r),
          2*(r^(n + 1)*Cos[t - (n + 1)*theta] + 
         (n + 1)*(Cos[t]*Log[r] + theta*Sin[t]))/(n + 1)}

    ennepercatenoidpolar[n_,a_,b_][t_][r_,theta_]:=
        {(n*(2*n + 1)*r^2*Cos[t - theta] + 
          a^2*n*(2*n + 1)*Cos[t + theta] - 
          b*r^(n + 1)*(2*a*Cos[t - n*theta] + 
          4*a*n*Cos[t - n*theta] + 
          b*n*r^(n + 1)*Cos[t - (2*n + 1)*theta]))/
         (n*(2*n + 1)*r),
         (n*(2*n + 1)*r^2*Sin[t - theta] - 
          a^2*n*(2*n + 1)*Sin[t + theta] + 
          b*r^n*(b*n*r^(n + 2)*Sin[t - theta - 2*n*theta] + 
          2*a*r*Sin[t + theta - (n + 1)*theta] + 
          4*a*n*r*Sin[t + theta - (n + 1)*theta]))/
         (n*(2*n + 1)*r),
         (2*b*r^(n + 1)*Cos[t - (n + 1)*theta] + 
          2*a*(n + 1)*Cos[t]*Log[r] + 2*a*(n + 1)*theta*Sin[t])
          /(n + 1)}

    ennepercatenoidpolar[n_,a_,b_][t_:0][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	ennepercatenoidpolar[n,a,b][t][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    enneperpolar[n_][t_][r_,theta_]:=
        {r*Cos[t - theta] - r^(2*n + 1)*
           Cos[t - (2*n + 1)*theta]/(2*n + 1),
         r*Sin[t - theta] + r^(2*n + 1)*
           Sin[t - (2*n + 1)*theta]/(2*n + 1),
         2*r^(n + 1)*Cos[t - (n + 1)*theta]/(n + 1)}

    enneperpolar[n_][t_:0][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	enneperpolar[n][t][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    exptwist[a_,c_][u_,v_]:=
             {a*v*Cos[u],a*v*Sin[u],a*E^(c*u)}

    exptwist[a_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	exptwist[a,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    exptwistasym[a_,c_][p_,q_]:=
            {a*E^(c*(p - q)/2)*Cos[p],
             a*E^(c*(p - q)/2)*Sin[p],
             a*E^(c*p)}

    exptwistasym[a_,c_][pmin_,pmax_,qmin_,qmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 exptwistasym[a,c][p,q]],{p,pmin,pmax},{q,qmin,qmax},
	opt,Axes->None];

    flattorus[a_,b_][u_,v_]:= {a*Cos[u],a*Sin[u],
                               b*Cos[v],b*Sin[v]}

    flattorus[a_,b_,c_,d_][u_,v_]:= {a*Cos[u],b*Sin[u],
                                     c*Cos[v],d*Sin[v]}

    funnel[a_,c_][u_,v_]:= {a*v*Cos[u],
                               c*v*Sin[u],c*Log[v]}

    funnel[a_,b_,c_][u_,v_]:= {a*v*Cos[u],
                               b*v*Sin[u],c*Log[v]}

    funnel[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	funnel[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    funnelasym[c_][p_,q_]:= {c*E^(p - q)*Cos[p + q],
                     -c*E^(p - q)*Sin[p + q],c*(p - q)}

    funnelasym[c_][pmin_,pmax_,qmin_,qmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 funnelasym[c][p,q]],{p,pmin,pmax},{q,qmin,qmax},
	opt,Axes->None];

    gencube[n_,a_,b_,c_][u_,v_]:=
        {a*(Cos[u]*Sqrt[Cos[u]^2]^(n - 1) +
            Sin[u]*Sqrt[Sin[u]^2]^(n - 1))*
           (Cos[v]*Sqrt[Cos[v]^2]^(n - 1) +
            Sin[v]*Sqrt[Sin[v]^2]^(n - 1)),
         b*(-Cos[u]*Sqrt[Cos[u]^2]^(n - 1) +
            Sin[u]*Sqrt[Sin[u]^2]^(n - 1))*
           (Cos[v]*Sqrt[Cos[v]^2]^(n - 1) +
            Sin[v]*Sqrt[Sin[v]^2]^(n - 1)),
         c*(-Cos[v]*Sqrt[Cos[v]^2]^(n - 1) +
            Sin[v]*Sqrt[Sin[v]^2]^(n - 1))}

    gencube[n_,a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	gencube[n,a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    genhypparab[a_,b_][u_,v_]:= {a*(u + v),b*v,u^2 + 2*u*v}

    genhypparab[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	genhypparab[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
	
    genoctahedron[n_,a_,b_,c_][u_,v_]:=
        {a*Cos[u]*Cos[v]*Sqrt[Cos[u]^2*Cos[v]^2]^(n - 1),
         a*Sin[u]*Cos[v]*Sqrt[Sin[u]^2*Cos[v]^2]^(n - 1),
         a*Sin[v]*Sqrt[Sin[v]^2]^(n - 1)}

    genoctahedron[n_,a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	genoctahedron[n,a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    handkerchief[a_][u_,v_]:= {u,v,
                    (1/3)*u^3 + u*v^2 +a*(u^2 - v^2)}

    handkerchief[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	handkerchief[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
	
    heart[a_][u_,v_]:= {2*a*Cos[u]*(1 + Cos[v])*Sin[v],
                        2*a*Sin[u]*(1 + Cos[v])*Sin[v],
                        -2*a*Cos[v]*(1 + Cos[v])}
 
    heart[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	heart[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
	
    helicoid[a_,c_][u_,v_]:= {a*v*Cos[u],a*v*Sin[u],c*u}

    helicoid[a_,b_,c_][u_,v_]:=
         {a*v*Cos[u],b*v*Sin[u],c*u}

    helicoid[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	helicoid[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
	
    helicoidprin[b_][p_,q_]:=  {b*Cos[p + q]*Sinh[-p + q],
         b*Sin[p + q]*Sinh[-p + q],b*(p + q)}

    helicoidprin[b_][pmin_,pmax_,qmin_,qmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 helicoidprin[b][p,q]],{p,pmin,pmax},{q,qmin,qmax},
	opt,Axes->None];

    heltocat[t_][u_,v_]:=
        {Cos[u]*Cosh[v]*Sin[t] + Cos[t]*Sin[u]*Sinh[v],
         Cosh[v]*Sin[t]*Sin[u] - Cos[t]*Cos[u]*Sinh[v],
         u*Cos[t] + v*Sin[t]}

    heltocat[t_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	heltocat[t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    henneberg[u_,v_]:=
        {2*Sinh[u]*Cos[v] - (2/3)*Sinh[3*u]*Cos[3*v],
         2*Sinh[u]*Sin[v] + (2/3)*Sinh[3*u]*Sin[3*v],
         2*Cosh[2*u]*Cos[2*v]}

    henneberg[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	henneberg[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    hyperbolicparaboloid[u_,v_]:= {u,v,u*v}

    hyperbolicparaboloid[a_][u_,v_]:= {u,v,a*u*v}
    
    hyperbolicparaboloid[umin_,umax_,vmin_,vmax_,opt___]:=
	ParametricPlot3D[Evaluate[
	hyperbolicparaboloid[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    hyperboloid[a_,c_][u_,v_]:= {a*Cosh[v]*Cos[u],
                                 a*Cosh[v]*Sin[u],
                                 c*Sinh[v]}

    hyperboloid[a_,b_,c_][u_,v_]:= {a*Cosh[v]*Cos[u],
                                    b*Cosh[v]*Sin[u],
                                    c*Sinh[v]}

    hyperboloid[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	hyperboloid[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    hyperboloidbis[a_,c_][u_,v_]:= {a*Sec[v]*Cos[u],
                                    a*Sec[v]*Sin[u],
                                    c*Tan[v]}

    hyperboloidbis[a_,b_,c_][u_,v_]:= {a*Sec[v]*Cos[u],
                                       b*Sec[v]*Sin[u],
                                       c*Tan[v]}

    hyperboloidbis[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	hyperboloidbis[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    hypocycloidinvolutemin[a_,b_][t_][u_,v_]:=
        {Cos[t]*((a - b)*Cos[u]*Cosh[v] + 
         b*Cos[(a - b)*u/b]*Cosh[(a - b)*v/b]) + 
         Sin[t]*((-a + b)*Sin[u]*Sinh[v] - 
         b*Sin[(a - b)*u/b]*Sinh[(a - b)*v/b]),
         Cos[t]*((a - b)*Cosh[v]*Sin[u] - 
         b*Cosh[(a - b)*v/b]*Sin[(a - b)*u/b]) + 
         Sin[t]*((a - b)*Cos[u]*Sinh[v] - 
         b*Cos[((a - b)*u)/b]*Sinh[(a - b)*v/b]),
         (-4*(a - b)*b*Cos[(a*u)/(2*b)]*Cosh[(a*v)/(2*b)]*Sin[t])/a - 
        (4*(a - b)*b*Cos[t]*Sin[(a*u)/(2*b)]*Sinh[(a*v)/(2*b)])/a}

    hypocycloidinvolutemin[a_,b_][t_:0][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	hypocycloidinvolutemin[a,b][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    hy2sheet[a_,c_][u_,v_]:= {a*Cosh[u]*Cosh[v],
                              a*Sinh[u]*Cosh[v],
                              c*Sinh[v]}

    hy2sheet[a_,b_,c_][u_,v_]:= {a*Cosh[u]*Cosh[v],
                                 b*Sinh[u]*Cosh[v],
                                 c*Sinh[v]}

    hy2sheet[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	hy2sheet[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    invertedcircularcone[{p1_,p2_,p3_},rho_,a_,b_][u_,v_]:=
        {p1 + rho^2*(-p1 + a*v*Cos[u])/
         ((-p3 + b*v)^2 + (-p1 + a*v*Cos[u])^2 + 
          (-p2 + a*v*Sin[u])^2),
         p2 + rho^2*(-p2 + a*v*Sin[u])/
         ((-p3 + b*v)^2 + (-p1 + a*v*Cos[u])^2 + 
          (-p2 + a*v*Sin[u])^2),
         p3 + rho^2*(-p3 + b*v)/
         ((-p3 + b*v)^2 + (-p1 + a*v*Cos[u])^2 + 
          (-p2 + a*v*Sin[u])^2)}

    invertedcircularcone[rho_,a_,b_][u_,v_]:=
        {a*rho^2*Cos[u]/(a^2*v + b^2*v),
         a*rho^2*Sin[u]/(a^2*v + b^2*v),
         b*rho^2/(a^2*v + b^2*v)}

    invertedcircularcone[{p1_,p2_,p3_},rho_,a_,b_][umin_,umax_,
                                                   vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	invertedcircularcone[{p1,p2,p3},rho,a,b][u,v]],{u,umin,umax},
                {v,vmin,vmax},opt,Axes->None];
     
    invertedcircularcylinder[{p1_,p2_,p3_},rho_,a_][u_,v_]:=
        {p1 + rho^2*(-p1 + a*Cos[u])/
         ((-p3 + v)^2 + (-p1 + a*Cos[u])^2 + 
          (-p2 + a*Sin[u])^2),
          p2 + rho^2*(-p2 + a*Sin[u])/
          ((-p3 + v)^2 + (-p1 + a*Cos[u])^2 + 
           (-p2 + a*Sin[u])^2),
          p3 + rho^2*(-p3 + v)/
          ((-p3 + v)^2 + (-p1 + a*Cos[u])^2 + 
           (-p2 + a*Sin[u])^2)}

    invertedcircularcylinder[rho_,a_][u_,v_]:=
        {a*rho^2*Cos[u]/(a^2 + v^2),
         a*rho^2*Sin[u]/(a^2 + v^2),
         rho^2*v/(a^2 + v^2)}

    invertedcircularcylinder[{p1_,p2_,p3_},rho_,a_][umin_,umax_,
                                                    vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    invertedcircularcylinder[{p1,p2,p3},rho,a][u,v]],{u,umin,umax},
                {v,vmin,vmax},opt,Axes->None];
     
    invertedtorus[{p1_,p2_,p3_},rho_,a_,b_][u_,v_]:=
        {p1 + (rho^2*(-p1 + Cos[u]*(a + b*Cos[v])))/
         ((-p1 + Cos[u]*(a + b*Cos[v]))^2 + 
          (-p2 + (a + b*Cos[v])*Sin[u])^2 + (-p3 + b*Sin[v])^2),
         p2 + (rho^2*(-p2 + (a + b*Cos[v])Sin[u]))/
         ((-p1 + Cos[u]*(a + b*Cos[v]))^2 + 
          (-p2 + (a + b*Cos[v])*Sin[u])^2 + (-p3 + b*Sin[v])^2),
         p3 + (rho^2*(-p3 + b*Sin[v]))/
         ((-p1 + Cos[u]*(a + b*Cos[v]))^2 + 
          (-p2 + (a + b*Cos[v])*Sin[u])^2 + (-p3 + b*Sin[v])^2)}

    invertedtorus[rho_,a_,b_][u_,v_]:=
        {rho^2*Cos[u]*(a + b*Cos[v])/(a^2 + b^2 + 2*a*b*Cos[v]),
         rho^2*(a + b*Cos[v])*Sin[u]/(a^2 + b^2 + 2*a*b*Cos[v]),
         b*rho^2*Sin[v]/(a^2 + b^2 + 2*a*b*Cos[v])}

    invertedtorus[{p1_,p2_,p3_},rho_,a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    invertedtorus[{p1,p2,p3},rho,a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    jorgemeekspolar[n_,a_][r_,theta_]:=
        Re[jorgemeeksmincurve[n,a][r*E^(I*theta)]]

    jorgemeekspolar[2,a_][r_,theta_]:=
        {(a/8)*Log[(1 + r^2 + 2*r*Cos[theta])/(1 + r^2 - 2*r*Cos[theta])],
          -a*r*(1 + r^2)*Sin[theta]/(2*(1 + r^4 - 2*r^2*Cos[2*theta])),
          (a - a*r^4)/(4*(1 + r^4 - 2*r^2*Cos[2*theta]))}

    jorgemeekspolar[2,a_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	jorgemeekspolar[2,a][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    jorgemeekspolar[3,a_][r_,theta_]:=
    {(a/18)*(3*r^2/(1 + r^2 + r^4 + 2*(r + r^3)*Cos[theta] + 
        2*r^2*Cos[2*theta]) + (3*(r + r^3)*Cos[theta])/
       (1 + r^2 + r^4 + 2*(r + r^3)*Cos[theta] + 2*r^2*Cos[2*theta]) - 
       2*Log[1 + r^2 - 2*r*Cos[theta]] + 
       Log[1 + r^2 + r^4 + 2*(r + r^3)*Cos[theta] + 2*r^2*Cos[2*theta]]),
     -(a*(-(Sqrt[3]*Log[1 + r^2 + r*Cos[theta] - Sqrt[3]*r*Sin[theta]]) - 
       Sqrt[3]*r^6*Log[1 + r^2 + r*Cos[theta] - 
       Sqrt[3]*r*Sin[theta]] + 2*Sqrt[3]*r^3*Cos[3*theta]*
       (Log[1 + r^2 + r*Cos[theta] - Sqrt[3]*r*Sin[theta]] - 
        Log[1 + r^2 + r*Cos[theta] + Sqrt[3]*r*Sin[theta]]) + 
        Sqrt[3]*Log[1 + r^2 + r*Cos[theta] + Sqrt[3]*r*Sin[theta]] + 
        Sqrt[3]*r^6*Log[1 + r^2 + r*Cos[theta] + Sqrt[3]*r*Sin[theta]] + 
        3*(r + r^5)*Sin[theta] + 3*r^2*Sin[2*theta] + 
        3*r^4*Sin[2*theta]))/(18*(1 + r^6 - 2*r^3*Cos[3*theta])),
     (a - a*r^6)/(6*(1 + r^6 - 2*r^3*Cos[3*theta]))}

    jorgemeekspolar[3,a_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	jorgemeekspolar[3,a][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    kidney[a_][u_,v_]:=
        {a*Cos[u]*(3*Cos[v] - Cos[3*v]),
         a*Sin[u]*(3*Cos[v] - Cos[3*v]),
         a*(3*Sin[v] - Sin[3*v])}

    kidney[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    kidney[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    kleinbottle[a_][u_,v_]:=
        {(a + Cos[u/2]*Sin[v] - Sin[u/2]*Sin[2*v])*Cos[u],
         (a + Cos[u/2]*Sin[v] - Sin[u/2]*Sin[2*v])*Sin[u],
          Sin[u/2]*Sin[v] + Cos[u/2]*Sin[2*v]}

    kleinbottle[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    kleinbottle[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    kleinbottlebis[e_,f_][u_,v_]:=If[u<0,
        {f*Sin[u]/2,(-e - f/2 + f*Cos[u]/2)*Cos[v],
         (e + f/2 - f*Cos[u]/2)*Sin[v]},
        {10*Sin[0.105*u] + 2*(e + 0.033*f*u)*
         (Cos[0.105*u] + 4*Cos[0.21*u])*Cos[v]*
         (Sin[0.105*u] + 2*Sin[0.21*u])/
         Sqrt[100*Cos[0.105*u]^2 + 4*(Cos[0.105*u] + 
          4*Cos[0.21*u])^2*(Sin[0.105*u] + 2*Sin[0.21*u])^2],
         (Sin[0.105*u] + 2*Sin[0.21*u])^2 - 
         (10*(e + 0.033*f*u)*Cos[0.105*u]*Cos[v])/
         Sqrt[100*Cos[0.105*u]^2 + 4*(Cos[0.105*u] + 
          4*Cos[0.21*u])^2*(Sin[0.105*u] + 2*Sin[0.21*u])^2],
         (e + 0.033*f*u)*Sin[v]}]

    kleinbottlebis[e_,f_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    kleinbottlebis[e,f][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    kleinbottlebis[e_,f_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    kleinbottlewri[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    kleinbottlewri[u_,v_]:=If[Pi < u <= 2 Pi,
        {6*Cos[u]*(1 + Sin[u]) + 4*(1 - Cos[u]/2)*Cos[v + Pi],
        16*Sin[u],4*(1 - Cos[u]/2)*Sin[v]},
        {6*Cos[u]*(1 + Sin[u]) + 4*(1 - Cos[u]/2)*Cos[u]*Cos[v],
 	 16*Sin[u] + 4*(1 - Cos[u]/2)*Sin[u]*Cos[v],
 	 4*(1 - Cos[u]/2)*Sin[v]}]

    kleinbottlewri[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[
    kleinbottlewri[u,v],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    klein4[a_,b_][u_,v_]:=
        {(a + b*Cos[v])*Cos[u],(a + b*Cos[v])Sin[u],
          b*Sin[v]*Cos[u/2],b*Sin[v]*Sin[u/2]}

    kuen[a_][u_,v_]:=
         {2*a*(Cos[u] + u*Sin[u])*Sin[v]/(1 + u^2*Sin[v]^2),
          2*a*(Sin[u] - u*Cos[u])Sin[v]/(1 + u^2*Sin[v]^2),
          a*(Log[Tan[v/2]] + 2*Cos[v]/(1 + u^2*Sin[v]^2))}

    kuen[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    kuen[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    kuentcheb[a_][u_,v_]:=
        {2*a*Cosh[u]*(Cos[v] + v*Sin[v])/(v^2 + Cosh[u]^2),
         2*a*Cosh[u]*(Sin[v] - v*Cos[v])/(v^2 + Cosh[u]^2),
          a*(u - Sinh[2*u]/(v^2 + Cosh[u]^2))}

    kuentcheb[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    kuentcheb[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    kuentcheb[a_,b_][u_,v_]:=
        {2*a*Cosh[u]*(Cos[v] + v*Sin[v])/(v^2 + Cosh[u]^2),
         2*a*Cosh[u]*(Sin[v] - v*Cos[v])/(v^2 + Cosh[u]^2),
          b*(u - Sinh[2*u]/(v^2 + Cosh[u]^2))}

    kuentcheb[a_,b_,c_][u_,v_]:=
        {2*a*Cosh[u]*(Cos[v] + v*Sin[v])/(v^2 + Cosh[u]^2),
         2*b*Cosh[u]*(Sin[v] - v*Cos[v])/(v^2 + Cosh[u]^2),
          c*(u - Sinh[2*u]/(v^2 + Cosh[u]^2))}

    lemniscatemin[a_][t_][u_,v_]:=
        {(-16*a*Cos[u]*Cosh[v]*Sin[u]*Sinh[v]*
         (Cos[u]*Cosh[v]*Sin[t] + Cos[t]*Sin[u]*Sinh[v]))/
         (18 + Cos[4*u] - 12*Cos[2*u]*Cosh[2*v] + Cosh[4*v]) +
         (4*a*(3 - Cos[2*u]*Cosh[2*v])*
         (Cos[t]*Cos[u]*Cosh[v] - Sin[t]*Sin[u]*Sinh[v]))/
         (18 + Cos[4*u] - 12*Cos[2*u]*Cosh[2*v] + Cosh[4*v]),
         (a*(-Sin[t - 4*u] + 6*Cosh[2*v]*Sin[t - 2*u] - 
          6*Cosh[2*v]*Sin[t + 2*u] + Sin[t + 4*u] - 
          6*Sin[t - 2*u]*Sinh[2*v] -  6*Sin[t + 2*u]*Sinh[2*v] +
          2*Sin[t]*Sinh[4*v]))/(2*(-18 - Cos[4*u] +
          12*Cos[2*u]*Cosh[2*v] - Cosh[4*v])),
          a*(-Im[EllipticF[u + I*v,-1]]Cos[t] + 
          Re[EllipticF[u + I*v,-1]]*Sin[t])}

    lemniscatemin[a_][t_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    lemniscatemin[a][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    logminpolar[n_][t_][r_,theta_]:=
        analytictominimalpolar[logmincurve[n]][t][r,theta]

    logminpolar[n_][t_:0][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	logminpolar[n][t][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    menn[a_][u_,v_]:= {u,v,a*u^4 + u^2*v - v^2}

    menn[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    menn[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    mercatorellipsoid[a_,b_,c_][u_,v_]:=
        {a*Sech[v]*Cos[u],b*Sech[v]*Sin[u],c*Tanh[v]}

    mercatorellipsoid[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    mercatorellipsoid[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    mercatorsphere[a_][u_,v_]:=
        {a*Sech[v]*Cos[u],a*Sech[v]*Sin[u],a*Tanh[v]}

    mercatorsphere[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    mercatorsphere[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    mercatorsphere[a_,b_][u_,v_]:=
        {a*(Cos[u] + b*Cosh[v])/(Cosh[v] + b*Cos[u]),
         a*Sqrt[1 - b^2]*Sin[u]/(Cosh[v] + b*Cos[u]),
         a*Sqrt[1 - b^2]*Sinh[v]/(Cosh[v] + b*Cos[u])}

    mercatorsphere[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    mercatorsphere[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    moebcirc[a_][u_,v_]:={a*(-2*Cos[2*v]*Sin[u])/
     (-2 + Sqrt[2]*Cos[u]*Sin[v] + Sqrt[2]*Sin[u]*Sin[2*v]),
      a*(Sqrt[2]*(Cos[u]*Sin[v] - Sin[u]*Sin[2*v]))/
     (-2 + Sqrt[2]*Cos[u]*Sin[v] + Sqrt[2]*Sin[u]*Sin[2*v]),
       a*(-2*Cos[u]*Cos[v])/
     (-2 + Sqrt[2]*Cos[u]*Sin[v] + Sqrt[2]*Sin[u]*Sin[2*v])}

    moebcirc[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	moebcirc[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    moebcirc[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    moebcirc[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    moebiusstrip[a_][u_,v_]:= {a*(Cos[u] + v*Cos[u/2]*Cos[u]),
        a*(Sin[u] + v*Cos[u/2]*Sin[u]),a*v*Sin[u/2]}

    moebiusstrip[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    moebiusstrip[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    monkeysaddle[u_,v_]:= {u,v,u^3 - 3*u*v^2}
    
    monkeysaddle[umin_,umax_,vmin_,vmax_,opt___]:=
	ParametricPlot3D[Evaluate[
	monkeysaddle[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    monkeysaddle[n_][u_,v_]:={u,v,
        Simplify[ComplexExpand[Re[(u + I*v)^n],
           TargetFunctions->{Re,Im}]]}

    monkeysaddlepolar[n_][r_,theta_]:=
        {r*Cos[theta],r*Sin[theta],r^n*Cos[n*theta]}

    monkeysaddlepolar[n_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	monkeysaddlepolar[n][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    nodoidnds[vmin_,vmax_][a_,b_,c_,e_][u_,v_]:=
         Module[{tmpintr},
         tmpintr=
             intrinsic[delaunaykappa2[c,e],0,{0,0,-Pi/2},
             MaxSteps->2000,{vmin,vmax}][v];
         {Cos[u]*(a + b*tmpintr[[1]]),Sin[u]*(a + b*tmpintr[[1]]),
            b*tmpintr[[2]]}]
            
    nodoid[vmin_,vmax_][a_,b_,c_,e_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    nodoidnds[vmin,vmax][a,b,c,e][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
     parabolamin[a_][t_][u_,v_]:=
        {2*a*(u*Cos[t] + v*Sin[t]),
         a*((u^2 - v^2)*Cos[t] + 2*u*v*Sin[t]),
         (a*(-2*(u*v + ArcTan[u + Sqrt[1 + u^2],2*v])*Cos[t] + 
         (-2*v^2 + Log[(u + Sqrt[1 + u^2])^2 + 4*v^2])*Sin[t] + 
          2*Sqrt[1 + u^2]*(-(v*Cos[t]) + u*Sin[t])))/2}

    parabolamin[a_][t_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    parabolamin[a][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    paraboliccoor[d_,e_,f_][a_,b_,c_][w_][u_,v_]:=
         {d*Sqrt[(a - u)*(a - v)*(a - w)/(b - a)],
          e*Sqrt[(b - u)*(b - v)*(b - w)/(a - b)],
          f*(c - a - b + u + v + w)/2}

    paraboliccoor[d_,e_,f_][a_,b_,c_][w_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    paraboliccoor[d,e,f][a,b,c][w][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    paraboliccyclide[a_,k_][u_,v_]:=
        {u*(8*a^2 + k + v^2)/(16*a^2 + u^2 + v^2),
         v*(8*a^2 - k + u^2)/(16*a^2 + u^2 + v^2),
        (16*a^2*(k - u^2 + v^2) - k*(u^2 + v^2))/
        (8*a*(16*a^2 + u^2 + v^2))}

    paraboliccyclide[a_,k_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    paraboliccyclide[a,k][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    paraboloid[a_][u_,v_]:= {u,v,a*(u^2 + v^2)}

    paraboloid[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    paraboloid[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    paraboloid[a_,b_][u_,v_]:= {u,v,a*u^2 + b*v^2}

    paraboloid[a_,b_,c_][u_,v_]:=
    {u,v,a*u^2 + b*v^2 + c*u*v}

    paraboloidpolar[n_,a_,b_][r_,theta_]:=
    {a*r*Cos[theta],b*r*Sin[theta],b*r^n}

    paraboloidpolar[n_,a_,b_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	paraboloidpolar[n,a,b][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    pear[a_,b_][u_,v_]:= {b*Cos[u]*Cos[v]*(1 + Sin[v]),
                          b*Sin[u]*Cos[v]*(1 + Sin[v]),
                            -a*(1 + Sin[v])}
                            
    pear[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
	ParametricPlot3D[Evaluate[
	pear[a,b][u,v]],
	{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    pearbis[u_,v_]:={Cos[u]*(4 + 3.8*Cos[v]),
                     Sin[u]*(4 + 3.8*Cos[v]),
        ((Cos[v] - 3.5)*(1 + Sin[v])*Log[1 - Pi*v/10] +
          10*Sin[v])*(1 - 0.4*Cos[v/2]*(1 + 0.2*Sin[v^2]))}
          
    pearbis[umin_,umax_,vmin_,vmax_,opt___]:=
	ParametricPlot3D[Evaluate[
	pearbis[u,v]],
	{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    perturbedms[a_][u_,v_]:=
            {u,v,u^3 - 3*u*v^2 + a*(u^2 + v^2)}

    perturbedms[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    perturbedms[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    perturbedms[a_,b_,c_][u_,v_]:=
            {u,v,c*(u^3 - 3*u*v^2) + a*u^2 + b*v^2}

    perturbedms[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    perturbedms[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    perturbedmspolar[n_,a_][r_,theta_]:=
        {r*Cos[theta],r*Sin[theta],r^n*Cos[n*theta] + a*r^2}

    perturbedmspolar[n_,a_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	perturbedmspolar[n,a][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    pillow[a_][u_,v_]:= {a*Cos[u],a*Cos[v],a*Sin[u]*Sin[v]}

    pillow[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    pillow[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    plane[a1_,a2_,b1_,b2_,c1_,c2_][u_,v_]:=
        {a1*u + a2*v,b1*u + b2*v,c1*u + c2*v}

    planepolar[a_][r_,theta_]:= {a*r*Cos[theta],a*r*Sin[theta],0}

    planepolar[a_,b_][r_,theta_]:= {a*r*E^(b*theta)*Cos[theta],
                                    a*r*E^(b*theta)*Sin[theta],0}

    planepolar[a_,b_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	planepolar[a,b][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    planetocat[t_][u_,v_]:=
        {Cos[t]*u + Sin[u]*Cosh[v],Sin[t]*Cos[u]*Cosh[v],
         v + Cos[t]*Cos[u]*Sinh[v]}
            
    planetocat[t_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    planetocat[t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    plucker[u_,v_]:= {u,v,2*u*v/(u^2 + v^2)}
    
    plucker[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    plucker[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    plucker[n_][u_,v_]:= {u,v,Sin[n*ArcTan[u,v]]}
        
    plucker[m_,n_,a_][u_,v_]:= {a*u,a*v,
                    a*(u + v)^(m/2)*Sin[n*ArcTan[u,v]]}

    pluckerpolar[n_][r_,theta_]:=
                    {r*Cos[theta],r*Sin[theta],Sin[n*theta]}

    pluckerpolar[m_,n_,a_][r_,theta_]:=
        {a*r*Cos[theta],a*r*Sin[theta],a*r^m*Sin[n*theta]}

    pluckerpolar[m_,n_,a_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	pluckerpolar[m,n,a][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    pseudocrosscap[u_,v_]:=  
                 {(1 - u^2)*Sin[v],(1 - u^2)*Sin[2*v],u}

    pseudocrosscap[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    pseudocrosscap[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    pseudosphere[a_][u_,v_]:= {a*Cos[u]*Sin[v],
                               a*Sin[u]*Sin[v],
                               a*(Cos[v] + Log[Tan[v/2]])}

    pseudosphere[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    pseudosphere[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    pseudosphere[a_,b_][u_,v_]:= {a*Cos[u]*Sin[v],
                                  a*Sin[u]*Sin[v],
                                  b*(Cos[v] + Log[Tan[v/2]])}

    pseudosphere[a_,b_,c_][u_,v_]:= {a*Cos[u]*Sin[v],
                                     b*Sin[u]*Sin[v],
                                     c*(Cos[v] + Log[Tan[v/2]])}

    pseudospherebis[a_][u_,v_]:= {a*Cos[u]*Tanh[v],
                                  a*Sin[u]*Tanh[v],
                                  a*(Sech[v] + Log[Tanh[v/2]])}

    pseudospherebis[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    pseudospherebis[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    pseudospheretcheb[a_][u_,v_]:= {a*Cos[v]/Cosh[u],
                                    a*Sin[v]/Cosh[u],
                                    a*(u - Tanh[u])}

    pseudospheretcheb[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    pseudospheretcheb[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    pseudospheretcheb[a_,b_][u_,v_]:= {a*Cos[v]/Cosh[u],
                                       a*Sin[v]/Cosh[u],
                                       b*(u - Tanh[u])}

    pseudospheretcheb[a_,b_,c_][u_,v_]:= {a*Cos[v]/Cosh[u],
                                          b*Sin[v]/Cosh[u],
                                          c*(u - Tanh[u])}

    pseudospheretcheb[a_,b_,c_,d_][u_,v_]:= {a*Cos[v]/Cosh[u],
                                             b*Sin[v]/Cosh[u],
                                             c*u - d*Tanh[u]}

    rconoid[u_,v_]:= {v*Cos[u],v*Sin[u],2*Sin[u]}

    rconoid[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    rconoid[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    richmond[u_,v_]:=
        {(-3*u - u^5 + 2*u^3*v^2 + 3*u*v^4)/(6*(u^2 + v^2)),
         (-3*v - 3*u^4*v - 2*u^2*v^3 + v^5)/(6*(u^2 + v^2)),u}

    richmond[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    richmond[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    richmond[n_][t_][u_,v_]:= Simplify[ComplexExpand[
        analytictominimal[richmondmincurve[n]][t][u,v],
        TargetFunctions->{Re,Im}]]

    richmondpolar[n_][t_][r_,theta_]:=
        {-Cos[t + theta]/(2*r) - 
          r^(2*n + 1)*Cos[t - (2*n + 1)*theta]/(4*n + 2),
         -Sin[t + theta]/(2*r) + 
          r^(2*n + 1)*Sin[t - (2*n + 1)*theta]/(4*n + 2),
          r^n*Cos[t - n*theta]/n}

    richmondpolar[n_][t_:0][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	richmondpolar[n][t][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    roman[a_][u_,v_]:= {(a^2/2)*Sin[2*u]Cos[v]^2,
                        (a^2/2)*Sin[u]*Sin[2*v],
                        (a^2/2)*Cos[u]*Sin[2*v]}

    roman[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    roman[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    roman[n_,a_][u_,v_]:= {(a^2/2)*(Cos[u]*Cos[v]^2*Sin[u])^n,
                           (a^2/2)*(Cos[v]*Sin[u]*Sin[v])^n,
                           (a^2/2)*(Cos[u]*Cos[v]*Sin[v])^n}

    rosemin[n_,a_][t_][u_,v_]:=
    {Cos[n*u]*Cosh[n*v]*(a*Cos[t]*Cos[u]*Cosh[v] - 
       a*Sin[t]*Sin[u]*Sinh[v]) + 
         Sin[n*u]*(-a*Cos[u]*Cosh[v]*Sin[t] - 
         a*Cos[t]*Sin[u]*Sinh[v])*Sinh[n*v],
       -Cos[u]*Sinh[v]*(-a*Cos[n*u]*Cosh[n*v]*Sin[t] - 
         a*Cos[t]*Sin[n*u]*Sinh[n*v]) + 
         Cosh[v]*Sin[u]*(a*Cos[t]*Cos[n*u]*Cosh[n*v] - 
         a*Sin[t]*Sin[n*u]*Sinh[n*v]),
       -a*Cos[t]*Im[EllipticE[n*(u + I*v),1 - n^2]]/n + 
        a*Re[EllipticE[n*(u + I*v),1 - n^2]]*Sin[t]/n}

    rosemin[n_,a_][t_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    rosemin[n,a][t][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    scherk[a_][u_,v_]:= {u,v,(1/a)*Log[Cos[a*v]/Cos[a*u]]}

    scherk[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    scherk[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    scherk[a_,theta_][u_,v_]:=
        {Sec[2*theta]*(Cos[theta]*u + Sin[theta]*v),
         Sec[2*theta]*(Sin[theta]*u + Cos[theta]*v),
         (1/a)*Log[
         Cos[a*Sec[2*theta]*(Sin[theta]*u + Cos[theta]*v)]/
         Cos[a*Sec[2*theta]*(Cos[theta]*u + Sin[theta]*v)]]}
   
    scherk5[u_,v_]:= {ArcSinh[u],ArcSinh[v],ArcSin[u*v]}

    scherk5[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    scherk5[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    scherk5[a_,b_,c_][u_,v_]:= {a*ArcSinh[u],
                                b*ArcSinh[v],
                                c*ArcSin[u*v]}

    shoe[a_,b_][u_,v_]:= {u,v,a*u^3 + b*v^2} 

    shoe[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
	ParametricPlot3D[Evaluate[
	shoe[a,b][u,v]],
	{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
    
    shoeasym[a_,b_][p_,q_]:=
    {(-3*b/(4*a))^(1/3)*(p - q)^(2/3),
                      -p - q,(b/4)*(p^2 + 14*p*q + q^2)}

    shoeasym[a_,b_][pmin_,pmax_,qmin_,qmax_,opt___]:=
    ParametricPlot3D[Evaluate[Re[
	 shoeasym[a,b][p,q]]],{p,pmin,pmax},{q,qmin,qmax},
	opt,Axes->None];

    sievert[a_][u_,v_]:=
       {(2/(a + 1 - a*Sin[v]^2*Cos[u]^2))*
        (Sqrt[(a + 1)*(1 + a*Sin[u]^2)]*Sin[v]/Sqrt[a])*
         Cos[-u/Sqrt[a + 1] + ArcTan[Sqrt[a + 1]*Tan[u]]],
         (2/(a + 1 - a*Sin[v]^2*Cos[u]^2))*
         (Sqrt[(a + 1)*(1 + a*Sin[u]^2)]*Sin[v]/Sqrt[a])*
         Sin[-u/Sqrt[a + 1] + ArcTan[Sqrt[a + 1]*Tan[u]]],
         Log[Tan[v/2]]/Sqrt[a] + 2*(a + 1)*Cos[v]/
                   ((a + 1 - a*Sin[v]^2*Cos[u]^2)*Sqrt[a])}

    sievert[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    sievert[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    sinovaloid[n_,a_][u_,v_]:= {a*Cos[u]*Cos[v],
                                a*Sin[u]*Cos[v],
                                a*Nest[Sin,v,n]}

    sinovaloid[n_,a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    sinovaloid[n,a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    sinsurface[a_][u_,v_]:=
        {a*Sin[u],a*Sin[v],a*Sin[u + v]}

    sinsurface[a_][umin_,umax_,vmin_,vmax_,opt___]:=
	ParametricPlot3D[Evaluate[
	sinsurface[a][u,v]],
	{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
    
    snail[u_,v_]:=
        {u*Cos[v]*Sin[u],u*Cos[u]*Cos[v],-u*Sin[v]}

    snail[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    snail[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    sphere[a_][u_,v_]:=
        {a*Cos[v]*Cos[u],a*Cos[v]*Sin[u],a*Sin[v]}
        
    sphere[a_][umin_,umax_,vmin_,vmax_,opt___]:=
	ParametricPlot3D[Evaluate[
	sphere[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    stereographicellipsoid[a_,b_,c_][u_,v_]:=
        {2*a*u/(u^2 + v^2 + 1),2*b*v/(u^2 + v^2 + 1),
         c*(u^2 + v^2 - 1)/(u^2 + v^2 + 1)}

    stereographicellipsoid[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
	ParametricPlot3D[Evaluate[
	stereographicellipsoid[a,b,c][u,v]],
	{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
    
    stereographicsphere[a_][u_,v_]:=
        {2*a*u/(u^2 + v^2 + 1),2*a*v/(u^2 + v^2 + 1),
         a*(u^2 + v^2 - 1)/(u^2 + v^2 + 1)}

    stereographicsphere[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    stereographicsphere[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    stereographicspherepolar[a_][r_,theta_]:=
        {2*a*r*Cos[theta]/(1 + r^2),
         2*a*r*Sin[theta]/(1 + r^2),
         a*(-1 + r^2)/(1 + r^2)}

    stereographicspherepolar[a_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	stereographicspherepolar[a][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    stereographicspherepolar[a_,b_][r_,theta_]:=
        {(2*a*E^(b*theta)*r*Cos[theta])/(1 + E^(2*b*theta)*r^2),
         (2*a*E^(b*theta)*r*Sin[theta])/(1 + E^(2*b*theta)*r^2),
         (a*(-1 + E^(2*b*theta)*r^2))/(1 + E^(2*b*theta)*r^2)}

    stereographicspherepolar[a_,b_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	stereographicspherepolar[a,b][r,theta]],
	{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];
    
    swallowtail[u_,v_]:={3*v^4 + u*v^2,-4*v^3 - 2*u*v,u}

    swallowtail[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    swallowtail[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    tetrahedral[m_,n_][A_,B_,C_][a_,b_,c_][u_,v_]:=
              {A*(u - a)^m*(v - a)^n,
               B*(u - b)^m*(v - b)^n,
               C*(u - c)^m*(v - c)^n}

    thomsen[a_,b_][u_,v_]:=
        {b*u/a + Sqrt[1+b^2]*Sinh[a*u]*Cos[a*v]/a^2,
        Sqrt[1 + b^2]*v/a + b*Cosh[a*u]*Sin[a*v]/a^2,
        Sinh[a*u]*Sin[a*v]/a^2}

    thomsen[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    thomsen[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    torus[a_,b_][u_,v_]:= {(a + b*Cos[v])*Cos[u],
                           (a + b*Cos[v])*Sin[u],
                            b*Sin[v]}

    torus[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    torus[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    torus[a_,b_,c_][u_,v_]:= {(a + b*Cos[v])*Cos[u],
                              (a + b*Cos[v])*Sin[u],
                               c*Sin[v]}

    torus[a_,b_,c_,d_][u_,v_]:= 
        {(a + b*Cos[v])*Cos[u],(d + b*Cos[v])*Sin[u],c*Sin[v]}

    torus[a_,b_,c_,d_,e_][u_,v_]:=
        {(a  +b*Cos[v])*Cos[u],(d + e*Cos[v])*Sin[u],c*Sin[v]}

    torusasym[pm_,a_][p_,q_]:=
        {(2*a*Cos[(p + q)/2]/(1 + Cosh[(p - q)/Sqrt[8]]^2)),
         (2*a*Sin[(p + q)/2]/(1 + Cosh[(p - q)/Sqrt[8]]^2)),
         pm*(2*a*Cosh[(p - q)/Sqrt[8]]/(1 + Cosh[(p - q)/Sqrt[8]]^2))}

    torusasym[pm_,a_][pmin_,pmax_,qmin_,qmax_,opt___]:=
    ParametricPlot3D[Evaluate[
	 torusasym[pm,a][p,q]],{p,pmin,pmax},{q,qmin,qmax},
	opt,Axes->None];

    torusbis[pm_,a_,b_][r_,theta_]:={r*Cos[theta],r*Sin[theta],
                                     pm*Sqrt[b^2 - (r - a)^2]}

    torusbis[pm_,a_,b_][rmin_,rmax_,thetamin_,thetamax_,opt___]:=
    ParametricPlot3D[Evaluate[
	torusbis[pm,a,b][r,theta]],{r,rmin,rmax},{theta,thetamin,thetamax},
	opt,Axes->None];

    twiflat[pm_,a_,c_][sl_][u_,v_]:= {v*Cos[u],v*Sin[u],
        sl*u + ArcSin[c/(v*a)]*c + Sqrt[v^2*a^2 - c^2]*pm}

    twiflat[pm_,a_,c_][sl_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    twiflat[pm,a,c][sl][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    twisphere[a_,b_][u_,v_]:= {a*Cos[u]*Cos[v],
                               a*Sin[u]*Cos[v],
                               a*Sin[v] + b*u}

    twisphere[a_,b_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    twisphere[a,b][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];

    veronese[a_][u_,v_]:= {u,v,u^2,u v,v^2}

    wallis[a_,b_,c_][u_,v_]:= {v*Cos[u],v*Sin[u],
                               c*Sqrt[a^2 - b^2*Cos[u]^2]}
 
    wallis[a_,b_,c_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    wallis[a,b,c][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    whitneyumbrella[u_,v_]:= {u*v,u,v^2}

    whitneyumbrella[umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    whitneyumbrella[u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
    wrinkles[a_][u_,v_]:= {u,v,a*(u/v + v/u)}

    wrinkles[a_][umin_,umax_,vmin_,vmax_,opt___]:=
    ParametricPlot3D[Evaluate[
    wrinkles[a][u,v]],{u,umin,umax},{v,vmin,vmax},
	opt,Axes->None];
     
(* Implicitly Defined Surfaces *)

    cossurfaceimplicit[a_][x_,y_,z_]:=
               (z - a*x*y)^2 - a^2(1 - x^2)(1 - y^2)

    crosscapimplicit[a_,b_][x_,y_,z_]:=
        (a*x^2 + b*y^2)(x^2 + y^2 + z^2) - 2*z*(x^2 + y^2)

    ellipsoidimplicit[a_,b_,c_][x_,y_,z_]:=
               (x/a)^2 + (y/b)^2 + (z/c)^2 - 1

    ellipticparaboloidimplicit[a_,b_,c_][x_,y_,z_]:=
               (x/a)^2 + (y/b)^2 - c*z
    
    equihom1implicit[x_,y_,z_]:= x*y*z - 1

    equihom2implicit[x_,y_,z_]:= (x^2 + y^2)*z - 1

    equihom3implicit[x_,y_,z_]:= x^2*(z - y^2)^3 - 1

    equihom4implicit[x_,y_,z_]:= x^2*(z - y^2)^3 + 1

    equihom5implicit[x_,y_,z_]:= z - x*y - x^3/3

    equihom6implicit[x_,y_,z_]:= z - x*y + Log[x]

    goursatimplicit[a_,b_,c_][x_,y_,z_]:= x^4 + y^4 + z^4 +
                   a*(x^2 + y^2 + z^2)^2 + b*(x^2 + y^2 + z^2) + c

    hyperboloidimplicit[a_,b_,c_][x_,y_,z_]:=
               (x/a)^2 + (y/b)^2 - (z/c)^2 - 1

    hy2sheetimplicit[a_,b_,c_][x_,y_,z_]:=
               (x/a)^2 - (y/b)^2 - (z/c)^2 - 1
    
    hyperbolicparaboloidimplicit[a_,b_,c_][x_,y_,z_]:=
               (x/a)^2 - (y/b)^2 - c*z

    kazoolaimplicit[a_,b_,c_,d_][x_,y_,z_]:=
        d + c*x^2y^2z^2 - a*(x^2 + y^2 + z^2) - b*(x^4 + y^4 + z^4)

    kummerimplicit[x_,y_,z_]:= x^4 + y^4 + z^4 -
             (y^2*z^2 + z^2*x^2 + x^2*y^2) - (x^2 + y^2 + z^2) + 1

    scherkimplicit[a_][x_,y_,z_]:=
               E^(a*z)*Cos[a*x] -  Cos[a*y]

    scherkimplicit[a_,phi_][x_,y_,z_]:=
               Tan[x/(a*Cos[phi])]*Tan[y/(a*Sin[phi])] - Tanh[z/a]

    scherk5implicit[x_,y_,z_]:=
               Sinh[x]*Sinh[y] - Sin[z]

    sinsurfaceimplicit[a_][x_,y_,z_]:=
               4*x^2*y^2*(a^2 - z^2) - a^2*(x^2 + y^2 - z^2)^2

    sphereimplicit[a_][x_,y_,z_]:=
               x^2 + y^2 + z^2 - a^2

    sphereimplicit[n_,a_][x_,y_,z_]:=
               x^(2*n) + y^(2*n) + z^(2*n) - a^(2*n)

    torusimplicit[a_,b_][x_,y_,z_]:=
              z^2 + (Sqrt[x^2 + y^2] - a)^2 - b^2

    twocuspsimplicit[x_,y_,z_]:= (z - 1)^2*(x^2 - z^2) -
           (x^2 - z)^2 - y^4 - y^2*(2*x^2 + z^2 + 2*z - 1)

(* Drum Plots *)

    circularmembrane[m_,n_,t_][r_,theta_]:=
        Flatten[{r*Cos[theta],r*Sin[theta],
            BesselJ[m,r*BesselJZeros[m,{n,n}]]*
            Cos[m*theta]*
            Cos[t*BesselJZeros[m,{n,n}]]}]

(*
    circularmembrane[m_,n_,t_][r_,theta_]:=
        {r*Cos[theta],r*Sin[theta],
            BesselJ[m,r*If[ListQ[BesselJZeros[m,{n,n}]]==True,
                  First[BesselJZeros[m,{n,n}]],BesselJZeros[m,{n,n}]]]*
            Cos[m*theta]*
            Cos[t*If[ListQ[BesselJZeros[m,{n,n}]==True],
                  First[BesselJZeros[m,{n,n}]],BesselJZeros[m,{n,n}]]]}
*)

      circularmembraneplot[m_,n_,t_,opts___]:=
        Module[{zzz},
             zzz=circularmembrane[m,n,t][r,theta];
             ParametricPlot3D[Evaluate[zzz],
                 {r,0,1},{theta,0,2*Pi},opts]]

      circularmembranemovie[m_,n_,{t_,tmin_,tmax_,numf_},
          opts___]:=
          Do[Module[{zzz},
             zzz=circularmembrane[m,n,t][r,theta];
             ParametricPlot3D[Evaluate[zzz],
                 {r,0,1},{theta,0,2*Pi},opts]],
                 {t,tmin,tmax,numf - 1}]

      squaremembrane[m_,n_,t_][u_,v_]:=
             {u,v,Sin[m*u]*Sin[n*v]*Cos[Sqrt[m^2 + n^2]*t]}

      squaremembraneplot[m_,n_,t_,opts___]:=
             ParametricPlot3D[Evaluate[
                squaremembrane[m,n,t][u,v]],
                 {u,0,Pi},{v,0,Pi},opts]

      squaremembranemovie[m_,n_,{t_,tmin_,tmax_,numf_},
          opts___]:=
             Do[ParametricPlot3D[Evaluate[
                squaremembrane[m,n,t][u,v]],
                 {u,0,Pi},{v,0,Pi},opts],
                 {t,tmin,tmax,numf - 1}]

(* Minimal Curves *)

    astroidmincurve[a_][z_]:= {a*Cos[z]^3,a*Sin[z]^3,
                               -3*I/4*a*Cos[2*z]}

    bourmincurve[m_][z_]:=
        {z^(m - 1)/(m - 1) - z^(m + 1)/(m + 1),
         I*(z^(m - 1)/(m - 1) + z^(m + 1)/(m + 1)),2*z^m/m}

    bourmincurve[m_,n_][z_]:=
        {z^(m - 1)/(m - 1) - z^(m + 2*n - 1)/(m + 2*n - 1),
         I*(z^(m - 1)/(m - 1) + z^(m + 2*n - 1)/(m + 2*n - 1)),
         2*z^(m + n - 1)/(m + n - 1)}

    cardioidmincurve[a_][z_]:= {2*a*Cos[z]*(1 + Cos[z]),
                  2*a*(1 + Cos[z])*Sin[z],8*I*a*Sin[z/2]}

    catalanmincurve[a_][z_]:= {a*z - a*Sin[z],a - a*Cos[z],
                                    -4*I*a*Cos[z/2] + 4*I*a}

    catenoidmincurve[a_][z_]:= {a*Cosh[z/a],z,-I*a*Sinh[z/a]}

(*
different: harmonictoanalytic[catenoid[a]][z]
*)

    cayleysexticmincurve[a_][z_]:=
         {4*a*Cos[z/2]^4*(-1 + 2*Cos[z]),
          4*a*Cos[z/2]^3*Sin[3*z/2],3*I*a*(z + Sin[z])}

    circleinvolutemincurve[n_,a_][z_]:=
        Module[{tmp1,tmp2,tmp3},
            tmp1=E^(I*tt)Normal[Series[E^(-I*tt),{tt,0,n}]];
            tmp2=a*ComplexExpand[{Re[tmp1],Im[tmp1]}];
            tmp3=tmp2 /. tt->z;
            Append[tmp3,I*a*z^(n + 1)/(n + 1)!]]

    cissoidmincurve[a_][z_]:= {2*a*z^2/(1 + z^2),
         2*a*z^3/(1 + z^2),2*I*a*(Sqrt[4 + z^2] - 
          Sqrt[3]*ArcTanh[Sqrt[4 + z^2]/Sqrt[3]])}

    clothoidmincurve[a_][z_]:=
            {a*Sqrt[Pi]*FresnelS[z/Sqrt[Pi]],
             a*Sqrt[Pi]*FresnelC[z/Sqrt[Pi]],I*a*z}

    cnchyprofilemincurve[pm_,a_,b_][z_]:=
        {pm*b*Cosh[z/a],-I*a*EllipticE[I*z/a,-b^2/a^2],I*z}

    costamincurve[z_]:=
        {(1/2)*(-WeierstrassZeta[z,{189.07272,0}] + Pi*(z - I) +
        (1 + I)*0.35889 +
        0.228473*(WeierstrassZeta[z - 1/2,{189.07272,0}] -
             WeierstrassZeta[z - I/2,{189.07272,0}])),
        (I/2)*(-WeierstrassZeta[z,{189.07272,0}] + Pi*(-z + 1) -
        (1 + I)*0.35889 -
        0.228473*(WeierstrassZeta[z - 1/2,{189.07272,0}] -
             WeierstrassZeta[z - I/2,{189.07272,0}])),
        Sqrt[2*Pi]/4*(Log[(WeierstrassP[z,{189.07272,0}] - 6.87519)/
             (WeierstrassP[z,{189.07272,0}] + 6.87519)] - Pi*I)}

    cpcprofilemincurve[pm_,a_,b_][z_]:= {b*pm*Cos[z/a],
                          a*EllipticE[z/a,b^2/a^2],I*z}

    deltoidinvolutemincurve[a_][z_]:=
        {(a/3)*(8*Cos[z/2] + 2*Cos[z] - Cos[2*z]),
         (a/3)*(-8*Sin[z/2] + 2*Sin[z] + Sin[2*z]),
          (4*a*I/9)*(3*z - 2*Sin[3*z/2])}

    deltoidmincurve[a_][z_]:= {-a + 2*a*Cos[z]*(1 + Cos[z]),
             2*a*(1 - Cos[z])*Sin[z],-8*(I/3)*a*Cos[3*z/2]}
    
    ellipsemincurve[a_,b_][z_]:= {a*Cos[z],b*Sin[z],
                                  I*b*EllipticE[z,-a^2/b^2]}

    ennepercatenoidmincurve[n_][z_]:=
        {1/z + z - 2*z^n/n - z^(2*n + 1)/(2*n + 1),
         I*(2*z^(n + 1) + 2*n^2*(z^2 - 1) + n*(-1 + z^2 +
         4*z^(n + 1) + z^(2*(n + 1))))/(n*(2*n + 1)*z),
         2*z^(n + 1)/(n + 1) + 2*Log[z]}

    ennepercatenoidmincurve[n_,a_,b_][z_]:=
        {a^2/z + z - 2*a*b*z^n/n - b^2*z^(2*n + 1)/(2*n + 1),
         I*(-a^2*n*(2*n + 1) + 2*a*b*(2*n + 1)*z^(n + 1) + 
         n*z^2*(1 + 2*n + b^2*z^(2*n)))/(n*(2*n + 1)*z),
         2*b*z^(n + 1)/(n + 1) + 2*a*Log[z]}    

    ennepermincurve[n_][z_]:= {z - z^(2*n + 1)/(2*n + 1),
        I*(z + z^(2*n + 1)/(2*n + 1)),2*z^(n + 1)/(n + 1)}

    epicycloidmincurve[a_,b_][z_]:=
              {(a + b)*Cos[z] - b*Cos[(a + b)*z/b],
               (a + b)*Sin[z] - b*Sin[(a + b)*z/b],
               -4*I*b*(a + b)*Cos[a*z/(2*b)]/a}

    hennebergmincurve[z_]:= {-8*Sinh[z]^3/3,
          8*I/3 - 2*I*Cosh[z] - 2*I/3*Cosh[3*z],2*Cosh[2*z]}

    hyperbolamincurve[a_,b_][z_]:= {a*Cos[z],b*Sin[z],
                                    b*EllipticE[I*z + a^2/b^2]}

    hypocycloidinvolutemincurve[a_,b_][z_]:=
       {(a^2*Cos[z] - 3*a*b*Cos[z] + 2*b^2*Cos[z] - 
         a*b*Cos[z - a*z/b] + 2*b^2*Cos[z - a*z/b] + 
         4*a*b*Cos[z - a*z/(2*b)] - 4*b^2*Cos[z - a*z/(2*b)])/a,
        (a^2*Sin[z] - 3*a*b*Sin[z] + 2*b^2*Sin[z] - 
         a*b*Sin[z - a*z/b] + 2*b^2*Sin[z - a*z/b] + 
         4*a*b*Sin[z - a*z/(2*b)] -  4*b^2*Sin[z - a*z/(2*b)])/a,
        2*I*(a^2 - 3*a*b + 2*b^2)*(a*z - 2*b*Sin[a*z/(2*b)])/a^2}

    hypocycloidmincurve[a_,b_][z_]:=
           {(a - b)*Cos[z] + b*Cos[(a - b)*z/b],
            (a - b)*Sin[z] - b*Sin[(a - b)*z/b],
            -4*I*(a - b)*b*Cos[a*z/(2*b)]/a}

    jorgemeeksmincurve[n_,a_][z_]:=
        {Integrate[jorgemeeksmincurveprime[n,a][zz][[1]],{zz,0,z}],
         Integrate[jorgemeeksmincurveprime[n,a][zz][[2]],{zz,0,z}],
         a*(1 + z^n)/(2*n*(1 - z^n))}

    jorgemeeksmincurve[2,a_][z_]:=
    {-I/4*a*Pi  + (a/4)*Log[(z + 1)/(z - 1)],
     -I*a*z/(2*(z^2 - 1)),a*(1 + z^2)/(4*(1 - z^2))}

    jorgemeeksmincurve[3,a_][z_]:=
    {(a/18)*(4*I*Pi + 3*z/(1 + z + z^2) - 4*Log[z - 1] + 2*Log[1 + z + z^2]),
      -(I*a/54)*(2*Sqrt[3]*Pi + 3*(3*z*(z + 1) - 
       4*Sqrt[3]*(z^3 - 1)*ArcTan[(2*z + 1)/Sqrt[3]])/(z^3 - 1)),
     a*(z^3 + 1)/(6*(1 - z^3))}

    jorgemeeksmincurve[4,a_][z_]:=
    {3*I/16*a*Pi + (a*(2*z - 3*(1 + z^2)*Log[-1 + z] + 
        3*(1 + z^2)*Log[1 + z]))/(16*(1 + z^2)),
     -I*a*z/(8*(z^2 - 1)) + 3*I*a*ArcTan[z]/8,
      a*(1 + z^4)/(8*(1 - z^4))}

    jorgemeeksmincurveprime[n_,a_][z_]:=
        {(a*(1 - z^(2*n - 2)))/(2*(z^n - 1)^2),
         (a*I*(1 + z^(2*n - 2)))/(2*(z^n - 1)^2),
          a*z^(n - 1)/(z^n - 1)^2}

    lemniscatemincurve[a_][z_]:=
        {(a*Cos[z])/(1 + Sin[z]^2),
         (a*Cos[z]*Sin[z])/(1 + Sin[z]^2),I*a*EllipticF[z,- 1]}

    logmincurve[n_][z_]:=
        Module[{ii,nn,zz},
            ii[nn_][zz_]:=Nest[Integrate[#,zz]&,1/zz,nn];
            wei[#&,ii[n]][z]]

    logspiralmincurve[a_,b_][z_]:= {a*E^(b*z)*Cos[z],
          a*E^(b*z)*Sin[z],(I*a*Sqrt[1 + b^2]*E^(b*z))/b}

    nephroidmincurve[a_][z_]:= {a*(3*Cos[z] - Cos[3*z]),
                   a*(3*Sin[z] - Sin[3*z]),-6*I*a*Cos[z]}

    nielsenspiralmincurve[a_][z_]:= 
        {a*CosIntegral[z],a*SinIntegral[z],I*a*Log[z]}

    parabolamincurve[a_][z_]:= {2*a*z,a*z^2,
        I*a*(z*Sqrt[1 + z^2] + Log[z + Sqrt[1 + z^2]])}

    richmondmincurve[n_][z_]:=
            {-1/(2*z) - z^(2*n + 1)/(4*n + 2),
             -I/(2*z) + I*z^(2*n + 1)/(4*n + 2),z^n/n}

    rosemincurve[n_,a_][z_]:= {a*Cos[z]*Cos[n*z],a*Cos[n*z]*Sin[z],
                               I*a*EllipticE[n*z,1 - n^2]/n}

    scherktowermincurveprime[n_][z_]:=
        {(1 - z^(2*n - 2))/(2*(1 + z^(2*n))),
         (I*(1 + z^(2*n - 2)))/(2*(1 + z^(2*n))),
          z^(n - 1)/(1 + z^(2*n))}

    scmincurve[z_]:= {z^2,z^3,
         I*(4/(27*z) + z/3)*Sqrt[4*z^2 + 9*z^4]}

    tractrixmincurve[a_][z_]:= {a*Sin[z],
                                a*(Cos[z] + Log[Tan[z/2]]),
                        I*a*(Log[Cos[z/2]] + Log[Sin[z/2]])}

(* Surface Metrics *)

    ee[geopolar[f_]][u_,v_]:= 1

    ff[geopolar[f_]][u_,v_]:= 0

    gg[geopolar[f_]][u_,v_]:= f[u,v]^2

    ee[liouville[f_,g_]][u_,v_]:= f[u] + g[v]

    ff[liouville[f_,g_]][u_,v_]:= 0

    gg[liouville[f_,g_]][u_,v_]:= f[u] + g[v]

    ee[poincare[lambda_]][u_,v_]:= 1/(lambda*v)^2

    ff[poincare[lambda_]][u_,v_]:= 0

    gg[poincare[lambda_]][u_,v_]:= 1/(lambda*v)^2

    ee[poincare[m_,n_,lambda_]][u_,v_]:= 1/(lambda*v)^m

    ff[poincare[m_,n_,lambda_]][u_,v_]:= 0

    gg[poincare[m_,n_,lambda_]][u_,v_]:= 1/(lambda*v)^n

    ee[poindisk[lambda_]][u_,v_]:=
                  1/(1 - lambda^2*(u^2 + v^2)/4)^2

    ff[poindisk[lambda_]][u_,v_]:= 0

    gg[poindisk[lambda_]][u_,v_]:=
                  1/(1 - lambda^2*(u^2 + v^2)/4)^2

    ee[schwarzschildplane[m_]][u_,v_]:=  -(1 - 2*m/v)

    ff[schwarzschildplane[m_]][u_,v_]:= 0

    gg[schwarzschildplane[m_]][u_,v_]:= 1/(1 - 2*m/v)

    ee[schwarzschildplane[m_,f_]][u_,v_]:= -(1 - 2*m/v + f[v])

    ff[schwarzschildplane[m_,f_]][u_,v_]:= 0

    gg[schwarzschildplane[m_,f_]][u_,v_]:= 1/(1 - 2*m/v + f[v])

End[]

EndPackage[]

Index of Surfaces

Computations

Area

Integrate[infarea[s][u, v],  {u,  (Subscript[u, i ]), 
        (Subscript[u, f])},  {v,  (Subscript[v, i]), (Subscript[v, f])}]

NIntegrate[infarea[s][u, v],  {u,  (Subscript[u, i ]), 
        (Subscript[u, f])},  {v,  (Subscript[v, i]), (Subscript[v, f])}]

Gassian Curvature

gaussiancurvature[s][u, v]

Plot3D[gaussiancurvature[s][u, v], {v, (Subscript[v, i]), (Subscript[v, f])},  {u, (Subscript[u, i]), 
        (Subscript[u, f])}]

Mean Curvature

meancurvature[s][u, v]

Plot3D[meancurvature[s][u, v], {v, (Subscript[v, i]), (Subscript[v, f])},  {u, (Subscript[u, i]), 
        (Subscript[u, f])}]

Plots

Basic Plot and Options

ParametricPlot3D[s[u, v],  {u, (Subscript[u, i]), (Subscript[u, f])},  {v, (Subscript[v, i]), (Subscript[v, f])}, 
      Axes -> False]

ParametricPlot3D[
      Append[s[u, v], Red],  {u, (Subscript[u, i]), (Subscript[u, f])},  {v, (Subscript[v, i]), (Subscript[v, f])}, 
      Axes -> False, Lighting -> False]

Options[ParametricPlot3D]

Plot Colored by Gaussian Curvature

colorfunction[u_, v_] := 
        If[gaussiancurvature[s][u, v] > 0, 
          0.5 + gaussiancurvature[s][u, v]/(2  maxgauss), 
          0.5 - gaussiancurvature[s][u, 
                v]/(2  mingauss)];
maxgauss = 
        Max[Table[gaussiancurvature[s][u, v], {u, (Subscript[u, i]), 
              (Subscript[u, f]), ((Subscript[u, f]) - (Subscript[u, i]))/20}, {v, (Subscript[v, i]), 
              (Subscript[v, f]), ((Subscript[v, f]) - (Subscript[v, i]))/20}]];
mingauss = 
        Min[Table[gaussiancurvature[s][u, v], {u, (Subscript[u, i]), 
              (Subscript[u, f]), ((Subscript[u, f]) - (Subscript[u, i]))/20}, {v, (Subscript[v, i]), 
              (Subscript[v, f]), ((Subscript[v, f]) - (Subscript[v, i]))/20}]];
ParametricPlot3D[
        Append[s[u, v], Hue[1/3 + 2  colorfunction[u, v]/3]] // 
          Evaluate,  {u, (Subscript[u, i]), (Subscript[u, f])},  {v, (Subscript[v, i]), (Subscript[v, f])}, 
        AspectRatio -> Automatic, Boxed -> False, 
        Axes -> None, Lighting -> False, 
        PlotPoints -> {30, 30}, PlotRange -> All];

Plot Colored by Mean Curvature

colorfunction[u_, v_] := 
        If[meancurvature[s][u, v] > 0, 
          0.5 + meancurvature[s][u, v]/(2  maxmean), 
          0.5 - meancurvature[s][u, v]/(2  minmean)];
maxmean = 
        Max[Table[meancurvature[s][u, v], {u, (Subscript[u, i]), 
              (Subscript[u, f]), ((Subscript[u, f]) - (Subscript[u, i]))/20}, {v, (Subscript[v, i]), 
              (Subscript[v, f]), ((Subscript[v, f]) - (Subscript[v, i]))/20}]];
minmean = 
        Min[Table[meancurvature[s][u, v], {u, (Subscript[u, i]), 
              (Subscript[u, f]), ((Subscript[u, f]) - (Subscript[u, i]))/20}, {v, (Subscript[v, i]), 
              (Subscript[v, f]), ((Subscript[v, f]) - (Subscript[v, i]))/20}]];
ParametricPlot3D[
        Append[s[u, v], Hue[1/3 + 2  colorfunction[u, v]/3]] // 
          Evaluate,  {u, (Subscript[u, i]), (Subscript[u, f])},  {v, (Subscript[v, i]), (Subscript[v, f])}, 
        AspectRatio -> Automatic, Boxed -> False, 
        Axes -> None, Lighting -> False, 
        PlotPoints -> {20, 40}, PlotRange -> All];