(*Sat Aug 29 18:57:11 EDT 1998
*)

(*:Title: CURVES.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 parametrically defined plane curves, implicitly defined plane curves, polar parametrizations and space curves. Short plotting commands are also given.
*)

(*:Keywords:
agnesi, agnesiimplicit, airyspiral, alysoid, alysoidnds, alysoidprime, archimedesspiral, archimedesspiralpolar, astroid, astroidimplicit, ast3d, baseballseam, bellcurve, besselcurve, besseln, bernstein, beziercurve, bicorn, bicylinder, biquadratic, bow, bowpolar, bowtie, bulletnose, bulletnoseimplicit, cardioid, cardioidimplicit, cardioidpolar, cardioidunitspeed, cartesianimplicit, cartesianoval, cassini, cassiniimplicit, cassinipolar, catenary, catenaryunitspeed, cayleysextic, circle, circleimplicit, circleinvolute, circleinvoluteunitspeed, circleunitspeed, cissoid, cissoidimplicit, cissoidoblique, clelia, clothoid, clothoidnds, clothoidprime, cnccoprofile, cnchyprofile, cochleoid, cochleoidpolar, conchoid, conchoidimplicit,  conchoidpolar, conicalhelix, coshspiral, cosn, cothspiral, cpcprofile, crosscurve, crosscurveimplicit, cubicalellipse, cubicparabola, cubicparabolaimplicit, curves3d, cycloid, cycloidunitspeed, delaunay, delaunaykappa2, delaunaynds, deltoid, deltoidbis, deltoidimplicit, deltoidinvolute, deltoidunitspeed, devilimplicit, diamond, eight, eightknot, elasticainflect, elasticanoninflect, ellipse, ellipsebis, ellipseimplicit, ellipseinc, epicycloid, epicycloidunitspeed, epispiral, epispiralpolar, epitrochoid, epitrochoid3d, fermatspiral, fermatspiralpolar, folium, foliumimplicit, foliumpolar, freeth, genepicycloid, genfolium, genparabola, genhelix, hankelspiral, helix, helixunitspeed, hippias, hippopede, hippopedeimplicit, horopter, hyperbola, hyperbolabis, hyperbolaimplicit, hyperbolicspiral, hyperbolicspiralpolar, hypocycloid, hypocycloidinvolute, hypocycloidunitspeed, hypotrochoid, kampyle, kampylepolar, kappacurve, kappacurveimplicit, keplerimplicit, keplerorbit, keplerorbitpolar, kochsnowflake, lehr, lehrkappa2, lehrnds, lemniscate, lemniscatebis, lemniscateimplicit, lemniscatepolar, lemniscateunitspeed, limacon, limaconimplicit, limaconpolar, line, linearpursuit, line3d, lissajous, lissajous3d, lituus, lituuspolar, logistic, logspiral, logspiralpolar, logspiralunitspeed, nephroid, nephroidimplicit, nephroidpolar, nephroidunitspeed, newtonphillipsimplicit, ngon, nielsenspiral, oneparametersubgroup, pacman, parabola, parabolaimplicit, parabolicspiral, perseusimplicit, piriform, piriformimplicit, predatorpreyimplicit, pseudocatenary, pseudocatenaryprime, pseudospheregeodesic, pseudosphericalloxodrome, rectangle, regularpolygon, reuleauxpolygon, riemannfractal, rose, rosepolar, sc, scarab, scimplicit, seiffertspiral, semicubic, semicubicimplicit, serpentine, serpentineimplicit, sinhspiral, sinn, sinoval, sphericalcardioid, sphericalellipse, sphericalhelix, sphericalloxodrome, sphericalloxodromeunitspeed, sphericalnephroid, sphericalspiral, spring, strophoid, strophoidimplicit, tanhspiral, teardrop, teeth, tennisballseam, tooth, torusknot, tractrix, tractrixminus, tractrixplus, tractrixunitspeed, triangle, trident, tridentimplicit, trisectrix, tschirnhausen, tschirnhausenpolar, twicubic, twistedn, veryflat, viviani, watt, wattimplicit, wigglyellipse, zcprofile, zetaspiral
*)

BeginPackage["CURVES`","CSPROGS`","PLTPROGS`","Graphics`Graphics`","Graphics`ImplicitPlot`"]

    curves:= Print[Select[Names["CURVES`*"],#!="curves"&]]

(* Parametrically Defined Plane Curves *)

    agnesi::usage="t->agnesi[a][t] is the parametrized curve whose implicit equation is x^2*y==4*a^2*(2*a - y).  Maria Gaetana Agnesi studied the curve 1748, but it had already been studied by Fermat in 1666 and Grandi in in 1703. Agnesi called the curve ``versiera''; this word was mistranslated from Italian into English as ``witch''. To plot try agnesi[1][-1.1,1.1]."

    airyspiral::usage="t->airyspiral[a][t] is a parametrized curve formed using Airy functions.  To plot try airyspiral[1][-20,1]."

    alysoid::usage="s->alysoid[a,b,c][s] is a unit-speed curve which is a generalization of a catenary. Its curvature is -c/(a*b + s^2). s->alysoid[a,a,a][s] is the same as s->catenaryunitspeed[a][s].  See alysoidnds. To plot try alysoid[1,1,3][-1.1,1.1]."

    alysoidkappa2::usage="s->alysoidkappa2[a,b,c][s] is the signed curvature of an alysoid. It is used in the definition of s->alysoidnds[smin,smax][a,b,c][s]."

    alysoidnds::usage="s->alysoidnds[smin,smax][a,b,c][s] is an approximation over the interval smin<s<smax found by NDSolve to s->alysoid[a,b,c][s]." 

    alysoidprime::usage="t->alysoidprime[a,b,c][t] is the velocity of  t->alysoid[a,b,c][t]. It is used in the definitions of s->alysoid[a,b,c][s] and s->alysoidnds[smin,smax][a,b,c][s]."  

    archimedesspiral::usage="t->archimedesspiral[n,a][t] is the spiral of Archimedes of radius a and degree n.  To plot try archimedesspiral[3,1][0,4Pi]."

    astroid::usage="t->astroid[a][t] is the parametrized curve whose implicit equation is x^(2/3) + y^(2/3)==a^(2/3). t->astroid[a,b][t] is the parametrized curve whose implicit equation is (x/a)^(2/3) + (y/b)^(2/3)==1. t->astroid[n,a,b][t] is the parametrized curve whose implicit equation is (x/a)^n + (y/b)^n==1. Jakob Bernoulli studied the astroid in 1691. To plot try astroid[3,1,1][0,2Pi]."

    bellcurve::usage="t->bellcurve[mean,variance][t] is the bell curve of statistics.  The curve originated in essence with De Moivre in 1733. To plot try bellcurve[0,1][-3,3]."

    bicorn::usage="t->bicorn[a][t] is a curve with two horns.  The bicorn or cocked hat was discussed by Sylvester in 1864 and by Cayley in 1867. To plot try bicorn[1][-Pi,Pi]."

    bow::usage="t->bow[a][t] is a parametrized curve whose shape is a bow. The polar equation is r==a*(1 - Tan[theta]^2).  To plot try bow[1][-Pi,Pi]."

    bowtie::usage="t->bowtie[a,b][t] is a parametrized curve whose shape is a bowtie.  To plot try bowtie[1,0.04][0,2Pi]."

    bulletnose::usage="t->bulletnose[a,b][t] is the parametric form of the curve whose implicit equation is a^2/x^2 - b^2/y^2==1.  The bullet-nose curve was discussed by P.H. Schoute in 1885.  The implicit equation is that of a hyperbola with the 2's replaced by -2's. To plot try bulletnose[3,1][-2Pi,2Pi]."

    cardioid::usage="t->cardioid[a][t] is a cardioid or heart-shaped curve that is traced by a point on the circumference of a circle of radius 2*a rolling around a fixed circle of the same radius. It has the polar equation r==2*a*(1 + Cos[theta]). The implicit equation is (x^2 + y^2 - 2*a*x)^2 - 4*a^2*(x^2 + y^2)==0. A cardioid is a special case of an epicycloid.  The cardioid was first conceived by the Dutch mathematician J. Koersma in 1689 an more extensively by Ozanam in 1691.  The English name originated with Castillon on 1741. To plot try cardioid[1][0,2Pi]."

    cardioidunitspeed::usage="s->cardioidunitspeed[a][s] is a unit-speed parametrization of t->cardioid[a][t]. To plot try cardioidunitspeed[1][-8,8]."

    cartesianoval::usage="u->cartesianoval[a,b][u,v] and v->cartesianoval[a,b][u,v] are orthogonal Cartesian ovals. To plot try cartesianoval[1,189.07][-0.5,0.5,0.3]."

    cassini::usage="t->cassini[a,b,pm1,pm2][t] is a parametrization of an oval of Cassini. The parameters pm1 and pm2 are plus or minus 1. The implicit equation is (x^2 + y^2 + a^2)^2 - b^4 - 4*a^2*x^2==0.  Cassinian ovals were first studied by Gian Domenico Cassini in 1680 in connection with the relative motions of the earth and sun. By definition a Cassinian oval is the locus of a point p moving so that the product of the distances of p from two fixed points is constant. A Cassinian oval with a self intersection is a lemniscate. To plot try cassini[1.99,2,1,1][-Pi,Pi]."

    catenary::usage="t->catenary[a][t] is the curve formed by a perfectly flexible inextensible chain of uniform density hanging from two supports. It is given by t->{a*Cosh[t/a],t}.  The curve was found in 1691 by Huygens and Leibniz. To plot try catenary[-1][-2,2]."

    catenaryunitspeed::usage="s->catenaryunitspeed[a][s] is a unit-speed parametrization of t->catenary[a][t]. To plot try catenaryunitspeed[1][-2,2]."

    cayleysextic::usage="t->cayleysextic[a][t] is Cayley's sextic. It is the pedal of a cardioid with respect to its cusp point.  This is the definition of the curve given by Maclaurin in 1718. Cayley made a detailed study of the curves. The evolute of Cayley's sextic is a nephroid. To plot try cayleysextic[1][0,2Pi]."

    circle::usage="t->circle[a][t] is a circle centered at the origin of radius a. The nonparametric form is x^2 + y^2==a^2.  To plot try circle[1][0,2Pi]."

    circleinvolute::usage="t->circleinvolute[n,a][t] is the nth involute starting at {a,0} of the circle t->{a*Cos[t],a*Sin[t]}. t->circleinvolute[a][t] is the same as t->circleinvolute[1,a][t].  To plot try circleinvolute[1,1][0,2Pi]."

    circleinvoluteunitspeed::usage="s->circleinvoluteunitspeed[n,a][s] is a unit-speed parametrization of t->circleinvolute[n,a][t]. s->circleinvoluteunitspeed[a][s] is the same as s->circleinvoluteunitspeed[1,a][s]. To plot try circleinvoluteunitspeed[1,1][0,4Pi]."  

    circleunitspeed::usage="s->circleunitspeed[a][s] is a unit-speed parametrization of t->circle[a][t]. To plot try circleunitspeed[2][0,4Pi]."

    cissoid::usage="t->cissoid[a][t] is an ivy-shaped curve whose implicit equation is x^3 + x*y^2 - 2*a*y^2==0.  Diocles studied the cissoid about 100 B.C. in connection with the classic problem of doubling the cube. Cissoid means ``ivy-shaped''. To plot try cissoid[1][-2,2]."

    cissoidoblique::usage="t->cissoidoblique[a,alpha][t] is a generalization of
cissoid. cissoidoblique[a,0][t] is the same as cissoid[a][Tan[t]]."

    clothoid::usage="s->clothoid[n,a][s] is the nth clothoid. The curvature is proportional to the nth power of the arc length function. It is a unit-speed curve. s->clothoid[1,a][s] is the ordinary clothoid, also known as Euler's spiral. s->clothoid[a][s] is the same as s->clothoid[1,a][s].  See clothoidnds. To plot try clothoid[1,1][-5,5]."  

    clothoidnds::usage="s->clothoidnds[smin,smax][n,a][s] is an approximation over the interval smin<s<smax found by NDSolve to s->clothoid[n,a][s]."  

    clothoidprime::usage="t->clothoidprime[n,a][t] is the velocity of the nth clothoid. It is used in the definitions of s->clothoid[n,a][s] and s->clothoidnds[smin,smax][n,a][s]."

    cnccoprofile::usage="s->cnccoprofile[pm,a,b][s] is a profile curve for a surface of revolution of constant negative curvature -1/a^2 of conic type. It is a unit-speed curve. The parameter pm is plus or minus 1.  To plot try cnccoprofile[1,1,0.3][-1.87,1.87]."

    cnchyprofile::usage="s->cnchyprofile[pm,a,b][s] is a profile curve for a surface of revolution of constant negative curvature -1/a^2 of hyperboloid type. It is a unit-speed curve. The parameter pm is plus or minus 1.  To plot try cnchyprofile[1,1,0.7][-1.15,1.15]."

    cochleoid::usage="t->cochleoid[n,a][t] is the parametrized curve whose polar equation is r==a*(Sin[theta]/theta)^n.  To plot try cochleoid[1,1][-4Pi,4Pi]."

    conchoid::usage="t->conchoid[a,b][t] is a parametrized version of a conchoid of Nicodemes.  The implicit equation is (x^2 + y^2)*(x - b)^2==a^2*x^2.  The curve was invented by Nicodemes around 225 B.C. for the purpose of duplicating a cube. The name conchoid means ``looks like a shell''. The conchoid is more properly the conchoid of a straight line, since more generally the conchoid of any curve can be defined. To plot try conchoid[4,-2][-Pi/2,Pi/2]."

    coshspiral::usage="t->coshspiral[n,a][t]  is a hyperbolic cosine spiral or Poinsot cosh spiral of radius a.  To plot try coshspiral[1,1][-Pi,Pi]."

    cothspiral::usage="t->cothspiral[a][t] is a hyperbolic cotangent spiral of radius a.  To plot try cothspiral[2][-Pi/2,Pi/2]."

    cpcprofile::usage="s->cpcprofile[pm,a,b][s] is a profile curve for a surface of revolution of constant positive curvature 1/a^2. It is a unit-speed curve. The parameter pm is plus or minus 1. There are 3 cases: football type (a>b), sphere (a=b), barrel type (a<b).  To plot try cpcprofile[1,1,0.5][-Pi/2,Pi/2]." 

    crosscurve::usage="t->crosscurve[a,b][t] is the parametric form of the curve whose implicit equation is a^2/x^2 + b^2/y^2==1. The implicit equation is that of an ellipse with the 2's replaced by -2's. To plot try crosscurve[1,2][-Pi,Pi]."

    cubicparabola::usage="t->cubicparabola[a,b,c,d][t] is the general cubic parabola whose implicit equation is y==a*x^3 + b*x^2 + c*x + d.  To plot try cubicparabola[1,0,-1,0][-2,2]."

    cycloid::usage="t->cycloid[a,b][t] is the locus of points traced out by a point on a circle of radius b when a concentric circle of radius a rolls without slipping along a straight line. The cycloid is called prolate if b > a and curtate if a > b. t->cycloid[a][t] is the same as t->cycloid[a,a][t].  Roberval computed the area under a cycloid in 1634. Sir Christopher Wren discover in 1658 that the length of a single arch is four times the diameter of the generating circle. Huygens showed in 1673 that the evolute of a cycloid is also a cycloid; this led to the construction of the cycloidial pendulum. To plot try cycloid[1,2][-4Pi,2Pi]."  

    cycloidunitspeed::usage="s->cycloidunitspeed[a,a][s] is a unit-speed parametrization of t->cycloid[a,a][t]. To plot try cycloidunitspeed[1,1][0,8]."

    delaunay::usage="A Delaunay curve is the path of a focus of a conic with long axis a and eccentricity e as the conic rolls on a straight line. If 0<e<1, the conic is an ellipse; if e>1, the conic is a hyperbola.  Its numerical implementation is delaunaynds.  To plot try delaunay[1,2][-4Pi,4Pi]."

    delaunaykappa2::usage="s->delaunaykappa2[a,e][s] is the signed curvature of a Delaunay curve. It is used in the definition of s->delaunaynds[smin,smax][a,e][s]."

    delaunaynds::usage="s->delaunaynds[smin,smax][a,e][s] is an approximation over the interval smin<s<smax found by NDSolve to a unit-speed parametrization of a Delaunay curve. A Delaunay curve is the path of a focus of a conic with long axis a and eccentricity e as the conic rolls on a straight line. If 0<e<1, the conic is an ellipse; if e>1, the conic is a hyperbola."

    deltoid::usage="t->deltoid[a][t] is a parametrized curve which resembles a delta. A deltoid is a special case of an hypocycloid.  To plot try deltoid[1][0,2Pi]."

    deltoidbis::usage="t->deltoidbis[a][t] is a parametrized curve which resembles a delta. A deltoid is a special case of an hypocycloid.  To plot try deltoidbis[1][0,2Pi]."

    deltoidinvolute::usage="t->deltoidinvolute[a][t] is the involute of t->deltoid[a][t] starting at t=0.  The curve has constant width. To plot try deltoidinvolute[1][-2Pi,2Pi]."

    deltoidunitspeed::usage="s->deltoidunitspeed[a][s] is a unit-speed parametrization of t->deltoid[a][t]. To plot try deltoidunitspeed[1][0,Pi]."

    diamond::usage="t->diamond[a][t] is a parametrized curve whose shape is a diamond. t->diamond[n,a,b][t] is a generalization of t->diamond[a][t] that resembles t->astroid[n,a,b][t].  To plot try diamond[0.5,1,1.5][0.1,6.38]."

    eight::usage="t->eight[a,b][t] is a parametrized curve which resembles a figure eight.  Both t->eight[1,1][t] and t->lemniscate[1][t] resemble a figure eight. Although the definition of t->eight[1,1][t] is simpler than that of  t->lemniscate[1][t], the curvature of the latter curve has fewer peaks and valleys.  To plot try eight[1,1][-Pi,Pi]."

    elasticainflect::usage="s->elasticainflect[a,k][s] is a unit-speed inflectional elastica starting at {0,0}.  To plot try elasticainflect[1,0.8][-8,8]."

    elasticanoninflect::usage="s->elasticanoninflect[a,k][s] is a unit-speed noninflectional elastica starting at {0,0}.  To plot try elasticanoninflect[1,0.8][-8,8]."

    ellipse::usage="t->ellipse[a,b][t] is the standard parametrization of an ellipse centered at the origin using Cos and Sin. The nonparametric form is x^2/a^2 + y^2/b^2==1.  To plot try ellipse[1,2][0,2Pi]."

    ellipsebis::usage="t->ellipsebis[a,b][t] is a parametrization of an ellipse centered at the origin using Sech and Tanh. The nonparametric form is x^2/a^2 + y^2/b^2==1.  To plot try ellipsebis[1,2][-2,2]."

    ellipseinc::usage="t->ellipseinc[a,b,c,psi][t] is the locus of points at a distance c along a line meeting an t->ellipse[a,b] at the point ellipse[a,b][t] with an angle psi.  To plot try ellipseinc[1,2,-1,Pi/4][0,3Pi]."

    epicycloid::usage="t->epicycloid[a,b][t] is a parametrized curve traced out by a point P on the circumference of a circle (of radius b) rolling outside another circle (of radius a). Special cases are: cardioid (b=a), nephroid (b=a/2).  To plot try epicycloid[5,2][0,4Pi]."

    epicycloidunitspeed::usage="s->epicycloidunitspeed[a,b][s] is a unit-speed parametrization of t->epicycloid[a,b][t]. To plot try epicycloidunitspeed[5,2][0,7Pi]."

    epispiral::usage="t->epispiral[n,a][t] is an epispiral. The curve has n branches if n is odd and 2*n branches if n is even. Each branch has two asymptotes radiating out from the center. The polar equation is r==a/Sin[n*theta].  To plot try epispiral[5,1][0,2Pi]."

    epitrochoid::usage="t->epitrochoid[a,b,h][t] is a parametrized curve traced out by a point P attached to a circle (of radius b) rolling outside another circle (of radius a). Let h be the distance from P to the center of the moving circle. Special cases are: circle (a=0), limacon (a=b), epicycloid (h=b), cardioid (h=b=a), nephroid (h=b=a/2). There are a/b loops. For an epicycloid the loops degenerate to points. To plot try epitrochoid[7,1,2][0,2Pi]."

    fermatspiral::usage="t->fermatspiral[a][t] is a spiral of Fermat of radius a, where t is positive. t->fermatspiral[n,a][t] is a spiral of Fermat of radius a and degree n, where t may be positive or negative. Special cases are: fermatspiral (n=1/4), lituus (n=3/4), archimedesspiral (n=(m - 1)/(2*m)).  To plot try fermatspiral[0.25,1][-4Pi,4Pi]."

    folium::usage="t->folium[t] is the folium of Descartes. The nonparametric form is x^3 + y^3 - 3*x*y==0. To plot try Show[folium[-0.999,200],folium[-60,-1.001],PlotRange->{{-3,2},{-3,2}}]." 

    freeth::usage="t->freeth[n,a][t] is Freeth's nephroid of order n. To plot try freeth[2,1][0,4Pi]."
    
    genepicycloid::usage="t->genepicycloid[coef]][t] is a generalized epicycloid generated by the list coef genepicycloid[{n*a,0,...,0,-a}]][t]
coincides with epicycloid[(n - 1)*a,a][t]."

    genfolium::usage="t->genfolium[n,a,b,c][t] is a generalized folium of Descartes. genfolium[1,1,1,1][t] is the same as folium[t]. The nonparametric form is a*x^(2*n + 1) + b*y^(2*n + 1) - c*x^n*y^n==0. To plot try genfolium[2,1,1,1][-0.999,60]."

    genparabola::usage="t->genparabola[n][t] is the curve t->{t,t^n}.  To plot try genparabola[3][-1,1]."

    hankelspiral::usage="t->hankelspiral[n,a][t] is a Hankel spiral of order n.  To plot try hankelspiral[1,1][0.1,14]."

    hippias::usage="t->hippias[a,b][t] is a quadratrix of Hippias (290-450).  To plot try hippias[1,1][-Pi,Pi]."

    hippopede::usage="t->hippopede[a,b][t] is the parametric form of a curve whose shape resembles a horse fetter.  To plot try hippopede[1.1,1][-Pi,Pi]."

    hyperbola::usage="t->hyperbola[a,b][t] is the standard parametrization of a hyperbola centered at the origin using Cosh and Sinh. The nonparametric form is x^2/a^2 - y^2/b^2==1.  To plot try hyperbola[1,2][-2,2]."

    hyperbolabis::usage="t->hyperbolabis[a,b][t] is a parametrization of a hyperbola centered at the origin using Sec and Tan. The nonparametric form is x^2/a^2 - y^2/b^2==1.  To plot try hyperbolabis[1,2][-1,1]."

    hyperbolicspiral::usage="t->hyperbolicspiral[a][t] is a hyperbolic spiral of radius a. It is a special case of a spiral of Archimedes.  The hyperbolic spiral was first studied by Varigon in 1704 and Cotes in 1722. To plot try hyperbolicspiral[1][0.03Pi,6Pi]."

    hypocycloid::usage="t->hypocycloid[a,b][t] is a parametrized curve traced out by a point P on the circumference of a circle (of radius b) rolling inside another circle (of radius a). Special cases are: deltoid (b=a/3), astroid (b=a/4).  To plot try hypocycloid[8,1][0,2Pi]."

    hypocycloidinvolute::usage="t->hypocycloidinvolute[a,b][t] is the involute of t->hypocycloid[a,b][t] starting at t=0.  To plot try hypocycloidinvolute[5,1][0,4Pi]."

    hypocycloidunitspeed::usage="s->hypocycloidunitspeed[a,b][s] is a unit-speed parametrization of t->hypocycloid[a,b][t]. To plot try hypocycloidunitspeed[8,1][0,2Pi]."

    hypotrochoid::usage="t->hypotrochoid[a,b,h][t] is a parametrized curve traced out by a point P on one circle (of radius b) rolling inside another circle (of radius a). Let h be the distance from P to the center of the moving circle. Special cases are: ellipse (a=2*b), hypocycloid  (h=b), deltoid (h=b=a/3), astroid (h=b=a/4). There are a/b loops when a/b is an integer larger than 2 and h>b. For a hypocycloid the loops degenerate to points.  To plot try hypotrochoid[5,1,2][0,2Pi]."

    kampyle::usage="t->kampyle[a][t] is the kampyle or crooked staff of Eudoxus, who studied the curve around 365 B.C.  The polar equation is r==a*Sec[theta]^2.  To plot try kampyle[1][-Pi,Pi]."

    kappacurve::usage="t->kappacurve[a][t] is a kappa or Gutschoven curve. Gutschoven was student of Descartes; he studied the curve in 1662. To plot try kappacurve[1][0,2Pi]."

    keplerorbit::usage="t->keplerorbit[a,e][t] is a Kepler orbit, where e is the eccentricity. The orbit is part of an ellipse (0<=e<1), parabola (e=1) or hyperbola (e>1).  To plot try keplerorbit[2,0.7][0,2Pi]."
 
    lehr::usage="A Lehr curve is a unit-speed plane curve whose signed curvature is a + b*Cos[c*s], or more generally,  a + b1*Cos[c*s] + b2*Cos[2*c*s] +.... Lehr curves were studied systematically by Eduard Lehr in 1932. The curve is numerically by lehrnds.  To plot try lehr[1,5,3][0,2Pi] or lehr[1,{3,4},5][0,2Pi]."

    lehrkappa2::usage="s->lehrkappa2[a,b,c][s] is the signed curvature of a Lehr curve. It is used in the definition of s->lehrnds[smin,smax][a,b,c][s]."
    
    lehrnds::usage="s->lehrnds[smin,smax][a,b,c][s] is an approximation over the interval smin<s<smax found by NDSolve to a unit-speed parametrization of a Lehr curve. Its signed curvature is a + b*Cos[c*s]. s->lehrnds[smin,smax][a,coef,c][s] is an approximation over the interval smin<s<smax found by NDSolve to a unit-speed parametrization of a generalized Lehr curve, where coef is a list of coefficients."

    lemniscate::usage="t->lemniscate[a][t] is the parametrized curve whose implicit equation is (x^2 + y^2)^2==a^2*(x^2 - y^2).  Lemniscates were first studied by Jakob Bernoulli in 1694. A lemniscate is a Cassinian oval with a self intersection. To plot try lemniscate[1][0,2Pi]."

    lemniscatebis::usage="t->lemniscatebis[a][t] is a parametrization of a lemniscate in which theta measures the angle subtended by a point on the curve.  To plot try lemniscatebis[1][-Pi/4,Pi/4]."

    lemniscateunitspeed::usage="s->lemniscateunitspeed[a][s] is a unit-speed parametrization of t->lemniscate[a][t]. To plot try lemniscateunitspeed[1][0,2Pi]."

    limacon::usage="t->limacon[a,b][t] is the parametrized curve whose implicit equation is (x^2 + y^2 - 2*a*x)^2==b^2*(x^2 + y^2). t->limacon[a,2*a][t] coincides with t->cardioid[a][t].  t->limacon[a,b][t] coincides with t->pedal[{-2*a,0},circle[b]][t] + {2*a,0}.  The limacon was discovered by Etienne Pascal (father of Blaise) and named by Roberval in 1650. When the curve has a loop, there is a resemblance to a snail (=limacon in French). To plot try limacon[-1,1][0,2Pi]."

    line::usage="t->line[a1,b1,a2,b2][t] is a parametrized line in R^2 with line[a1,b1,a2,b2][0]={a2,b2} and line[a1,b1,a2,b2]'[0]={a1,b1}. To plot try line[1,1,1,1][-3,3]."

    linearpursuit::usage="t->linearpursuit[a,k][t] is a linear pursuit curve.  The curve was first discussed by Bouguer in 1732. To plot try linearpursuit[1,1.5][0,1]."

    lissajous::usage="t->lissajous[n,d,a,b][t] is a Lissajous or Bowditch curve.  The curve was studied by Bowditch in 1815 and by  Lissajous in 1857. To plot try lissajous[3/4,0,9,8][0,8Pi]."

    lituus::usage="t->lituus[a][t] is the lituus spiral of radius a. It is a trumpet-shaped curve that is the locus of points P such that the square of distance of P from the origin is inversely proportional to the angle theta that p makes with the horizontal axis. The polar equation is r==a*/Sqrt[theta].  The curve was first studied by Cotes in 1722. To plot try lituus[1][0,12Pi]."

    logistic::usage="t->logistic[k,a,b][t] is the logistic curve used to model growth in biological populations for which saturation occurs.  To plot try logistic[3,2,-2][-3,3]."

    logspiral::usage="t->logspiral[a,b][t] is the logarithmic spiral or equiangular spiral of radius a and slant b. The polar equation is r==a*E^(b*theta).  The radius vectors make equal angles with one another and are in geometric progression. The stereographic projection of a logspiral onto the sphere S^3 is a loxodrome. Descartes (in 1638) and Toricelli (in 1641) first studied. To plot try logspiral[1,0.04][0,12Pi]."

    logspiralunitspeed::usage="s->logspiralunitspeed[a,b][s] is a unit-speed parametrization of t->logspiral[a,b][t]. To plot try logspiralunitspeed[1,0.04][0,12Pi]."

    nephroid::usage="t->nephroid[a][t] is a nephroid or kidney-shaped curve that is traced by a point on the circumference of a circle of radius a rolling around a fixed circle of radius 2*a. It has the polar equation r==2*a*(Sin[theta/2]^(2/3) + Cos[theta/2]^(2/3))^(3/2). The implicit equation is (x^2 + y^2 - 4*a^2)^3 - 108*a^4*y^2==0. A nephroid is a special case of an epicycloid.  Huygens, Tschirnhausen and Bernoulli studied the nephroid. To plot try nephroid[1][0,2Pi]."

    nephroidunitspeed::usage="s->nephroidunitspeed[a][s] is a unit-speed parametrization of t->nephroid[a][t]. To plot try nephroidunitspeed[1][-6,6]."

    ngon::usage="t->ngon[n,a][t] is a parametrization of a regular polygon with n sides, where a is the length of the apothem.  To plot try ngon[5,1][0,1]."

    nielsenspiral::usage="t->nielsenspiral[a][t] is Nielsen's spiral of radius a.  To plot try nielsenspiral[1][0.1,12Pi]."

    parabola::usage="t->parabola[a][t] is the vertical parabola with vertex at the origin, focus at {0,a} and directrix t->{t,-a}. The nonparametric form is x^2==4*a*y.  To plot try parabola[1][-1.5,1.5]."

    parabolicspiral::usage="t->parabolicspiral[n,a,b][t] is a parabolic spiral of radius a, offset b and degree n, where t may be positive or negative. Special cases are: fermatspiral (b=0), lituus (n=3/4,b=0), archimedesspiral (n=(m - 1)/(2*m),b=0).  To plot try parabolicspiral[0.25,1,1][-6Pi,2Pi]."

    piriform::usage="t->piriform[a,b][t] is a pear-shaped curve whose implicit equation is  a^4*y^2 - b^2*x^3*(2*a - x)==0.  The curve was studied by De Longchamps in 1886. To plot try piriform[1,1][0,2Pi]."

    pseudocatenary::usage="s->pseudocatenary[a,b][s] is a unit-speed curve which is a generalization of a catenary. Its curvature is -a/(a*b + s^2). s->pseudocatenary[a,a][s] is the same as s->catenaryunitspeed[a][s].  See pseudocatenarynds. To plot try pseudocatenary[1,0.3][-1,1]."

    pseudocatenarykappa2::usage="s->pseudocatenarykappa2[a,b][s] is the signed curvature of a pseudocatenary. It is used in the definition of s->pseudocatenarynds[smin,smax][a,b,c][s]."

    pseudocatenarynds::usage="s->pseudocatenarynds[a,b][s] is an approximation over the interval smin<s<smax found by NDSolve to s->pseudocatenary[a,b][s]."

    pseudocatenaryprime::usage="t->pseudocatenaryprime[a,b][t] is the velocity of t->pseudocatenary[a,b][t]."

    rectangle::usage="t->rectangle[n,a,b][t] is a generalized rectangle. t->rectangle[1,a,b][t] is an ellipse. t->rectangle[2,a,b][t] is a parametrization of an ordinary rectangle. t->rectangle[3,a,b][t] is an astroid rotated by -Pi/4.  To plot try rectangle[2,1,2][0.1,6.38]."

    reuleauxpolygon::usage="t->reuleauxpolygon[n,a][t] is a parametrization of a Reuleaux constant width curve with 2*n +1 sides, where a is the distance between any vertex and the center.  To plot try reuleauxpolygon[3,1][0,1]."

    riemannfractal::usage="t->riemannfractal[k,a,m][t] is a Riemann fractal.  To plot try riemannfractal[2,2,15][0,1]."

    rose::usage="t->rose[n,a][t] is a rose of radius a with n petals if n is odd and 2*n petals if n is even. The polar equation is r==a*Cos[n*theta].  To plot try rose[3,1][0,Pi]."

    sc::usage="t->sc[a][t] is the semicubical parabola t->{(a*t)^2,(a*t)^3}. The nonparametric form is y^2==x^3.  To plot try sc[1][-2,2]."

    scarab::usage="t->scarab[a,b][t] is a parametrization of a scarab.  To plot try scarab[1,1][0,2Pi]."

    scunitspeed::usage="s->scunitspeed[s] is a unit-speed parametrization of t->sc[1][t]. To plot try scunitspeed[0,2]."

    semicubic::usage="z->semicubic[a,b,c,d][z] is a parametrization of the semicubic y^2 = a*x^3 + b*x^2 + c*x + d using the Weierstrass P function and its derivative.  To plot try semicubic[2,2,0,1][-5,5,PlotRange->{{-1.5,1.5},{-2,2}}]."

    serpentine::usage="t->serpentine[a,b][t] is a serpentine curve.  Serpentines were first studied by Newton in 1701 in his great paper on cubics. To plot try serpentine[2,2][-3,3]."

    sinhspiral::usage="t->sinhspiral[n,a][t] is a hyperbolic sine spiral or Poinsot sinh spiral of radius a.  To plot try sinhspiral[2,1][0.01,Pi-0.01]."

    sinoval::usage="t->sinoval[n,a][t] is the oval t->{a*Cos[t],a*Sin[Sin[...Sin[t]]]}.  To plot try sinoval[3,1][0,2Pi]."

    spring::usage="t->spring[a,b,v][t] is a plane curve which resembles a spring. It is the projection of t->helix[a,b][t] onto a plane whose normal makes an angle v with the z-axis.  To plot try spring[1,0.03,Pi/4][0,30Pi]."

    strophoid::usage="t->strophoid[a][t] is the parametrized curve whose implicit equation is y^2*(a - x)==x^3. t->strophoid[m,a][t] is the parametrized curve whose implicit equation is y^2*(a - x)==x^2*(a*m + x). The case m=0 is the cissoid t->cissoid[a/2][t] and the case m=1 is the strophoid t->strophoid[a][t].  To plot try strophoid[1,1][-5,5]."

    talbot::usage="t->talbot[a,b,f][t] is Talbot's curve.  To plot try talbot[1,1,0.9][-Pi,Pi]."

    tanhspiral::usage="t->tanhspiral[a][t] is the hyperbolic tangent spiral of radius a.  To plot try tanhspiral[1][-Pi/2,Pi/2]."

    teardrop::usage="t->teardrop[a,b,c,d,n][t] is a curve resembling a teardrop that has high order contact with the horizontal axis.  To plot try teardrop[8,1,2,1,2][-Pi,Pi]."

    tractrix::usage="t->tractrix[a][t] is a tractrix. It represents the path of a particle P pulled by an inextensible string whose end moves along a line. The involute of a catenary is a tractrix. Huygens studied the tractrix in 1692. A tractrix is the profile curve of a surface of revolution of constant Gaussian curvature -1/a^2. t->tractrix[a,b][t] is a generalization of a tractrix in the same way that an ellipse is a generalization of a circle. In the terminology of Bianchi, t->tractrix[a,b][t] is a shortened tractrix for b<a or extended tractrix for b>a.  To plot try tractrix[1][0.01,Pi-0.1]."

    tractrixminus::usage="t->tractrixminus[a][t] is a parametrization of the bottom half of t->tractrix[a][t].  To plot try tractrixminus[1][0.01,1]."

    tractrixplus::usage="t->tractrixplus[a][t] is a parametrization of the top half of t->tractrix[a][t].  To plot try tractrixplus[1][0.01,1]."

    tractrixunitspeed::usage="s->tractrixunitspeed[a][s] is a unit-speed parametrization of t->tractrix[a][t].  To plot try tractrixunitspeed[1][0,2]."

    triangle::usage="t->triangle[t] is a parametrization of an equilateral triangle.  To plot try triangle[0,3]."

    trident::usage="t->trident[a,b,c,d][t] is the parametric form of the curve whose implicit equation is x*y==a*x^3 + b*x^2 + c*x + d.  In  classical mythology a trident is a 3-pronged spear serving in as the attribute of a sea god. The curve was named by Newton in 1701, having been previously studied by Descartes in 1637. To plot try trident[1,0,0,1][-1,1]."

    trisectrix::usage="t->trisectrix[a][t] is the trisectrix of Maclaurin, who studied the curve in 1742.  To plot try trisectrix[1][-0.49Pi,0.49Pi]."

    tschirnhausen::usage="t->tschirnhausen[n,a][t] is a parametrization of Tschirnhausen's cubic.  To plot try tschirnhausen[1,1][-1.4Pi,1.4Pi]."

    watt::usage="t->watt[a,b,c][t] is a parametrization of Watt's curve.  To plot try watt[1,2,1][-2Pi,2Pi]."

    zcprofile::usage="t->zcprofile[pm,a,c][t] is the profile curve of a generalized helicoid of zero Gaussian curvature. The parameter pm is +1 or -1.  To plot try zcprofile[1,1,1][1,3]."

    zetaspiral::usage="t->zetaspiral[a][t] is the spiral t->{a*Re[Zeta[1/2 + I*t]],a*Im[Zeta[1/2 + I*t]]}."
    
(* Implicitly Defined Plane Curves *)

    agnesiimplicit::usage="agnesiimplicit[a][x,y]==0 is the nonparametric form of the parametrized curve t->{2*a*Tan[t],2*a*Cos[t]^2}. To plot try agnesiimplicit[1][{-5,5}]."

    astroidimplicit::usage="astroidimplicit[n,a,b][x,y]==0 is the nonparametric form of the parametrized curve t->{a*Cos[t]^n,b*Sin[t]^n}."

    bulletnoseimplicit::usage="bulletnoseimplicit[a,b][x,y]==0 is the nonparametric form of the parametrized curve t->{a*Cos[t],b*Cot[t]}. To plot try bulletnoseimplicit[2,2][{-2.01,2}]."

    cardioidimplicit::usage="cardioidimplicit[a][x,y]==0 is the nonparametric form of the parametrized curve t->{a*(1 + Cos[t])*Cos[t],a*(1 + Cos[t])*Sin[t]}.  To plot try cardioidimplicit[1][{-0.5,4.1}]."

    cartesianimplicit::usage="cartesianimplicit[a,c,m][x,y]==0 is an oval of Descartes: ((1 - m^2)*(x^2 + y^2) + 2*m^2*c*x + a^2 - m^2*c^2)^2 - 4*a^2*(x^2 + y^2)==0. To plot try cartesianimplicit[1,3,2][{-2,8}]."

    cassiniimplicit::usage="cassiniimplicit[a,b][x,y]==0 is an oval of Cassini: (x^2 + y^2 + a^2)^2 - b^4 - 4*a^2*x^2==0. To plot try cassiniimplicit[2.01,2][{-3,3}]."

    circleimplicit::usage="circleimplicit[a][x,y]==0 is the nonparametric form of the parametrized circle t->{a*Cos[t],a*Sin[t]}. To plot try circleimplicit[1][{-1,1}]."

    cissoidimplicit::usage="cissoidimplicit[a][x,y]==0 is the nonparametric form of the parametrized curve t->{2*a*t^2/(t^2 + 1),2*a*t^3/(t^2 + 1)}. To plot try cissoidimplicit[1][{0,1}]."

    conchoidimplicit::usage="conchoidimplicit[a,b][x,y]==0 is the nonparametric form of a conchoid of Nicodemes.  The parametrized version is t->{b + a*Cos[t],b*Tan[t] + a*Sin[t]}. To plot try conchoidimplicit[4,1][{-4,5}]."

    crosscurveimplicit::usage="crosscurveimplicit[a,b][x,y]==0 is the nonparametric form of the parametrized curve t->{a*Sec[t],b*Csc[t]}. To plot try crosscurveimplicit[1,1][{-4,4}]."

    cubicparabolaimplicit::usage="cubicparabolaimplicit[a,b,c,d][x,y]==0 is the nonparametric form of the parametrized curve t->{t,a*t^3 + b*t^2 + c*t + d}. To plot try cubicparabolaimplicit[1,-2,-1,0][{-2,3}]."

    deltoidimplicit::usage="deltoidimplicit[a][x,y]==0 is the nonparametric form of the parametrized curve t->{2*a*Cos[t]*(1 + Cos[t]) - a,2*a*Sin[t]*(1 - Cos[t])}. To plot try deltoidimplicit[1][{-1.5,3}]."

    devilimplicit::usage="devilimplicit[a,b][x,y]==0 is the devil's curve: y^2*(y^2 - b^2) - x^2*(x^2 - a^2)==0. To plot try devilimplicit[1.5,1][{-2,2}]."

    duererimplicit::usage="duererimplicit[a,b][x,y]==0 is Duerer's conchoid. Let q + r = b. The locus of points {x,y} whose distance from {q,0} is a and lie on the line through {q,0} and {r,0} constitute Duerer's conchoid. To plot try duererimplicit[1,1][{-10,10}]."

    ellipseimplicit::usage="ellipseimplicit[a,b][x,y]==0 is the nonparametric form of the parametrized curve t->{a*Cos[t],b*Sin[t]}. To plot try ellipseimplicit[2,1][{-2,2}]."

    foliumimplicit::usage="foliumimplicit[x,y]==0 is the nonparametric form of the parametrized curve t->{3*t/(1 +t^3),3*t^2/(t^3 + 1)}. foliumimplicit[n,a,b,c][x,y]==0 is the nonparametric form of the parametrized curve t->{c*t^n/(a + b*t^(2*n + 1)),c*t^(n + 1)/(a + b*t^(2*n + 1))}. To plot try foliumimplicit[5,2,-1,1][{-1,1}]."

    hippopedeimplicit::usage="hippopedeimplicit[a,b][x,y]==0 is the nonparametric form of a curve whose shape resembles a horse fetter. to plot try hippopedeimplicit[1.1,1][{-3,3}]."

    hyperbolaimplicit::usage="hyperbolaimplicit[a,b][x,y]==0 is the nonparametric form of the parametrized curve t->{a*Cosh[t],b*Sinh[t]}. To plot try hyperbolaimplicit[1,2][{-2,2}]."

    kappacurveimplicit::usage="kappacurveimplicit[a][x,y]==0 is the nonparametric form of a kappa or Gutschoven curve. To plot try kappacurveimplicit[2][{-2,2}]."

    keplerimplicit::usage="keplerimplicit[a,b][x,y]==0 is a folium of Kepler: ((x - b)^2 + y^2)*(x*(x - b) + y^2) - 4*a*(x - b)*y^2==0. To plot try keplerimplicit[2,1][{-1,3}]."

    lemniscateimplicit::usage="lemniscateimplicit[a][x,y]==0 is the nonparametric form of the parametrized curve t->{a*Cos[t]/(1 + Sin[t]^2),a*Sin[t]*Cos[t]/(1 + Sin[t]^2)}. lemniscateimplicit[b*Sqrt[2]][x,y] coincides with cassiniimplicit[b,b][x,y]. To plot try lemniscateimplicit[1][{-1,1}]."

    limaconimplicit::usage="limaconimplicit[a,b][x,y]==0 is the nonparametric form of the parametrized curve t->{(2*a*Cos[t] + b)*Cos[t]*(2*a*Cos[t] + b)*,Sin[t]}. To plot try limaconimplicit[2,1][{-1,5}]."

    nephroidimplicit::usage="nephroidimplicit[a][x,y]==0 is the nonparametric form of the parametrized curve t->{a*(3*Cos[t] - Cos[3*t]),a*(3*Sin[t] - Sin[3*t])}. To plot try nephroidimplicit[1][{-3,3}]."

    newtonphillipsimplicit::usage="newtonphillipsimplicit[m,n,a,b][x,y]==0 is the nonparametric curve studied by Newton and Phillips."

    parabolaimplicit::usage="parabolaimplicit[a][x,y]==0 is the nonparametric form of the parametrized curve t->{2*a*t,a*t^2}. To plot try parabolaimplicit[1][{-3,3}]."

    perseusimplicit::usage="perseusimplicit[a,b,c][x,y]==0 is the nonparametric form of a spiric of Perseus. perseusimplicit[a,b^2/(2*a),b^2/(2*a)][x,y]==0 coincides with cassiniimplicit[a,b][x,y]==0. A section of a torus by a plane parallel to the axis of revolution is a spiric of Perseus. To plot try perseusimplicit[5,3,1.9][{-11,11}]."

    piriformimplicit::usage="piriformimplicit[a,b][x,y]==0 is the nonparametric form of the parametrized curve t->{(1 + Sin[t])*a,(1 + Sin[t])*b*Cos[t]}. To plot try piriformimplicit[1,1][{0,2}]."

    predatorpreyimplicit::usage="predatorpreyimplicit[a,b,p,q][x,y] represents the implicit solution of the predator-prey system x'==x*(a - b*y),y'==y*(-p + q*x) centered at {p/q,a/b}.  To plot try predatorpreyimplicit[2,1,2,1][0,6,0,6]."

    scimplicit::usage="scimplicit[x,y]==0 is the nonparametric form of the semicubical parabola t->{t^2,t^3}. To plot try scimplicit[{-1,1}]."

    semicubicimplicit::usage="semicubicimplicit[a,b,c,d][x,y]==0 is the nonparametric form of the curve y^2 = a*x^3 + b*x^2 + c*x + d. To plot try semicubicimplicit[1,0,-1,0][{-2,2}]."

    serpentineimplicit::usage="serpentineimplicit[a,b][x,y]==0 is the nonparametric form of the parametrized curve t->{t,a*t/(1 + b*t^2)}. To plot try serpentineimplicit[3,3][{-2,2}]."

    strophoidimplicit::usage="strophoidimplicit[a][x,y]==0 is the nonparametric form of the parametrized curve t->{a*(t^2 - 1)/(t^2 + 1),a*t*(t^2 - 1)/(t^2 + 1)}. strophoidimplicit[m,a][x,y]==0 is the nonparametric form of the parametrized curve t->{a*(t^2 - m)/(t^2 + 1),a*t*(t^2 - m)/(t^2 + 1)}. To plot try strophoidimplicit[1,1][{-2,2}]."

    tridentimplicit::usage="tridentimplicit[a,b,c,d][x,y]==0 is the nonparametric form of the curve x*y==a*x^3 + b*x^2 + c*x + d. To plot try tridentimplicit[1,1,1,1][{-2.001,2}]."

    wattimplicit::usage="wattimplicit[a,b,c][x,y]==0 is the nonparametric form of Watt's curve."

(* Polar Defined Plane Curves *)

    archimedesspiralpolar::usage="theta->archimedesspiralpolar[n,a][theta] is the polar form of a spiral of Archimedes of radius a and degree n. To plot try archimedesspiralpolar[1,1][-3Pi,3Pi]."

    bowpolar::usage="theta->bowpolar[a][theta] is the polar form of a parametrized curve whose shape is a bow. To plot try bowpolar[1][-Pi,Pi]."

    cardioidpolar::usage="theta->cardioidpolar[a][theta] is the polar form of a heart-shaped curve that is traced by a point on the circumference of a circle of radius 2*a rolling around a fixed circle of the same radius. To plot try cardioidpolar[1][0,2Pi]."

    cassinipolar::usage="theta->cassinipolar[a,b,pm1,pm2][theta] is the polar form of an oval of Cassini. The parameters pm1 and pm2 are plus or minus 1. The implicit equation is (x^2 + y^2 + a^2)^2 - b^4 - 4*a^2*x^2==0. To plot try cassinipolar[1,2,1,1][-Pi,Pi]."

    cochleoidpolar::usage="theta->cochleoidpolar[n,a][theta] is the polar form of a cochleoid of degree n. To plot try cochleoidpolar[2,1][-4Pi,4Pi]."

    conchoidpolar::usage="theta->conchoidpolar[a,b][theta] is the polar form of a conchoid of Nicodemes. To plot try conchoidpolar[4,-2][-Pi/2,Pi/2]."

    epispiralpolar::usage="theta->epispiralpolar[n,a][theta] is the polar form of a curve that has n branches if n is odd and 2*n branches if n is even. Each branch has two asymptotes radiating out from the center. To plot try epispiralpolar[4,1][0,2Pi]."

    fermatspiralpolar::usage="theta->fermatspiralpolar[a][theta] is the polar form of the spiral of Fermat of radius a, which may be positive or negative. To plot try fermatspiralpolar[1][0,4Pi]."

    foliumpolar::usage="theta->foliumpolar[theta] is the polar form of a folium of Descartes. The nonparametric form is x^3 + y^3 - 3*x*y==0. theta->->foliumpolar[n,a,b,c][theta->] is the polar form of a generalized folium of Descartes. The nonparametric form is a*x^(2*n + 1) + b*y^(2*n + 1) - c*x^n*y^n==0. To plot try foliumpolar[2,1,1,1][-Pi,Pi]."

    hyperbolicspiralpolar::usage="theta->hyperbolicspiralpolar[a][theta] is the polar form of the hyperbolic spiral of radius a. It is a special case of a spiral of Archimedes. To plot try hyperbolicspiralpolar[1][0.03Pi,9Pi]."

    kampylepolar::usage="theta->kampylepolar[a][theta] is the polar form of the kampyle or crooked staff of Eudoxus. To plot try kampylepolar[1][-Pi,Pi]."

    keplerorbitpolar::usage="theta->keplerorbitpolar[a,e][theta] is the polar form of a Kepler orbit, where e is the eccentricity. The orbit is part of an ellipse (0<=e<1), parabola (e=1) or hyperbola (e>1). To plot try keplerorbitpolar[1,1.2][0,2Pi]."

    lemniscatepolar::usage="theta->lemniscatepolar[a][theta] is the polar form of the curve whose implicit equation is (x^2 + y^2)^2==a^2*(x^2 - y^2). The parameter a is positive or negative. To plot try lemniscatepolar[1][0,2Pi]."

    limaconpolar::usage="theta->limaconpolar[a,b][theta] is the polar form of the curve whose implicit equation is (x^2 + y^2 - 2*a*x)^2==b^2*(x^2 + y^2). theta->limaconpolar[n,a,b][theta_] is a generalization. To plot try limaconpolar[2,-1,1][0,2Pi]."

    lituuspolar::usage="theta->lituuspolar[a][theta] is the polar form of the lituus spiral of radius a, which may be positive or negative. To plot try lituuspolar[-1][0,12Pi]."

    logspiralpolar::usage="theta->logspiralpolar[a,b][theta] is the polar form of the log spiral of radius a and slant b. To plot try logspiralpolar[1,-0.04][0,12Pi]."

    nephroidpolar::usage="theta->nephroidpolar[a][theta] is the polar form of a kidney-shaped curve that is traced by a point on the circumference of a circle of radius a rolling around a fixed circle of radius 2*a. To plot try nephroidpolar[-1][0,Pi]."

    pacman::usage="theta->pacman[n][theta] is the polar form of a pacman-shaped curve of degree n. To plot try pacman[555][0,2Pi]."

    rosepolar::usage="theta->rosepolar[n,a][theta] is the polar form of a rose of radius a with n petals if n is odd and 2*n petals if n is even. To plot try rosepolar[3,-1][0,Pi]."

    tschirnhausenpolar::usage="theta->tschirnhausenpolar[n,a][theta] is a polar parametrization of Tschirnhausen's cubic.  To plot try tschirnhausenpolar[1,1][-1.4Pi,1.4Pi]."

(* Parametrically Defined Space Curves *)

    ast3d::usage="t->ast3d[n,a,b][t] is a space curve analog of the plane curve astroid.  To plot try ast3d[3,1,1][0,2Pi]."

    baseballseam::usage="t->baseballseam[a,c][t] is an approximation to the seam on a baseball of radius a.  To plot try baseballseam[1,0.4][0,4Pi]."

    besselcurve::usage="t->besselcurve[a,b,c][t] is a space curve whose components are Bessel functions.  To plot try besselcurve[-1/2,1/2,3/2][1,30]."

    bicylinder::usage="t->bicylinder[a,b,pm][t] is a parametrization of the intersection of the circular cylinders x^2 + y^2 = a^2 and y^2 + z^2 = b^2. One component is obtained by taking pm = +1, the other by taking pm = -1.  To plot try bicylinder[1,1.01,1][0,2Pi]."

    biquadratic::usage="t->biquadratic[a,k][t] is a parametrization of the intersection of the circular cylinder x^2 + y^2 = a^2 and the elliptic cylinder x^2 + k^2*z^2 = a^2.  To plot try biquadratic[1,0.9][0,12]."

    clelia::usage="t->clelia[a][m,n][t] is a clelia of type (m,n) on a sphere of radius a. The number of petals is n.  To plot try clelia[1][7/8,5][0,Pi]."

    conicalhelix::usage="t->conicalhelix[a,b,c,d][t] is a parametrization of a helix on a cone over t->ellipse[a,b][t]. The slant of the cone is determined by c and the amount of twisting is determined by d.  To plot try conicalhelix[1,1,2,20][-2Pi,2Pi]."

    cubicalellipse::usage="t->cubicalellipse[a,b,c,d,e,f,pm][t] is a cubical ellipse for pm = 1 and a cubical hyperbola for pm = -1.  To plot try cubicalellipse[1,1,3,1,0,1/4,1][-2Pi,2Pi]."

    eightknot::usage="t->eightknot[t] is a parametrization of an eight knot.  To plot try eightknot[0,2Pi]."

    epitrochoid3d::usage="t->epitrochoid3d[a,b,h,omega][t] is the parametrized curve traced out by a point P attached to a circle of radius b rolling along a fixed circle of radius a. The rolling circle lies in a plane whose inclination with the plane of the fixed circle is omega, and h is the distance of P to the center of the rolling circle. The curve t->epitrochoid3d[a,b,h,0][t] is a (planar) epitrochoid and the curve epitrochoid3d[a,b,h,Pi][t] is a (planar) hypotrochoid.  To plot try epitrochoid3d[8,1,4,0.4Pi][0,2Pi]."

    genhelix::usage="t->genhelix[a,b][f][t] is a generalized elliptical helix with slant function f. To plot try genhelix[1,1][3#^0.3&][0,9Pi]."

    helix::usage="t->helix[a,b][t] is a circular helix of radius a and slant b. t->helix[a,b,c][t] is an elliptical helix of slant c.  To plot try helix[1,3,0.5][0,4Pi]."

    helixunitspeed::usage="s->helixunitspeed[a,b][s] is a unit-speed parametrization of a circular helix of radius a and slant b. To plot try helixunitspeed[1,-0.3][0,4Pi]."

    horopter::usage="t->horopter[a,b][t] is a horopter. It is the intersection of the hyperbolic paraboloid y==b*x*z and the cylinder over the curve y^2==x*(2*a - x).  To plot try horopter[1,1][-2,2]."

    line3d::usage="t->line3d[a1,b1,c1,a2,b2,c2][t] is a parametrized line in R^3 with line3d[a1,b1,c1,a2,b2,c2][0]={a2,b2,c2} and line3d[a1,b1,c1,a2,b2,c2]'[0]={a1,b1,c1}."

    lissajous3d::usage="t->lissajous3d[n,d,a,b,c][t] is a space curve analog of Lissajous curve.  To plot try lissajous3d[3/4,0,9,8,4][0,8Pi]."

    oneparametersubgroup::usage="t->oneparametersubgroup[mat,vec][t] is the curve in R^n formed by the one parameter subgroup t->exp(t*mat) acting on a vector v in R^n, where mat is an n by n matrix To plot try oneparametersubgroup[{{0,1,0},{-1,0,0},{1,0,1}},{1,-1,-1/10}][-3Pi,5Pi/4]."

    pseudospheregeodesic::usage="t->pseudospheregeodesic[a,c,u0,pm][t] is a geodesic of slant c on a pseudosphere of constant negative curvature -1/a^2. The parameter pm is plus or minus 1.  To plot try pseudospheregeodesic[1,0.25,Pi,1][0.5Pi,0.919569Pi]."

    pseudosphericalloxodrome::usage="t->pseudosphericalloxodrome[a,b][t] is a loxodrome on a pseudosphere of constant negative curvature -1/a^2.  To plot try pseudosphericalloxodrome[1,4][-2.8,2.8]."

    seiffertspiral::usage="t->seiffertspiral[a,k][t] is Seiffert's spiral of slant k on an sphere of radius a.  To plot try seiffertspiral[1.0,0.9][-3Pi,5Pi/4]."

    sphericalcardioid::usage="t->sphericalcardioid[a][t] is a spherical curve on a sphere of radius 3|a|. It is a curve of constant slope that projects to t->cardioid[-a][t + Pi].  To plot try sphericalcardioid[1/3][0,4Pi]."

    sphericalellipse::usage="t->sphericalellipse[a,b,c,pm][t] is a parametrization of the intersection of the sphere x^2 + y^2 + z^2= a^2 and the elliptical cylinder x^2/b^2 + y^2/c^2 = 1. One component is obtained by taking pm = +1, the other by taking pm = -1.  To plot try sphericalellipse[1,0.8,0.1,1][-3,3]."

    sphericalhelix::usage="t->sphericalhelix[a,b][t] is a spherical helix on a sphere of radius |a + 2*b|. It is a curve of constant slope that projects to t->epicycloid[a,b][t].  To plot try sphericalhelix[3/5,1/5][0,6Pi]."

    sphericalloxodrome::usage="t->sphericalloxodrome[a,b,c][t] is a loxodrome on a sphere of radius a. It is the image under the stereographic injection {u,v}->stereographicsphere[a][u,v] of the logarithmic spiral t->logspiral[b,c][t]. A sphericalloxodrome meets each meridian of the sphere at the same angle.  Note that t->sphericalloxodrome[a,E^d,c][t] is the same as t->mercatorsphere[a][t,d + c*t]. To plot try sphericalloxodrome[1,1,0.2][-20,20]."

    sphericalloxodromeunitspeed::usage="s->sphericalloxodromeunitspeed[a,b,c][s] is a unit-speed parametrization of a loxodrome on a sphere of radius a.  To plot try sphericalloxodromeunitspeed[1,0,0.2][-8,8]."

    sphericalnephroid::usage="t->sphericalnephroid[a][t] is a spherical helix on a sphere of radius 4|a|. It is a curve of constant slope that projects to t->nephroid[a][t].  To plot try sphericalnephroid[1][0,2Pi]."

    sphericalspiral::usage="t->sphericalspiral[a][m,n][t] is a spherical spiral of type (m,n) on a sphere of radius a.  To plot try sphericalspiral[1][24,1][-Pi/2,Pi/2]."

    tennisballseam::usage="t->tennisballseam[n,a,b,c][t] is a curve on the unit sphere in R^3.  It is the normal to t->wigglyellipse[n,a,b,c][t]. For those values of n, a, b and c for which the curve has no self-intersections, the curve (parametrized on 0<=t<=2*Pi) is the boundary of subsets of the unit sphere of equal area.  To plot try tennisballseam[3,1,1,0.1][0,2Pi]."

    torusknot::usage="t->torusknot[a,b,c][p,q][t] is a torus knot of type (p,q) on an elliptical torus of type (a,b,c).  To plot try torusknot[8,3,5][2,5][0,3Pi]."

    twicubic::usage="t->twicubic[t] is the twisted cubic t->{t,t^2,t^3}.  To plot try twicubic[0,1]."

    veryflat::usage="t->veryflat[1][t] is the curve t->{t,0,5*E^(-1/t^2)}. t->veryflat[2][t] is the curve t->If[t<0,{t,5*E^(-1/t^2),0},If[t==0,{0,0,0},{t,0,5*E^(-1/t^2)}]]. These curves both have curvature and torsion zero, but are not congruent by a Euclidean motion.  To plot try veryflat[1][-1,1]."

    viviani::usage="t->viviani[a][t] is a Viviani curve on a sphere of radius 2*a. It is the intersection of the cylinder (x - a)^2 + y^2==a^2 and the sphere x^2 + y^2 + z^2==4*a^2.  To plot try viviani[1][-2Pi,2Pi]."

    wigglyellipse::usage="t->wigglyellipse[n,a,b,c][t] is an ellipse perturbed vertically by a sin curve. The normal to t->wigglyellipse[n,a,b,c][t] divides the unit sphere into two regions of equal area.  To plot try wigglyellipse[3,1,1,0.1][0,2Pi]."


(* Parametrically Defined Curves in R^n*)

    besseln::usage="t->besseln[n,a][t] is the curve t->a*{BesselJ[0,t],...,BesselJ[n - 1,t]}."

    bernstein::usage="t->bernstein[p,q[t]] is the (p,q)-th Bernstein polynomial."

    beziercurve::usage="t->beziercurve[pts][t] is the Bezier curve in R^n determined by pts, where pts is a list of points in R^n."

    cosn::usage="t->cosn[n,a][t] is the curve t->a*{Cos[t],Cos[2*t],...,Cos[n*t]}."

    lagrangecurve::usage="t->hermitecurve[p0,m0,m1,p1][t] is the Hermite curve in R^n determined by p0 and p1 with directions m0 and m1."

    sinn::usage="t->sinn[n,a][t] is the curve t->a*{Sin[t],Sin[2*t],...,Sin[n*t]}."

    twistedn::usage="t->twistedn[n,a][t] is the curve t->a*{t,t^2,...,t^n}."

(* Approximations *)

   regularpolygon::usage="regularpolygon[n,a] is a regular polygon with n vertices circumscribed by a circle of radius a. To plot try Show[Graphics[regularpolygon[7,1]],AspectRatio->Automatic]."

   teeth::usage="teeth[m,t][shortlines] is the mth iterate of tooth of height t."

   tooth::usage="tooth[t][Line[{{a,b},{c,d}}]] constructs from the line Line[{{a,b},{c,d}}] a broken line with a tooth of height t in the middle."

   kochsnowflake::usage="kochsnowflake[m] is a Koch snow flake constructed as the mth iterate of an equilateral triangle using a tooth of height (3 + Sqrt[3])/6. To plot try Show[Graphics[kochsnowflake[2]],AspectRatio->Automatic]."

Begin["`Private`"]

(* Parametrically Defined Plane Curves *)

    agnesi[a_][t_]:= {2*a*Tan[t],2*a*Cos[t]^2}

    agnesi[a_][tmin_,tmax_,opt___]:=
        pplot[agnesi[a][t],{t,tmin,tmax},opt]

    agnesi[{a_,theta_}][t_]:=
    {-2*a*Cos[t]^2*Sin[theta] + 2*a*Cos[theta]*Tan[t],
      2*a*(Cos[t]^2*Cos[theta] + Sin[theta]*Tan[t])}

    agnesi[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[agnesi[{a,theta}][t],{t,tmin,tmax},opt]

    airyspiral[a_][t_]:= {a*AiryAi[t],a*AiryBi[t]}

    airyspiral[a_][tmin_,tmax_,opt___]:=
        pplot[airyspiral[a][t],{t,tmin,tmax},opt]
            
    airyspiral[{a_,theta_}][t_]:=
    {a*(AiryAi[t]*Cos[theta] - AiryBi[t]*Sin[theta]),
     a*(AiryBi[t]*Cos[theta] + AiryAi[t]*Sin[theta])}

    airyspiral[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[airyspiral[{a,theta}][t],{t,tmin,tmax},opt]
            
    alysoid[a_,b_,c_][s_]:=
        Integrate[alysoidprime[a,b,c][ss],{ss,0,s}] + {a,0}

    alysoid/: Derivative[m_?IntegerQ][alysoid[a_,b_,c_]] = 
        Derivative[m - 1][alysoidprime[a,b,c]]

    alysoid[a_,b_,c_][smin_,smax_,opt___]:=
        pplot[alysoidnds[smin,smax][a,b,c][s],{s,smin,smax},opt]

    alysoidkappa2[a_,b_,c_][s_]:= -c/(a*b + s^2)

(*
    alysoidnds[smin_,smax_][a_,b_,c_][s_]:=
         Module[{x,y,tmp,ss},
             tmp=NDSolve[{x'[ss]==alysoidprime[a,b,c][ss][[1]],
               y'[ss]==alysoidprime[a,b,c][ss][[2]],
               x[0]==a,y[0]==0},{x,y},{ss,smin,smax},
               MaxSteps->2000];
               First[{x[s],y[s]} /. tmp]]
*)

    alysoidnds[smin_,smax_][a_,b_,c_][s_]:=
      intrinsic[alysoidkappa2[a,b,c],0,{a,0,Pi/2},
      MaxSteps->2000,{smin,smax}][s]

    alysoidprime[a_,b_,c_][t_]:=
        {Sin[(c/Sqrt[a*b])*ArcTan[t/Sqrt[a*b]]],
         Cos[(c/Sqrt[a*b])*ArcTan[t/Sqrt[a*b]]]}

    archimedesspiral[n_,a_][t_]:=
        {a*t^(1/n)*Cos[t],a*t^(1/n)*Sin[t]}

    archimedesspiral[n_,a_][tmin_,tmax_,opt___]:=
        pplot[archimedesspiral[n,a][t],{t,tmin,tmax},opt]

    archimedesspiral[{n_,a_,theta_}][t_]:=
    {a*t^(1/n)*Cos[t + theta],a*t^(1/n)*Sin[t + theta]}

    archimedesspiral[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[archimedesspiral[{n,a,theta}][t],{t,tmin,tmax},opt]

    astroid[a_][t_]:= {a*Cos[t]^3,a*Sin[t]^3}              

    astroid[a_][tmin_,tmax_,opt___]:=
        pplot[astroid[a][t],{t,tmin,tmax},opt]
            
    astroid[{a_,theta_}][t_]:=
    {a*(Cos[theta]*Cos[t]^3 - Sin[theta]*Sin[t]^3),
     a*(Cos[theta]*Sin[t]^3 + Sin[theta]*Cos[t]^3)}

    astroid[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[astroid[{a,theta}][t],{t,tmin,tmax},opt]
            
    astroid[a_,b_][t_]:= {a*Cos[t]^3,b*Sin[t]^3}

    astroid[a_,b_][tmin_,tmax_,opt___]:=
        pplot[astroid[a,b][t],{t,tmin,tmax},opt]
            
    astroid[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[astroid[{a,b,theta}][t],{t,tmin,tmax},opt]
            
    astroid[{a_,b_,theta_}][t_]:=
    {a*Cos[theta]*Cos[t]^3 - b*Sin[t]^3*Sin[theta],
     b*Cos[theta]*Sin[t]^3 + a*Cos[t]^3*Sin[theta]}

    astroid[n_,a_,b_][t_]:= {a*Cos[t]^n,b*Sin[t]^n}

    astroid[n_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[astroid[n,a,b][t],{t,tmin,tmax},opt]
            
    astroid[{n_,a_,b_,theta_}][t_]:=
    {a*Cos[theta]*Cos[t]^n - b*Sin[t]^n*Sin[theta],
     b*Cos[theta]*Sin[t]^n + a*Cos[t]^n*Sin[theta]}

    astroid[{n_,a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[astroid[{n,a,b,theta}][t],{t,tmin,tmax},opt]
            
    bellcurve[mean_,variance_][t_]:= {t,
        (1/(variance*Sqrt[2*Pi]))*E^(-(t - mean)^2/(2*variance^2))}

    bellcurve[mean_,variance_][tmin_,tmax_,opt___]:=
        pplot[bellcurve[mean,variance][t],
            {t,tmin,tmax},opt,Axes->True,AspectRatio->1/GoldenRatio]

    bellcurve[{mean_,variance_,theta_}][t_]:= 
    {t*Cos[theta] - Sin[theta]/
     (E^((mean - t)^2/(2*variance^2))*Sqrt[2*Pi]* Sqrt[variance^2]),
     Cos[theta]/
     (E^((mean - t)^2/(2*variance^2))*Sqrt[2*Pi]*Sqrt[variance^2]) +
      t*Sin[theta]}

    bellcurve[{mean_,variance_,theta_}][tmin_,tmax_,opt___]:=
        pplot[bellcurve[{mean,variance,theta}][t],
            {t,tmin,tmax},opt,Axes->True,AspectRatio->1/GoldenRatio]

    bicorn[a_][t_]:= {a*Sin[t],a*Cos[t]^2*(2 + Cos[t])/(3 + Sin[t]^2)}

    bicorn[a_][tmin_,tmax_,opt___]:=
        pplot[bicorn[a][t],{t,tmin,tmax},opt]
            
    bicorn[{a_,theta_}][t_]:=
    {a*Cos[theta]*Sin[t] + a*Sin[theta]*Cos[t]^2/(Cos[t] - 2),
     -a*Cos[theta]*Cos[t]^2/(Cos[t] - 2) + a*Sin[theta]*Sin[t]}

    bicorn[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[bicorn[{a,theta}][t],{t,tmin,tmax},opt]
            
    bow[a_][t_]:= {a*(1 - Tan[t]^2)*Cos[t],a*(1 - Tan[t]^2)*Sin[t]}

    bow[a_][tmin_,tmax_,opt___]:=
        pplot[bow[a][t],{t,tmin,tmax},opt]

    bow[{a_,theta_}][t_]:=
    {a*Cos[t + theta]*Cos[2*t]*Sec[t]^2,
     a*Sin[t + theta]*Cos[2*t]*Sec[t]^2}

    bow[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[bow[{a,theta}][t],{t,tmin,tmax},opt]

    bowtie[a_,b_][t_]:= {a*(1 + Cos[t]^2)*Sin[t],
                          (b + Sin[t]^2)*Cos[t]}
  
    bowtie[a_,b_][tmin_,tmax_,opt___]:=
        pplot[bowtie[a,b][t],{t,tmin,tmax},opt]

    bowtie[{a_,b_,theta_}][t_]:=
    {a*(1 + Cos[t]^2)*Cos[theta]*Sin[t] - 
     Cos[t]*(b + Sin[t]^2)*Sin[theta],
     Cos[t]*Cos[theta]*(b + Sin[t]^2) + 
     a*(1 + Cos[t]^2)*Sin[t]*Sin[theta]}

    bowtie[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[bowtie[{a,b,theta}][t],{t,tmin,tmax},opt]

    bulletnose[a_,b_][t_]:= {a*Cos[t],b*Cot[t]}
    
    bulletnose[a_,b_][tmin_,tmax_,opt___]:=
        pplot[bulletnose[a,b][t],
            {t,tmin,tmax},opt,PlotRange->{Automatic,{-5,5}}]

    bulletnose[{a_,b_,theta_}][t_]:=
    {a*Cos[t]*Cos[theta] - b*Cot[t]*Sin[theta],
     b*Cos[theta]*Cot[t] + a*Cos[t]*Sin[theta]}

    bulletnose[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[bulletnose[{a,b,theta}][t],
            {t,tmin,tmax},opt,PlotRange->{Automatic,{-5,5}}]

    cardioid[a_][t_]:= {2*a*Cos[t]*(1 + Cos[t]),
                        2*a*Sin[t]*(1 + Cos[t])} 
    
    cardioid[a_][tmin_,tmax_,opt___]:=
        pplot[cardioid[a][t],{t,tmin,tmax},opt]

    cardioid[{a_,theta_}][t_]:=
    {4*a*Cos[t/2]^2*Cos[t + theta],
     4*a*Cos[t/2]^2*Sin[t + theta]} 
    
    cardioid[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cardioid[{a,theta}][t],{t,tmin,tmax},opt]

    cardioidunitspeed[a_][s_]:=
        {2*a*Cos[2*ArcSin[s/(8*a)]]*(1 + Cos[2*ArcSin[s/(8*a)]]),
         2*a*(1 + Cos[2*ArcSin[s/(8*a)]])*Sin[2*ArcSin[s/(8*a)]]}

    cardioidunitspeed[a_][smin_,smax_,opt___]:=
        pplot[cardioidunitspeed[a][s],{s,smin,smax},opt]

    cardioidunitspeed[{a_,theta_}][s_]:=
    {4*a*(1 - s^2/(64*a^2))*Cos[theta + 2*ArcSin[s/(8*a)]],
     4*a*(1 - s^2/(64*a^2))*Sin[theta + 2*ArcSin[s/(8*a)]]}

    cardioidunitspeed[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[cardioidunitspeed[{a,theta}][s],{s,smin,smax},opt]

    cartesianoval[a_,b_][u_,v_]:=
      {a*Re[WeierstrassP[u + I*v,b,0]],
       a*Im[WeierstrassP[u + I*v,b,0]]}

    cartesianoval[a_,b_][umin_,umax_,v_,opt___]:=
        pplot[cartesianoval[a,b][u,v],{u,umin,umax},opt]

    cassini[a_,b_,pm1_,pm2_][t_]:= 
      {pm1*Sqrt[a^2*Cos[2*t] +
           pm2*Sqrt[b^4 - a^4*Sin[2*t]^2]]*Cos[t],
       pm1*Sqrt[a^2*Cos[2*t] +
           pm2*Sqrt[b^4 - a^4*Sin[2*t]^2]]*Sin[t]}

    cassini[a_,b_,pm1_,pm2_][tmin_,tmax_,opt___]:=
        pplot[cassini[a,b,pm1,pm2][t],{t,tmin,tmax},opt]

    cassini[{a_,b_,pm1_,pm2_,theta_}][t_]:= 
    {pm1*Cos[t + theta]*Sqrt[a^2*Cos[2*t] + 
      pm2*Sqrt[b^4 - a^4*Sin[2*t]^2]],
     pm1*Sin[t + theta]*Sqrt[a^2*Cos[2*t] + 
      pm2*Sqrt[b^4 - a^4*Sin[2*t]^2]]}

    cassini[{a_,b_,pm1_,pm2_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cassini[{a,b,pm1,pm2,theta}][t],{t,tmin,tmax},opt]

    catenary[a_][t_]:= {a*Cosh[t/a],t} 

    catenary[a_][tmin_,tmax_,opt___]:=
        pplot[catenary[a][t],{t,tmin,tmax},opt]

    catenary[{a_,theta_}][t_]:=
    {a*Cos[theta]*Cosh[t/a] - t*Sin[theta],
     t*Cos[theta] + a*Cosh[t/a]*Sin[theta]} 

    catenary[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[catenary[{a,theta}][t],{t,tmin,tmax},opt]

    catenaryunitspeed[a_][s_]:=
        {a*Sqrt[1 + s^2/a^2],a*ArcSinh[s/a]}

    catenaryunitspeed[a_][smin_,smax_,opt___]:=
        pplot[catenaryunitspeed[a][s],{s,smin,smax},opt]

    catenaryunitspeed[{a_,theta_}][s_]:=
    {a*(Sqrt[1 + s^2/a^2]*Cos[theta] - ArcSinh[s/a]*Sin[theta]),
     a*(ArcSinh[s/a]*Cos[theta] + Sqrt[1 + s^2/a^2]*Sin[theta])}

    catenaryunitspeed[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[catenaryunitspeed[{a,theta}][s],{s,smin,smax},opt]

    cayleysextic[a_][t_]:=
    {4*a*Cos[t/2]^3*Cos[3*t/2],4*a*Cos[t/2]^3*Sin[3*t/2]}

    cayleysextic[a_][tmin_,tmax_,opt___]:=
        pplot[cayleysextic[a][t],{t,tmin,tmax},opt]

    cayleysextic[{a_,theta_}][t_]:=
    {4*a*Cos[t/2]^3*Cos[3*t/2 + theta],4*a*Cos[t/2]^3*Sin[3*t/2 + theta]}

    cayleysextic[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cayleysextic[{a,theta}][t],{t,tmin,tmax},opt]

    circle[a_][t_]:= {a*Cos[t],a*Sin[t]}

    circle[a_][tmin_,tmax_,opt___]:=
        pplot[circle[a][t],{t,tmin,tmax},opt]

    circle[{a_,theta_}][t_]:= {a*Cos[t + theta],a*Sin[t + theta]}

    circle[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[circle[{a,theta}][t],{t,tmin,tmax},opt]

    circleinvolute[n_,a_][t_]:= Module[{tmp1,tmp2},
         tmp1=E^(I*tt)Normal[Series[E^(-I*tt),{tt,0,n}]];
         tmp2=a*ComplexExpand[{Re[tmp1],Im[tmp1]}];
         tmp2 /. tt->t]

    circleinvolute[n_,a_][tmin_,tmax_,opt___]:=
        pplot[circleinvolute[n,a][t],{t,tmin,tmax},opt]

    circleinvolute[a_][t_]:=
         {a*(Cos[t] + t*Sin[t]),a*(-t*Cos[t] + Sin[t])}

    circleinvolute[a_][tmin_,tmax_,opt___]:=
        pplot[circleinvolute[a][t],{t,tmin,tmax},opt]

    circleinvolute[{a_,theta_}][t_]:=
         {a*(Cos[t] + t*Sin[t]),a*(-t*Cos[t] + Sin[t])}

    circleinvolute[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[circleinvolute[{a,theta}][t],{t,tmin,tmax},opt]

    circleinvoluteunitspeed[n_,a_][s_]:=
         circleinvolute[n,a][((n + 1)!*s/a)^(1/(n +1))]

    circleinvoluteunitspeed[n_,a_][smin_,smax_,opt___]:=
        pplot[circleinvoluteunitspeed[n,a][s],{s,smin,smax},opt]

    circleinvoluteunitspeed[a_][s_]:=
    {a*(Cos[Sqrt[2]*Sqrt[s/a]] + 
      Sqrt[2]*Sqrt[s/a]*Sin[Sqrt[2]*Sqrt[s/a]]),
     a*(-Sqrt[2]*Sqrt[s/a]*Cos[Sqrt[2]*Sqrt[s/a]] + 
      Sin[Sqrt[2]*Sqrt[s/a]])}

    circleinvoluteunitspeed[a_][smin_,smax_,opt___]:=
        pplot[circleinvoluteunitspeed[a][s],{s,smin,smax},opt]

    circleinvoluteunitspeed[{a_,theta_}][s_]:=
    {a*(Cos[Sqrt[2]*Sqrt[s/a] + theta] + 
      Sqrt[2]*Sqrt[s/a]*Sin[Sqrt[2]*Sqrt[s/a] + theta]),
     a*(-Sqrt[2]*Sqrt[s/a]*Cos[Sqrt[2]*Sqrt[s/a] + theta] + 
      Sin[Sqrt[2]*Sqrt[s/a] + theta])}

    circleinvoluteunitspeed[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[circleinvoluteunitspeed[{a,theta}][s],{s,smin,smax},opt]

    circleunitspeed[a_][s_]:= {a*Cos[s/a],a*Sin[s/a]}

    circleunitspeed[a_][smin_,smax_,opt___]:=
        pplot[circleunitspeed[a][s],{s,smin,smax},opt]

    circleunitspeed[{a_,theta_}][s_]:=
    {a*Cos[s/a + theta],a*Sin[s/a + theta]}

    circleunitspeed[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[circleunitspeed[{a,theta}][s],{s,smin,smax},opt]

    cissoid[a_][t_]:= {2*a*t^2/(t^2 + 1),2*a*t^3/(t^2 + 1)}

    cissoid[a_][tmin_,tmax_,opt___]:=
        pplot[cissoid[a][t],
            {t,tmin,tmax},opt,PlotRange->{{-0.2,1},Automatic}]

    cissoid[{a_,theta_}][t_]:=
    {(2*a*t^2*(Cos[theta] - t*Sin[theta]))/(t^2 + 1),
     (2*a*t^2*(t*Cos[theta] + Sin[theta]))/(t^2 + 1)}

    cissoid[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cissoid[{a,theta}][t],
            {t,tmin,tmax},opt,PlotRange->{{-0.2,1},Automatic}]

    cissoidoblique[a_,alpha_][t_]:=
    {2*a*Sec[alpha - t]*Sin[t]^3, 
     2*a*Sec[alpha - t]*Cos[t]*Sin[t]^2}

    cissoidoblique[a_,alpha_][tmin_,tmax_,opt___]:=
        pplot[cissoidoblique[a,alpha][t],
            {t,tmin,tmax},opt]

    cissoidoblique[{a_,alpha_,theta_}][t_]:=
    {2*a*Sec[alpha - t]*Sin[t]^2*Cos[t + theta], 
     2*a*Sec[alpha - t]*Sin[t]^2*Sin[t + theta]}

    cissoidoblique[{a_,alpha_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cissoidoblique[{a,alpha,theta}][t],
            {t,tmin,tmax},opt]

    clothoid[n_,a_][s_]:=
         Integrate[clothoidprime[n,a][ss],{ss,0,s}]

    clothoid/: Derivative[m_?IntegerQ][clothoid[n_,a_]] = 
        Derivative[m - 1][clothoidprime[n,a]]

    clothoid[n_,a_][smin_,smax_,opt___]:=
        pplot[clothoidnds[smin,smax][n,a][s],
            {s,smin,smax},opt,PlotPoints->300]

    clothoid[a_][s_]:=
         {a*Sqrt[Pi]*FresnelS[s/Sqrt[Pi]],
          a*Sqrt[Pi]*FresnelC[s/Sqrt[Pi]]}

    clothoid[a_][smin_,smax_,opt___]:=
        pplot[clothoidnds[smin,smax][a][s],
            {s,smin,smax},opt,PlotPoints->300]

    clothoid[{a_,theta_}][s_]:=
    {a*Sqrt[Pi]*(Cos[theta]*FresnelS[s/Sqrt[Pi]] - 
      FresnelC[s/Sqrt[Pi]]*Sin[theta]),
     a*Sqrt[Pi]*(Cos[theta]*FresnelC[s/Sqrt[Pi]] + 
      FresnelS[s/Sqrt[Pi]]*Sin[theta])}

    clothoid[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[clothoid[{a,theta}][s],
            {s,smin,smax},opt,PlotPoints->300]

    clothoidnds[smin_,smax_][n_,a_][s_]:=
         Module[{x,y,tmp,ss},
             tmp=NDSolve[{x'[ss]==clothoidprime[n,a][ss][[1]],
               y'[ss]==clothoidprime[n,a][ss][[2]],
               x[0]==0,y[0]==0},{x,y},{ss,smin,smax},
               MaxSteps->2000];
               First[{x[s],y[s]} /. tmp]]

    clothoidprime[n_,a_][t_]:=
         {a*Sin[t^(n + 1)/(n + 1)],a*Cos[t^(n + 1)/(n + 1)]}
       
    cnccoprofile[pm_,a_,b_][s_]:=
        {pm*b*Sinh[s/a],-I*Sqrt[a^2 - b^2]*EllipticE[I*s/a,
                       -b^2/(a^2 - b^2)]}

    cnccoprofile[pm_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[cnccoprofile[pm,a,b][t],{t,tmin,tmax},opt]

    cnccoprofile[{pm_,a_,b_,theta_}][s_]:=
    {I*Sqrt[a^2 - b^2]*EllipticE[I*s/a,
     b^2/(-a^2 + b^2)]*Sin[theta] + b*pm*Cos[theta]*Sinh[s/a], 
     -I*Sqrt[a^2 - b^2]*Cos[theta]*EllipticE[I*s/a,b^2/(-a^2 + b^2)] + 
     b*pm*Sin[theta]*Sinh[s/a]}

    cnccoprofile[{pm_,a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cnccoprofile[{pm,a,b,theta}][t],{t,tmin,tmax},opt]

    cnchyprofile[pm_,a_,b_][s_]:=
        {pm*b*Cosh[s/a],-I*a*EllipticE[I*s/a,-b^2/a^2]}

    cnchyprofile[pm_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[cnchyprofile[pm,a,b][t],{t,tmin,tmax},opt]

    cnchyprofile[{pm_,a_,b_,theta_}][s_]:=
    {b*pm*Cos[theta]*Cosh[s/a] + 
     I*a*EllipticE[I*s/a,-b^2/a^2]*Sin[theta],
     -I*a*Cos[theta]*EllipticE[I*s/a,-b^2/a^2] + 
     b*pm*Cosh[s/a]*Sin[theta]}

    cnchyprofile[{pm_,a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cnchyprofile[{pm,a,b,theta}][t],{t,tmin,tmax},opt]

    cochleoid[n_,a_][t_]:=
       {a*(Sin[t]/t)^n*Cos[t],a*(Sin[t]/t)^n*Sin[t]}

    cochleoid[n_,a_][tmin_,tmax_,opt___]:=
        pplot[cochleoid[n,a][t],{t,tmin,tmax},opt,PlotRange->All]

    cochleoid[{n_,a_,theta_}][t_]:=
    {a*Cos[t + theta]*(Sin[t]/t)^n,a*Sin[t + theta]*(Sin[t]/t)^n}

    cochleoid[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cochleoid[{n,a,theta}][t],{t,tmin,tmax},opt,
        PlotRange->All]

    conchoid[a_,b_][t_]:= {b + a*Cos[t],b*Tan[t] + a*Sin[t]}

    conchoid[a_,b_][tmin_,tmax_,opt___]:=
        pplot[conchoid[a,b][t],{t,tmin,tmax},opt]

    conchoid[{a_,b_,theta_}][t_]:=
    {(b + a*Cos[t])*Sec[t]*Cos[t + theta],
     (b + a*Cos[t])*Sec[t]*Sin[t + theta]}

    conchoid[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[conchoid[{a,b,theta}][t],{t,tmin,tmax},opt]

    coshspiral[n_,a_][t_]:=
    {a*Cos[t]*Sech[n*t],a*Sin[t]*Sech[n*t]}

    coshspiral[n_,a_][tmin_,tmax_,opt___]:=
        pplot[coshspiral[n,a][t],{t,tmin,tmax},opt]

    coshspiral[{n_,a_,theta_}][t_]:=
    {a*Cos[t + theta]*Sech[n*t],a*Sin[t + theta]*Sech[n*t]}

    coshspiral[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[coshspiral[{n,a,theta}][t],{t,tmin,tmax},opt]

    cothspiral[a_][t_]:= {-Sinh[2*t]/(Cos[2*a*t] - Cosh[2*t]),
                            Sin[2*a*t]/(Cos[2*a*t] - Cosh[2*t])}

    cothspiral[a_][tmin_,tmax_,opt___]:=
        pplot[cothspiral[a][t],{t,tmin,tmax},opt,
        PlotRange->{Automatic,{-4,4}}]

    cothspiral[{a_,theta_}][t_]:=
    {-(Sin[2*a*t]*Sin[theta] + Cos[theta]*Sinh[2*t])/
      (Cos[2*a*t] - Cosh[2*t]),
     (Cos[theta]*Sin[2*a*t] - Sin[theta]*Sinh[2*t])/
     (Cos[2*a*t] - Cosh[2*t])}

    cothspiral[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cothspiral[{a,theta}][t],{t,tmin,tmax},opt,
        PlotRange->{Automatic,{-4,4}}]

    cpcprofile[pm_,a_,b_][s_]:=
       {pm*b*Cos[s/a],a*EllipticE[s/a,b^2/a^2]}

    cpcprofile[pm_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[cpcprofile[pm,a,b][t],{t,tmin,tmax},opt]

    cpcprofile[{pm_,a_,b_,theta_}][s_]:=
    {b*pm*Cos[s/a]*Cos[theta] - 
     a*EllipticE[s/a,b^2/a^2]*Sin[theta],
     a*Cos[theta]*EllipticE[s/a,b^2/a^2] + 
     b*pm*Cos[s/a]*Sin[theta]}

    cpcprofile[{pm_,a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cpcprofile[{pm,a,b,theta}][t],{t,tmin,tmax},opt]

    crosscurve[a_,b_][t_]:= {a*Sec[t],b*Csc[t]}
    
    crosscurve[a_,b_][tmin_,tmax_,opt___]:=
        pplot[crosscurve[a,b][t],{t,tmin,tmax},opt]

    crosscurve[{a_,b_,theta_}][t_]:=
    {a*Cos[theta]*Sec[t] - b*Sin[theta]*Csc[t],
     b*Cos[theta]*Csc[t] + a*Sin[theta]*Sec[t]}
    
    crosscurve[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[crosscurve[{a,b,theta}][t],{t,tmin,tmax},opt]

    cubicparabola[a_,b_,c_,d_][t_]:= {t,a*t^3 + b*t^2 + c*t + d}

    cubicparabola[a_,b_,c_,d_][tmin_,tmax_,opt___]:=
        pplot[cubicparabola[a,b,c,d][t],{t,tmin,tmax},opt]

    cubicparabola[{a_,b_,c_,d_,theta_}][t_]:=
    {t*Cos[theta] - (d + t*(c + t*(b + a*t)))*Sin[theta],
     (d + t*(c + t*(b + a*t)))*Cos[theta] + t*Sin[theta]}

    cubicparabola[{a_,b_,c_,d_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cubicparabola[{a,b,c,d,theta}][t],{t,tmin,tmax},opt]

    cycloid[a_,b_][t_]:= {a*t - b*Sin[t],a - b*Cos[t]}

    cycloid[a_,b_][tmin_,tmax_,opt___]:=
        pplot[cycloid[a,b][t],{t,tmin,tmax},opt]

    cycloid[{a_,b_,theta_}][t_]:=
    {a*t*Cos[theta] - b*Sin[t - theta] - a*Sin[theta],
     -b*Cos[t - theta] + a*(Cos[theta] + t*Sin[theta])}

    cycloid[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[cycloid[{a,b,theta}][t],{t,tmin,tmax},opt]

    cycloidunitspeed[a_,a_][s_]:=
         {4*a*ArcSin[Sqrt[s/(8*a)]] - (4*a - s)*Sqrt[(s/a)*(8 - s/a)]/8,
          s - s^2/(8*a)}

    cycloidunitspeed[a_,a_][smin_,smax_,opt___]:=
        pplot[cycloidunitspeed[a,a][s],{s,smin,smax},opt]

    cycloidunitspeed[{a_,a_,theta_}][s_]:=
    {(-a*(4*a - s)*Sqrt[(8*a - s)*s/a^2]*Cos[theta] + 
      32*a^2*ArcSin[Sqrt[s/a]/(2*Sqrt[2])]*Cos[theta] + 
      s*(-8*a + s)*Sin[theta])/(8*a),
     ((8*a - s)*s*Cos[theta] + 
      a*Sqrt[((8*a - s)*s)/a^2]*(-4*a + s)*Sin[theta] + 
      32*a^2*ArcSin[Sqrt[s/a]/(2*Sqrt[2])]*Sin[theta])/(8*a)}

    cycloidunitspeed[{a_,a_,theta_}][smin_,smax_,opt___]:=
        pplot[cycloidunitspeed[{a,a,theta}][s],{s,smin,smax},opt]

    delaunay[a_,e_][smin_,smax_,opt___]:=
        pplot[delaunaynds[smin,smax][a,e][s],
            {s,smin,smax},opt,PlotPoints->300]

    delaunaykappa2[a_,e_][s_]:=
         e*(e - Cos[s/a])/(a*(1 - 2*e*Cos[s/a] + e^2))

    delaunaynds[smin_,smax_][a_,e_][s_]:=
      intrinsic[delaunaykappa2[a,e],0,{0,0,0},
      MaxSteps->2000,{smin,smax}][s]

    deltoid[a_][t_]:= {2*a*Cos[t]*(1 + Cos[t]) - a,
                       2*a*Sin[t]*(1 - Cos[t])}

    deltoid[a_][tmin_,tmax_,opt___]:=
        pplot[deltoid[a][t],{t,tmin,tmax},opt]

    deltoid[{a_,theta_}][t_]:=
    {a*(Cos[2*t - theta] + 2*Cos[t + theta]),
     a*(-Sin[2*t - theta] + 2*Sin[t + theta])}

    deltoid[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[deltoid[{a,theta}][t],{t,tmin,tmax},opt]

    deltoidbis[a_][t_]:=
        {a*(2*Cos[t] - Sin[2*t]),a*(2*Sin[t] - Cos[2*t])}

    deltoidbis[a_][tmin_,tmax_,opt___]:=
        pplot[deltoidbis[a][t],{t,tmin,tmax},opt]

    deltoidbis[{a_,theta_}][t_]:=
    {a*(2*Cos[t + theta] - Sin[2*t - theta]),
     -a*(Cos[2*t - theta] - 2*Sin[t + theta])}

    deltoidbis[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[deltoidbis[{a,theta}][t],{t,tmin,tmax},opt]

    deltoidinvolute[a_][t_]:=
        {(a/3)*(8*Cos[t/2] + 2*Cos[t] - Cos[2*t]),
         (a/3)*(-8*Sin[t/2] + 2*Sin[t] + Sin[2*t])}

    deltoidinvolute[a_][tmin_,tmax_,opt___]:=
        pplot[deltoidinvolute[a][t],{t,tmin,tmax},opt]

    deltoidinvolute[{a_,theta_}][t_]:=
    {(a/3)*(8*Cos[t/2 - theta] + 2*Cos[t + theta] - Cos[2*t - theta]),
     (a/3)*(-8*Sin[t/2 - theta] + 2*Sin[t + theta] + Sin[2*t - theta])}

    deltoidinvolute[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[deltoidinvolute[{a,theta}][t],{t,tmin,tmax},opt]

    deltoidunitspeed[a_][s_]:=
        {2*a*Cos[4*ArcSin[Sqrt[3]*Sqrt[s/a]/4]/3] + 
         a*Cos[8*ArcSin[Sqrt[3]*Sqrt[s/a]/4]/3],
         2*a*Sin[4*ArcSin[Sqrt[3]*Sqrt[s/a]/4]/3] - 
         a*Sin[8*ArcSin[Sqrt[3]*Sqrt[s/a]/4]/3]}

    deltoidunitspeed[a_][smin_,smax_,opt___]:=
        pplot[deltoidunitspeed[a][s],{s,smin,smax},opt]

    deltoidunitspeed[{a_,theta_}][s_]:=
    {a*(Cos[theta - (8*ArcSin[Sqrt[3]*Sqrt[s/a]/4])/3] + 
      2*Cos[theta + (4*ArcSin[Sqrt[3]*Sqrt[s/a]/4])/3]),
     a*(Sin[theta - (8*ArcSin[Sqrt[3]*Sqrt[s/a]/4])/3] + 
      2*Sin[theta + (4*ArcSin[Sqrt[3]*Sqrt[s/a]/4])/3])}

    deltoidunitspeed[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[deltoidunitspeed[{a,theta}][s],{s,smin,smax},opt]

    diamond[n_,a_,b_][t_]:=
        {a*Sqrt[Cos[t]^2]^(n - 1)*Cos[t],
         b*Sqrt[Sin[t]^2]^(n - 1)*Sin[t]}

    diamond[n_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[diamond[n,a,b][t],{t,tmin,tmax},opt]

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

    diamond[{n_,a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[diamond[{n,a,b,theta}][t],{t,tmin,tmax},opt]

    diamond[a_][t_]:=
        {a*Cos[t]*Sqrt[Cos[t]^2],a*Sin[t]*Sqrt[Sin[t]^2]}

    diamond[a_][tmin_,tmax_,opt___]:=
        pplot[diamond[a][t],{t,tmin,tmax},opt]

    diamond[{a_,theta_}][t_]:=
    {a*(Cos[t]*Sqrt[Cos[t]^2]*Cos[theta] - 
      Sin[t]*Sqrt[Sin[t]^2]*Sin[theta]),
     a*(Cos[theta]*Sin[t]*Sqrt[Sin[t]^2] + 
      Cos[t]*Sqrt[Cos[t]^2]*Sin[theta])}

    diamond[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[diamond[{a,theta}][t],{t,tmin,tmax},opt]

    eight[a_,b_][t_]:= {a*Sin[t],b*Sin[t]*Cos[t]}

    eight[a_,b_][tmin_,tmax_,opt___]:=
        pplot[eight[a,b][t],{t,tmin,tmax},opt]

    eight[{a_,b_,theta_}][t_]:=
    {a*Sin[t]*(Cos[theta] - Cos[t]*Sin[theta]),
     b*Sin[t]*(Cos[t]*Cos[theta] + Sin[theta])}

    eight[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[eight[{a,b,theta}][t],{t,tmin,tmax},opt]

    elasticainflect[a_,k_][s_]:= 
       {(2/a)*EllipticE[JacobiAmplitude[a*s,k^2],k^2] - s,
          -(2*k/a)*(JacobiCN[a*s,k^2] - 1)}

    elasticainflect/: Derivative[m_?IntegerQ][elasticainflect[a_,k_]] = 
        Derivative[m - 1][{-1 + 2*JacobiDN[a*#,k^2]^2,
                  2*k*JacobiDN[a*#,k^2]*JacobiSN[a*#,k^2]}&]

    elasticainflect[a_,k_][smin_,smax_,opt___]:=
        pplot[elasticainflect[a,k][s],
            {s,smin,smax},opt,PlotPoints->300]

    elasticainflect[{a_,k_,theta_}][s_]:= 
    {Cos[theta]*(-s + (2*EllipticE[JacobiAmplitude[a*s,k^2],
     k^2])/a) + (2*k*(-1 + JacobiCN[a*s,k^2])*Sin[theta])/a,
     (2*k*Cos[theta] -  2*k*Cos[theta]*JacobiCN[a*s,k^2] - 
      a*s*Sin[theta] + 2*EllipticE[JacobiAmplitude[a*s,k^2],k^2]*
       Sin[theta])/a}

    elasticainflect[{a_,k_,theta_}][smin_,smax_,opt___]:=
        pplot[elasticainflect[{a,k,theta}][s],
            {s,smin,smax},opt,PlotPoints->300]

(*    elasticainflect[a_,k_][smin_,smax_,opt___]:=
        pplot[elasticainflectnds[smin,smax][a,k][s],
            {s,smin,smax},opt,PlotPoints->300]

    elasticainflectnds[smin_,smax_][a_,k_][s_]:=
         intrinsic[2*a*k*JacobiCN[a*#,k^2]&,0,{0,0,0},
         MaxSteps->2000,{smin,smax}][s]
*)

    elasticanoninflect[a_,k_][s_]:= 
        {(2/(k*a))*EllipticE[JacobiAmplitude[a*s/k,k^2],k^2] +
              (1 - 2/k^2)*s,
          -(2/(k*a))*(JacobiDN[a*s/k,k^2] - 1)}

    elasticanoninflect/: Derivative[m_?IntegerQ][elasticanoninflect[a_,k_]] = 
        Derivative[m - 1][{1 - 2*JacobiSN[a*#/k,k^2]^2,
              2*JacobiCN[a*#/k,k^2]*JacobiSN[a*#/k,k^2]}&]

    elasticanoninflect[a_,k_][smin_,smax_,opt___]:=
        pplot[elasticanoninflect[a,k][s],
            {s,smin,smax},opt,PlotPoints->300]

    elasticanoninflect[{a_,k_,theta_}][s_]:= 
        {(2/(k*a))*EllipticE[JacobiAmplitude[a*s/k,k^2],k^2] +
              (1 - 2/k^2)*s,
          -(2/(k*a))*(JacobiDN[a*s/k,k^2] - 1)}

    elasticanoninflect[{a_,k_,theta_}][smin_,smax_,opt___]:=
        pplot[elasticanoninflect[{a,k,theta}][s],
            {s,smin,smax},opt,PlotPoints->300]

(*    elasticanoninflectnds[smin_,smax_][a_,k_][s_]:=
         intrinsic[(2*a*JacobiDN[a*#/k,k^2])/k&,0,{0,0,0},
         MaxSteps->2000,{smin,smax}][s]

    elasticanoninflect[a_,k_][smin_,smax_,opt___]:=
        pplot[elasticanoninflectnds[smin,smax][a,k][s],
            {s,smin,smax},opt,PlotPoints->300]
*)
    ellipse[a_,b_][t_]:= {a*Cos[t],b*Sin[t]}

    ellipse[a_,b_][tmin_,tmax_,opt___]:=
        pplot[ellipse[a,b][t],{t,tmin,tmax},opt]

    ellipse[{a_,b_,theta_}][t_]:=
    {a*Cos[t]*Cos[theta] - b*Sin[t]*Sin[theta],
     b*Cos[theta]*Sin[t] + a*Cos[t]*Sin[theta]}

    ellipse[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[ellipse[{a,b,theta}][t],{t,tmin,tmax},opt]

    ellipsebis[a_,b_][t_]:= {a*Sech[t],b*Tanh[t]}

    ellipsebis[a_,b_][tmin_,tmax_,opt___]:=
        pplot[ellipsebis[a,b][t],{t,tmin,tmax},opt]

    ellipsebis[{a_,b_,theta_}][t_]:=
    {a*Cos[t]*Cos[theta] - b*Sin[t]*Sin[theta],
     b*Cos[theta]*Sin[t] + a*Cos[t]*Sin[theta]}

    ellipsebis[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[ellipsebis[{a,b,theta}][t],{t,tmin,tmax},opt]

    ellipseinc[a_,b_,c_,psi_][t_]:=
        Module[{k},k=1 - b^2/a^2;
           {a*JacobiSN[t,k] + (c/a)*(a*Cos[psi]*JacobiCN[t,k] +
                     b*Sin[psi]*JacobiSN[t,k])/JacobiDN[t,k],
            b*JacobiCN[t,k] + (c/a)*(a*Sin[psi]*JacobiCN[t,k] -
                   b*Cos[psi]*JacobiSN[t,k])/JacobiDN[t,k]}]

    ellipseinc[a_,b_,c_,psi_][tmin_,tmax_,opt___]:=
        pplot[ellipseinc[a,b,c,psi][t],{t,tmin,tmax},opt]

    ellipseinc[{a_,b_,c_,psi_,theta_}][t_]:=
    {(a*JacobiCN[t,1 - b^2/a^2]*(c*Cos[psi + theta] - 
         b*JacobiDN[t,1 - b^2/a^2]*Sin[theta]) + 
      JacobiSN[t,1 - b^2/a^2]*
       (a^2*Cos[theta]*JacobiDN[t,1 - b^2/a^2] + 
         b*c*Sin[psi + theta]))/
    (a*JacobiDN[t,1 - b^2/a^2]),
    (JacobiSN[t,1 - b^2/a^2]*
       (-b*c*Cos[psi + theta] + 
         a^2*JacobiDN[t,1 - b^2/a^2]*Sin[theta]) + 
      a*JacobiCN[t,1 - b^2/a^2]*
       (b*Cos[theta]*JacobiDN[t,1 - b^2/a^2] + 
         c*Sin[psi + theta]))/(a*JacobiDN[t,1 - b^2/a^2])}

    ellipseinc[{a_,b_,c_,psi_,theta_}][tmin_,tmax_,opt___]:=
        pplot[ellipseinc[{a,b,c,psi,theta}][t],{t,tmin,tmax},opt]

    epicycloid[a_,b_][t_]:=
            {(a + b)*Cos[t] - b*Cos[(a + b)*t/b],
             (a + b)*Sin[t] - b*Sin[(a + b)*t/b]}

    epicycloid[a_,b_][tmin_,tmax_,opt___]:=
        pplot[epicycloid[a,b][t],{t,tmin,tmax},opt]

    epicycloid[{a_,b_,theta_}][t_]:=
    {(a + b)*Cos[t + theta] - b*Cos[t + (a*t)/b + theta],
     (a + b)*Sin[t + theta] - b*Sin[t + (a*t)/b + theta]}

    epicycloid[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[epicycloid[{a,b,theta}][t],{t,tmin,tmax},opt]

    epicycloidunitspeed[a_,b_][s_]:=
    {(a + b)*Cos[4*b*ArcSin[Sqrt[a*s/(b*(8*a + 8*b))]]/a] - 
      b*Cos[4*(a + b)*ArcSin[Sqrt[a*s/(b*(8*a + 8*b))]]/a],
     (a + b)*Sin[4*b*ArcSin[Sqrt[a*s/(b*(8*a + 8*b))]]/a] -
      b*Sin[4*(a + b)*ArcSin[Sqrt[a*s/(b*(8*a + 8*b))]]/a]}

    epicycloidunitspeed[a_,b_][smin_,smax_,opt___]:=
        pplot[epicycloidunitspeed[a,b][s],{s,smin,smax},opt]

    epicycloidunitspeed[{a_,b_,theta_}][s_]:=
    {(a + b)*Cos[theta + (4*b*
          ArcSin[Sqrt[(a*s)/(8*a*b + 8*b^2)]])/a] - 
      b*Cos[theta + (4*(a + b)*
          ArcSin[Sqrt[(a*s)/(8*a*b + 8*b^2)]])/a],
     (a + b)*Sin[theta + (4*b*
          ArcSin[Sqrt[(a*s)/(8*a*b + 8*b^2)]])/a] - 
      b*Sin[theta + (4*(a + b)*
          ArcSin[Sqrt[(a*s)/(8*a*b + 8*b^2)]])/a]}

    epicycloidunitspeed[{a_,b_,theta_}][smin_,smax_,opt___]:=
        pplot[epicycloidunitspeed[{a,b,theta}][s],{s,smin,smax},opt]

    epispiral[n_,a_][t_]:=
    {a*Cos[t]*Csc[n*t],a*Csc[n*t]*Sin[t]}

    epispiral[n_,a_][tmin_,tmax_,opt___]:=
        pplot[epispiral[n,a][t],{t,tmin,tmax},opt,
        PlotRange->{{-3,3},{-3,3}}]

    epispiral[{n_,a_,theta_}][t_]:=
    {a*Cos[t + theta]*Csc[n*t],a*Csc[n*t]*Sin[t + theta]}

    epispiral[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[epispiral[{n,a,theta}][t],{t,tmin,tmax},opt,
        PlotRange->{{-3,3},{-3,3}}]

    epitrochoid[a_,b_,h_][t_]:=
            {(a + b)*Cos[t] - h*Cos[(a + b)*t/b],
             (a + b)*Sin[t] - h*Sin[(a + b)*t/b]}

    epitrochoid[a_,b_,h_][tmin_,tmax_,opt___]:=
        pplot[epitrochoid[a,b,h][t],{t,tmin,tmax},opt]

    epitrochoid[{a_,b_,h_,theta_}][t_]:=
    {(a + b)*Cos[t + theta] - h*Cos[t + (a*t)/b + theta],
     (a + b)*Sin[t + theta] - h*Sin[t + (a*t)/b + theta]}

    epitrochoid[{a_,b_,h_,theta_}][tmin_,tmax_,opt___]:=
        pplot[epitrochoid[{a,b,h,theta}][t],{t,tmin,tmax},opt]

    fermatspiral[a_][t_]:= {a*Sqrt[t]*Cos[t],a*Sqrt[t]*Sin[t]}

    fermatspiral[a_][tmin_,tmax_,opt___]:=
        pplot[fermatspiral[a][t],{t,tmin,tmax},opt]

    fermatspiral[{a_,theta_}][t_]:=
    {a*Sqrt[t]*Cos[t + theta],a*Sqrt[t]*Sin[t + theta]}

    fermatspiral[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[fermatspiral[{a,theta}][t],{t,tmin,tmax},opt]

    fermatspiral[n_,a_][t_]:= {a*t/(t^2)^n*Cos[Sqrt[t^2]],
                               a*t/(t^2)^n*Sin[Sqrt[t^2]]}

    fermatspiral[{n_,a_,theta_}][t_]:=
    {(a*t*Cos[Sqrt[t^2] + theta])/(t^2)^n,
     (a*t*Sin[Sqrt[t^2] + theta])/(t^2)^n}

    fermatspiral[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[fermatspiral[{n,a,theta}][t],{t,tmin,tmax},opt]

    folium[t_]:= {3*t/(t^3 + 1),3*t^2/(t^3 + 1)}

    folium[tmin_,tmax_,opt___]:=
        pplot[folium[t],{t,tmin,tmax},opt]

    freeth[n_,a_][t_]:=
    {a*Cos[t]*(1 + n*Sin[t/n]),a*Sin[t]*(1 +n*Sin[t/n])}
    
    freeth[n_,a_][tmin_,tmax_,opt___]:=
        pplot[freeth[n,a][t],{t,tmin,tmax},opt]
        
    freeth[{n_,a_,theta_}][t_]:=
    {a*Cos[t + theta]*(1 + n*Sin[t/n]), 
     a*Sin[t + theta]*(1 + n*Sin[t/n])}
       
    freeth[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[freeth[{a,theta}][t],{t,tmin,tmax},opt]
            genfolium[n_,a_,b_,c_][t_]:= {c*t^n/(a + b*t^(2*n + 1)),
                               c*t^(n + 1)/(a + b*t^(2*n + 1))}

    genepicycloid[coef_List][t_]:=
        ComplexExpand[{Re[ExpToTrig[
        Sum[coef[[k]]*Exp[I*k*t],{k,1,Length[coef]}]]],
        Im[ExpToTrig[
        Sum[coef[[k]]*Exp[I*k*t],{k,1,Length[coef]}]]]}]

    genepicycloid[{coef_List,theta}][t_]:=
        ComplexExpand[{Re[ExpToTrig[
        Sum[coef[[k]]*Exp[I*k*(t + theta)],{k,1,Length[coef]}]]],
        Im[ExpToTrig[
        Sum[coef[[k]]*Exp[I*k*(t + theta)],{k,1,Length[coef]}]]]}]

    genfolium[n_,a_,b_,c_][tmin_,tmax_,opt___]:=
        pplot[genfolium[n,a,b,c][t],{t,tmin,tmax},opt]

    genfolium[n_,a_,b_,c_][tmin_,tmax_,opt___]:=
        pplot[genfolium[n,a,b,c][t],{t,tmin,tmax},opt]

    genfolium[{n_,a_,b_,c_,theta_}][t_]:=
    {(c*t^n*(Cos[theta] - t*Sin[theta]))/(a + b*t^(2*n + 1)),
     (c*t^n*(t*Cos[theta] + Sin[theta]))/(a + b*t^(2*n + 1))}

    genfolium[{n_,a_,b_,c_,theta_}][tmin_,tmax_,opt___]:=
        pplot[genfolium[{n,a,b,c,theta}][t],{t,tmin,tmax},opt]

    genparabola[n_][t_]:= {t,t^n}

    genparabola[n_][tmin_,tmax_,opt___]:=
        pplot[genparabola[n][t],
            {t,tmin,tmax},opt,AspectRatio->1]

    genparabola[{n_,theta_}][t_]:= 
    {t*Cos[theta] - t^n*Sin[theta],
     t^n*Cos[theta] + t*Sin[theta]}

    genparabola[{n_,theta_}][tmin_,tmax_,opt___]:=
        pplot[genparabola[{n,theta}][t],
            {t,tmin,tmax},opt,AspectRatio->1]

    hankelspiral[n_,a_][t_]:= {a*BesselJ[n,t],a*BesselY[n,t]}

    hankelspiral[n_,a_][tmin_,tmax_,opt___]:=
        pplot[hankelspiral[n,a][t],{t,tmin,tmax},opt]

    hankelspiral[{n_,a_,theta_}][t_]:=
    {a*(BesselJ[n,t]*Cos[theta] - BesselY[n,t]*Sin[theta]),
     a*(BesselY[n,t]*Cos[theta] + BesselJ[n,t]*Sin[theta])}

    hankelspiral[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hankelspiral[{n,a,theta}][t],{t,tmin,tmax},opt]

    hippias[a_,b_][t_]:= {a*t,a*t*Cot[Pi*t/(2*b)]}

    hippias[a_,b_][tmin_,tmax_,opt___]:=
        pplot[hippias[a,b][t],{t,tmin,tmax},opt,
        PlotRange->{Automatic,{-4,4}}]

    hippias[{a_,b_,theta_}][t_]:=
    {a*t*(Cos[theta] - Cot[(Pi*t)/(2*b)]*Sin[theta]),
     a*t*(Cos[theta]*Cot[(Pi*t)/(2*b)] + Sin[theta])}

    hippias[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hippias[{a,b,theta}][t],{t,tmin,tmax},opt,
        PlotRange->{Automatic,{-4,4}}]

    hippopede[a_,b_][t_]:=
              {2*Cos[t]*Sqrt[a*b - b^2*Sin[t]^2],
               2*Sin[t]*Sqrt[a*b - b^2*Sin[t]^2]}

    hippopede[a_,b_][tmin_,tmax_,opt___]:=
        pplot[hippopede[a,b][t],{t,tmin,tmax},opt]

    hippopede[{a_,b_,theta_}][t_]:=
    {Sqrt[2]*Sqrt[b*(2*a - b + b*Cos[2*t])]*Cos[t + theta],
     Sqrt[2]*Sqrt[b*(2*a - b + b*Cos[2*t])]*Sin[t + theta]}

    hippopede[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hippopede[{a,b,theta}][t],{t,tmin,tmax},opt]

    hyperbola[a_,b_][t_]:= {a*Cosh[t],b*Sinh[t]}

    hyperbola[a_,b_][tmin_,tmax_,opt___]:=
        pplot[hyperbola[a,b][t],{t,tmin,tmax},opt]

    hyperbola[{a_,b_,theta_}][t_]:=
    {a*Cos[theta]*Cosh[t] - b*Sin[theta]*Sinh[t],
     a*Sin[theta]*Cosh[t] + b*Cos[theta]*Sinh[t]}

    hyperbola[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hyperbola[{a,b,theta}][t],{t,tmin,tmax},opt]

    hyperbolabis[a_,b_][t_]:= {a*Sec[t],b*Tan[t]}

    hyperbolabis[a_,b_][tmin_,tmax_,opt___]:=
        pplot[hyperbolabis[a,b][t],{t,tmin,tmax},opt]

    hyperbolabis[{a_,b_,theta_}][t_]:=
    {a*Cos[theta]*Sec[t] - b*Sin[theta]*Tan[t],
     a*Sin[theta]*Sec[t] + b*Cos[theta]*Tan[t]}

    hyperbolabis[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hyperbolabis[{a,b,theta}][t],{t,tmin,tmax},opt]

    hyperbolicspiral[a_][t_]:= {a*Cos[t]/t,a*Sin[t]/t}

    hyperbolicspiral[a_][tmin_,tmax_,opt___]:=
        pplot[hyperbolicspiral[a][t],{t,tmin,tmax},opt]

    hyperbolicspiral[{a_,theta_}][t_]:=
    {a*Cos[t + theta]/t,a*Sin[t + theta]/t}

    hyperbolicspiral[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hyperbolicspiral[{a,theta}][t],{t,tmin,tmax},opt]

    hypocycloid[a_,b_][t_]:=
            {(a - b)*Cos[t] + b*Cos[(a - b)*t/b],
             (a - b)*Sin[t] - b*Sin[(a - b)*t/b]}
    
    hypocycloid[a_,b_][tmin_,tmax_,opt___]:=
        pplot[hypocycloid[a,b][t],{t,tmin,tmax},opt]

    hypocycloid[{a_,b_,theta_}][t_]:=
    {(a - b)*Cos[t + theta] + b*Cos[t - (a*t)/b + theta],
     (a - b)*Sin[t + theta] + b*Sin[t - (a*t)/b + theta]}
    
    hypocycloid[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hypocycloid[{a,b,theta}][t],{t,tmin,tmax},opt]

    hypocycloidinvolute[a_,b_][t_]:=
        {(a^2*Cos[t] - 3*a*b*Cos[t] + 2*b^2*Cos[t] - 
          a*b*Cos[t - a*t/b] + 2*b^2*Cos[t - a*t/b] + 
          4*a*b*Cos[t - a*t/(2*b)] - 4*b^2*Cos[t - a*t/(2*b)])/a,
         (a^2*Sin[t] - 3*a*b*Sin[t] + 2*b^2*Sin[t] - 
          a*b*Sin[t - a*t/b] + 2*b^2*Sin[t - a*t/b] + 
          4*a*b*Sin[t - a*t/(2*b)] - 4*b^2*Sin[t - a*t/(2*b)])/a}

    hypocycloidinvolute[a_,b_][tmin_,tmax_,opt___]:=
        pplot[hypocycloidinvolute[a,b][t],{t,tmin,tmax},opt]

    hypocycloidinvolute[{a_,b_,theta_}][t_]:=
    {((a^2 - 3*a*b + 2*b^2)*Cos[t + theta] + 
       b*(-((a - 2*b)*Cos[t - (a*t)/b + theta]) + 
         4*(a - b)*Cos[t - (a*t)/(2*b) + theta]))/a,
     ((a^2 - 3*a*b + 2*b^2)*Sin[t + theta] + 
       b*(-((a - 2*b)*Sin[t - (a*t)/b + theta]) + 
         4*(a - b)*Sin[t - (a*t)/(2*b) + theta]))/a}

    hypocycloidinvolute[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hypocycloidinvolute[{a,b,theta}][t],{t,tmin,tmax},opt]

    hypocycloidunitspeed[a_,b_][s_]:=
    {(a - b)*Cos[4*b*ArcSin[Sqrt[a*s/((8*a - 8*b)*b)]]/a] +
      b*Cos[4*(a - b)*ArcSin[Sqrt[a*s/((8*a - 8*b)*b)]]/a],
     (a - b)*Sin[4*b*ArcSin[Sqrt[a*s/((8*a - 8*b)*b)]]/a]
      -b*Sin[4*(a - b)*ArcSin[Sqrt[a*s/((8*a - 8*b)*b)]]/a]}

    hypocycloidunitspeed[a_,b_][smin_,smax_,opt___]:=
        pplot[hypocycloidunitspeed[a,b][s],{s,smin,smax},opt]

    hypocycloidunitspeed[{a_,b_,theta_}][s_]:=
    {(a - b)*Cos[theta + (4*b*ArcSin[Sqrt[(a*s)/(8*a*b - 8*b^2)]])/a] + 
      b*Cos[theta + (-4 + (4*b)/a)*ArcSin[Sqrt[(a*s)/(8*a*b - 8*b^2)]]],
     (a - b)*Sin[theta + (4*b*ArcSin[Sqrt[(a*s)/(8*a*b - 8*b^2)]])/a] + 
      b*Sin[theta + (-4 + (4*b)/a)*ArcSin[Sqrt[(a*s)/(8*a*b - 8*b^2)]]]}

    hypocycloidunitspeed[{a_,b_,theta_}][smin_,smax_,opt___]:=
        pplot[hypocycloidunitspeed[{a,b,theta}][s],{s,smin,smax},opt]

    hypotrochoid[a_,b_,h_][t_]:=
            {(a - b)*Cos[t] + h*Cos[(a - b)*t/b],
             (a - b)*Sin[t] - h*Sin[(a - b)*t/b]}

    hypotrochoid[a_,b_,h_][tmin_,tmax_,opt___]:=
        pplot[hypotrochoid[a,b,h][t],{t,tmin,tmax},opt]

    hypotrochoid[{a_,b_,h_,theta_}][t_]:=
    {(a - b)*Cos[t + theta] + h*Cos[t - (a*t)/b + theta],
     (a - b)*Sin[t + theta] + h*Sin[t - (a*t)/b + theta]}

    hypotrochoid[{a_,b_,h_,theta_}][tmin_,tmax_,opt___]:=
        pplot[hypotrochoid[{a,b,h,theta}][t],{t,tmin,tmax},opt]

    kampyle[a_][t_]:= {a*Sec[t],a*Sin[t]*Sec[t]^2}

    kampyle[a_][tmin_,tmax_,opt___]:=
        pplot[kampyle[a][t],{t,tmin,tmax},opt]

    kampyle[{a_,theta_}][t_]:=
    {a*Cos[t + theta]*Sec[t]^2,a*Sec[t]^2*Sin[t + theta]}

    kampyle[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[kampyle[{a,theta}][t],{t,tmin,tmax},opt]

    kappacurve[a_][t_]:= {a*Cos[t]^2/Sin[t],a*Cos[t]}

    kappacurve[a_][tmin_,tmax_,opt___]:=
        pplot[kappacurve[a][t],
            {t,tmin,tmax},opt,PlotRange->{{-2,2},Automatic}]

    kappacurve[{a_,theta_}][t_]:=
    {a*Cos[t + theta]*Cot[t],a*Sin[t + theta]*Cot[t]}

    kappacurve[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[kappacurve[{a,theta}][t],
            {t,tmin,tmax},opt,PlotRange->{{-2,2},Automatic}]

    keplerorbit[a_,e_][t_]:= {a*Cos[t]/(1 - e*Cos[t]),
                              a*Sin[t]/(1 - e*Cos[t])}

    keplerorbit[a_,e_][tmin_,tmax_,opt___]:=
        pplot[keplerorbit[a,e][t],{t,tmin,tmax},opt,Axes->True]

    keplerorbit[{a_,e_,theta_}][t_]:=
    {a*Cos[t + theta]/(1 - e*Cos[t]),
     a*Sin[t + theta]/(1 - e*Cos[t])}

    keplerorbit[{a_,e_,theta_}][tmin_,tmax_,opt___]:=
        pplot[keplerorbit[{a,e,theta}][t],{t,tmin,tmax},opt,Axes->True]

    lehr[a_,b_,c_][smin_,smax_,opt___]:=
        pplot[lehrnds[smin,smax][a,b,c][s],
            {s,smin,smax},opt,PlotPoints->300]

    lehr[a_,coef_List,c_][smin_,smax_,opt___]:=
        pplot[lehrnds[smin,smax][a,coef,c][s],
            {s,smin,smax},opt,PlotPoints->300]

    lehr/: Derivative[m_?IntegerQ][lehr[a_,b_,c_]] = 
        Derivative[m - 1][{Cos[a*# + (b/c)*Sin[c*#]],
        Sin[a*# + (b/c)*Sin[c*#]]}&]

    lehrkappa2[a_,b_,c_][s_]:= a + b*Cos[c*s]

    lehrkappa2[a_,coef_List,c_][s_]:= 
         a + coef.Table[Cos[k*c*s],{k,1,Length[coef]}]
    
    lehrnds[smin_,smax_][a_,b_,c_][s_]:=
         intrinsic[lehrkappa2[a,b,c],0,{0,0,0},
         MaxSteps->2000,{smin,smax}][s]

    lehrnds[smin_,smax_][a_,coef_List,c_][s_]:=
         intrinsic[lehrkappa2[a,coef,c],0,{0,0,0},
         MaxSteps->2000,{smin,smax}][s]

    lemniscate[a_][t_]:=
         {a*Cos[t]/(1 + Sin[t]^2),
          a*Sin[t]*Cos[t]/(1 + Sin[t]^2)}

    lemniscate[a_][tmin_,tmax_,opt___]:=
        pplot[lemniscate[a][t],{t,tmin,tmax},opt]

    lemniscate[{a_,theta_}][t_]:=
    {a*Cos[t]*(Cos[theta] - Sin[t]*Sin[theta])/(1 + Sin[t]^2),
     a*Cos[t]*(Cos[theta]*Sin[t] + Sin[theta])/(1 + Sin[t]^2)}

    lemniscate[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[lemniscate[{a,theta}][t],{t,tmin,tmax},opt]

    lemniscatebis[a_][t_]:=
         {a*Sqrt[Cos[2*t]]*Cos[t],a*Sqrt[Cos[2*t]]*Sin[t]}

    lemniscatebis[a_][tmin_,tmax_,opt___]:=
        pplot[lemniscatebis[a][t],{t,tmin,tmax},opt]

    lemniscatebis[{a_,theta_}][t_]:=
    {a*Sqrt[Cos[2*t]]*Cos[t + theta],a*Sqrt[Cos[2*t]]*Sin[t + theta]}

    lemniscatebis[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[lemniscatebis[{a,theta}][t],{t,tmin,tmax},opt]

    lemniscateunitspeed[a_][s_]:=
         {a*JacobiCN[s/a,-1]/JacobiDN[s/a,-1]^2,
          a*JacobiCN[s/a,-1]*JacobiSN[s/a,-1]/JacobiDN[s/a,-1]^2}

    lemniscateunitspeed[a_][smin_,smax_,opt___]:=
        pplot[lemniscateunitspeed[a][s],{s,smin,smax},opt]

    lemniscateunitspeed[{a_,theta_}][s_]:=
    {a*JacobiCN[s/a,-1]*(Cos[theta] - 
        JacobiSN[s/a,-1]*Sin[theta])/JacobiDN[s/a,-1]^2,
     a*JacobiCN[s/a,-1]*
      (Cos[theta]*JacobiSN[s/a,-1] + Sin[theta])/JacobiDN[s/a,-1]^2}

    lemniscateunitspeed[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[lemniscateunitspeed[{a,theta}][s],{s,smin,smax},opt]

    limacon[a_,b_][t_]:= {(2*a*Cos[t] + b)*Cos[t],
                          (2*a*Cos[t] + b)*Sin[t]}

    limacon[a_,b_][tmin_,tmax_,opt___]:=
        pplot[limacon[a,b][t],{t,tmin,tmax},opt]

    limacon[{a_,b_,theta_}][t_]:=
    {(b + 2*a*Cos[t])*Cos[t + theta],
     (b + 2*a*Cos[t])*Sin[t + theta]}

    limacon[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[limacon[{a,b,theta}][t],{t,tmin,tmax},opt]

    line[a1_,b1_,a2_,b2_][t_]:= {a1*t + a2,b1*t + b2}

    line[a1_,b1_,a2_,b2_][tmin_,tmax_,opt___]:=
        pplot[line[a1,b1,a2,b2][t],{t,tmin,tmax},opt]

    line[{a1_,b1_,a2_,b2_,theta_}][t_]:=
    {(a2 + a1*t)*Cos[theta] - (b2 + b1*t)*Sin[theta],
     (b2 + b1*t)*Cos[theta] + (a2 + a1*t)*Sin[theta]}

    line[{a1_,b1_,a2_,b2_,theta_}][tmin_,tmax_,opt___]:=
        pplot[line[{a1,b1,a2,b2,theta}][t],{t,tmin,tmax},opt]

    linearpursuit[a_,k_][t_]:= {t,a*k/(k^2 - 1) +
        k*(a - t)^(1 + 1/k)/(2*a^(1/k)*(k + 1)) -
        k*(a - t)^(1 - 1/k)*a^(1/k)/(2*(k - 1))}

    linearpursuit[a_,1][t_]:= {t,(a/4)*(((a - t)/a)^2 -
        1 -2*Log[(a - t)/a])}

    linearpursuit[a_,k_][tmin_,tmax_,opt___]:=
        pplot[linearpursuit[a,k][t],{t,tmin,tmax},opt,Axes->True]

    linearpursuit[{a_,k_,theta_}][t_]:=
    {t*Cos[theta] + (k*((-2*a)/(-1 + k^2) - 
         (a - t)^(1 + k^(-1))/(a^k^(-1)*(1 + k)) + 
         (a^k^(-1)*(a - t)^((-1 + k)/k))/(-1 + k))*Sin[theta])/2,
     (k*((2*a)/(-1 + k^2) +(a - t)^(1 + k^(-1))/(a^k^(-1)*(1 + k)) - 
         (a^k^(-1)*(a - t)^((-1 + k)/k))/(-1 + k))*
       Cos[theta])/2 + t*Sin[theta]}

    linearpursuit[{a_,1,theta_}][t_]:=
    {t*Cos[theta] + (((2*a - t)*t + 2*a^2*Log[1 - t/a])*Sin[theta])/(4*a),
    (a*Cos[theta]*((t*(-2*a + t))/a^2 - 2*Log[1 - t/a]))/4 + t*Sin[theta]}

    linearpursuit[{a_,k_,theta_}][tmin_,tmax_,opt___]:=
        pplot[linearpursuit[{a,k,theta}][t],{t,tmin,tmax},opt,Axes->True]

    lissajous[n_,d_,a_,b_][t_]:= {a*Sin[n*t + d],b*Sin[t]}

    lissajous[n_,d_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[lissajous[n,d,a,b][t],{t,tmin,tmax},opt]

    lissajous[{n_,d_,a_,b_,theta_}][t_]:=
    {a*Cos[theta]*Sin[d + n*t] - b*Sin[t]*Sin[theta],
     b*Cos[theta]*Sin[t] + a*Sin[d + n*t]*Sin[theta]}

    lissajous[{n_,d_,a_,b_,theta_}n_,d_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[lissajous[{n,d,a,b,theta}][t],{t,tmin,tmax},opt]

    lituus[a_][t_]:= {a*t/(t^2)^(3/4)*Cos[Sqrt[t^2]],
                      a*t/(t^2)^(3/4)*Sin[Sqrt[t^2]]}

    lituus[a_][tmin_,tmax_,opt___]:=
        pplot[lituus[a][t],{t,tmin,tmax},opt]

    lituus[{a_,theta_}][t_]:=
    {a*(t^2)^(1/4)*Cos[Sqrt[t^2] + theta]/t,
     a*(t^2)^(1/4)*Sin[Sqrt[t^2] + theta]/t}

    lituus[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[lituus[{a,theta}][t],{t,tmin,tmax},opt]

    logistic[k_,a_,b_][t_]:= {t,k/(1 + E^(a + b*t))}

    logistic[k_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[logistic[k,a,b][t],{t,tmin,tmax},opt]

    logistic[{k_,a_,b_,theta_}][t_]:=
    {a*(t^2)^(1/4)*Cos[Sqrt[t^2] + theta]/t,
     a*(t^2)^(1/4)*Sin[Sqrt[t^2] + theta]/t}

    logistic[{k_,a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[logistic[{k,a,b,theta}][t],{t,tmin,tmax},opt]

    logspiral[b_][t_]:= {E^(b*t)*Cos[t],E^(b*t)*Sin[t]}

    logspiral[b_][tmin_,tmax_,opt___]:=
        pplot[logspiral[b][t],{t,tmin,tmax},opt]

    logspiral[{b_,theta_}][t_]:=
    {E^(b*t)*Cos[t + theta],E^(b*t)*Sin[t + theta]}

    logspiral[{b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[logspiral[{b,theta}][t],{t,tmin,tmax},opt]

    logspiral[a_,b_][t_]:= {a*E^(b*t)*Cos[t],a*E^(b*t)*Sin[t]}

    logspiral[a_,b_][tmin_,tmax_,opt___]:=
        pplot[logspiral[a,b][t],{t,tmin,tmax},opt]

    logspiral[{a_,b_,theta_}][t_]:=
    {a*E^(b*t)*Cos[t + theta],a*E^(b*t)*Sin[t + theta]}

    logspiral[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[logspiral[{a,b,theta}][t],{t,tmin,tmax},opt]

    logspiralunitspeed[a_,b_][s_]:= 
         {b*s*Cos[Log[b*s/(a*Sqrt[1 + b^2])]/b]/Sqrt[1 + b^2],
          b*s*Sin[Log[b*s/(a*Sqrt[1 + b^2])]/b]/Sqrt[1 + b^2]}

    logspiralunitspeed[a_,b_][smin_,smax_,opt___]:=
        pplot[logspiralunitspeed[a,b][s],{s,smin,smax},opt]

    logspiralunitspeed[{a_,b_,theta_}][s_]:= 
    {b*s*Cos[theta + Log[(b*s)/(a*Sqrt[1 + b^2])]/b]/Sqrt[1 + b^2],
     b*s*Sin[theta + Log[(b*s)/(a*Sqrt[1 + b^2])]/b]/Sqrt[1 + b^2]}

    logspiralunitspeed[{a_,b_,theta_}][smin_,smax_,opt___]:=
        pplot[logspiralunitspeed[{a,b,theta}][s],{s,smin,smax},opt]

    nephroid[a_][t_]:= {a*(3*Cos[t] - Cos[3*t]),
                        a*(3*Sin[t] - Sin[3*t])}

    nephroid[a_][tmin_,tmax_,opt___]:=
        pplot[nephroid[a][t],{t,tmin,tmax},opt]

    nephroid[{a_,theta_}][t_]:=
    {a*(3*Cos[t + theta] - Cos[3*t + theta]),
     a*(3*Sin[t + theta] - Sin[3*t + theta])}

    nephroid[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[nephroid[{a,theta}][t],{t,tmin,tmax},opt]

    nephroidunitspeed[a_][s_]:=
      {s/2 - a*Cos[3*ArcCos[s/(6*a)]],
       a*(3*Sqrt[1 - s^2/(36*a^2)] - Sin[3*ArcCos[s/(6*a)]])}

    nephroidunitspeed[a_][smin_,smax_,opt___]:=
        pplot[nephroidunitspeed[a][s],{s,smin,smax},opt]
        
    nephroidunitspeed[{a_,theta_}][s_]:=
    {(s*Cos[theta] - a*(2*Cos[theta + 3*ArcCos[s/(6*a)]] + 
         Sqrt[36 - s^2/a^2]*Sin[theta]))/2,
     (a*Sqrt[36 - s^2/a^2]*Cos[theta] + s*Sin[theta] - 
      2*a*Sin[theta + 3*ArcCos[s/(6*a)]])/2}

    nephroidunitspeed[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[nephroidunitspeed[{a,theta}][s],{s,smin,smax},opt]
        
    ngon[n_,a_][t_]:=
        Which@@Flatten[Table[{(k - 1)/n<=t<=k/n,
        a*{Cos[2*Pi*t],Sin[2*Pi*t]}/Cos[2*Pi*t - (2*k - 1)*Pi/n]},
        {k,1,n}],1]

    ngon[n_,a_][tmin_,tmax_,opt___]:=
        ParametricPlot[Evaluate[ngon[n,a][t]],
        {t,tmin,tmax},opt,Axes->None,AspectRatio->Automatic,
         PlotStyle->{{RGBColor[1,0,0],AbsoluteThickness[2]}}]

    nielsenspiral[a_][t_]:= {a*CosIntegral[t],a*SinIntegral[t]}

    nielsenspiral[a_][tmin_,tmax_,opt___]:=
        pplot[nielsenspiral[a][t],
            {t,tmin,tmax},opt,PlotRange->All]

    nielsenspiral[{a_,theta_}][t_]:=
    {a*(Cos[theta]*CosIntegral[t] - Sin[theta]*SinIntegral[t]),
     a*(CosIntegral[t]*Sin[theta] + Cos[theta]*SinIntegral[t])}

    nielsenspiral[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[nielsenspiral[{a,theta}][t],
            {t,tmin,tmax},opt,PlotRange->All]

    parabola[a_][t_]:= {2*a*t,a*t^2}

    parabola[a_][tmin_,tmax_,opt___]:=
        pplot[parabola[a][t],{t,tmin,tmax},opt]

    parabola[{a_,theta_}][t_]:=
    {a*t*(2*Cos[theta] - t*Sin[theta]),
     a*t*(t*Cos[theta] + 2*Sin[theta])}

    parabola[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[parabola[{a,theta}][t],{t,tmin,tmax},opt]

    parabola[a_,b_][t_]:= {2*b*t,a*t^2}

    parabola[a_,b_][tmin_,tmax_,opt___]:=
        pplot[parabola[a,b][t],{t,tmin,tmax},opt]

    parabola[{a_,b_,theta_}][t_]:=
    {t*(2*b*Cos[theta] - a*t*Sin[theta]),
     t*(a*t*Cos[theta] + 2*b*Sin[theta])}

    parabola[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[parabola[{a,b,theta}][t],{t,tmin,tmax},opt]

    parabolicspiral[n_,a_,b_][t_]:=
        {(a*t/(t^2)^n + b)*Cos[Sqrt[t^2]],
         (a*t/(t^2)^n + b)*Sin[Sqrt[t^2]]}

    parabolicspiral[n_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[parabolicspiral[n,a,b][t],{t,tmin,tmax},opt]

    parabolicspiral[{n_,a_,b_,theta_}][t_]:=
    {(a*t + b*(t^2)^n)*Cos[Sqrt[t^2] + theta]/(t^2)^n,
     (a*t + b*(t^2)^n)*Sin[Sqrt[t^2] + theta]/(t^2)^n}

    parabolicspiral[{n_,a_,b_,theta}][tmin_,tmax_,opt___]:=
        pplot[parabolicspiral[{n,a,b,theta}][t],{t,tmin,tmax},opt]

    piriform[a_,b_][t_]:=
           {a*(1 + Sin[t]),b*Cos[t]*(1 + Sin[t])}

    piriform[a_,b_][tmin_,tmax_,opt___]:=
        pplot[piriform[a,b][t],{t,tmin,tmax},opt]

    piriform[{a_,b_,theta_}][t_]:=
    {(1 + Sin[t])*(a*Cos[theta] - b*Cos[t]*Sin[theta]),
     (1 + Sin[t])*(b*Cos[t]*Cos[theta] + a*Sin[theta])}

    piriform[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[piriform[{a,b,theta}][t],{t,tmin,tmax},opt]

    pseudocatenary[a_,b_][tmin_,tmax_,opt___]:=
        pplot[pseudocatenarynds[tmin,tmax][a,b][t],
            {t,tmin,tmax},opt]

    pseudocatenary[a_,b_][s_]:= alysoid[a,b,a][s]

    pseudocatenary/: Derivative[m_?IntegerQ][pseudocatenary[a_,b_]] = 
        Derivative[m - 1][pseudocatenaryprime[a,b]]

    pseudocatenarykappa2[a_,b_][s_]:= -a/(a*b + s^2)

    pseudocatenarynds[smin_,smax_][a_,b_][s_]:=
         alysoidnds[smin,smax][a,b,a][s]

    pseudocatenaryprime[a_,b_][t_]:= alysoidprime[a,b,a][t]

    rectangle[n_,a_,b_][t_]:=
            {a*(Cos[t]*Sqrt[Cos[t]^2]^(n - 1) +
                Sin[t]*Sqrt[Sin[t]^2]^(n - 1)),
             b*(-Cos[t]*Sqrt[Cos[t]^2]^(n - 1) +
                Sin[t]*Sqrt[Sin[t]^2]^(n - 1))}

    rectangle[n_,a_,b_][tmin_,tmax_,opt___]:=
        pplot[rectangle[n,a,b][t],{t,tmin,tmax},opt]

    rectangle[{n_,a_,b_,theta_}][t_]:=
    {Sin[t]*Sqrt[Sin[t]^2]^(n - 1)*(a*Cos[theta] - b*Sin[theta]) + 
     Cos[t]*Sqrt[Cos[t]^2]^(n - 1)*(a*Cos[theta] + b*Sin[theta]),
     -Cos[t]*Sqrt[Cos[t]^2]^(n - 1)*(b*Cos[theta] - a*Sin[theta]) + 
       Sin[t]*Sqrt[Sin[t]^2]^(n - 1)*(b*Cos[theta] + a*Sin[theta])}

    rectangle[{n_,a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[rectangle[{n,a,b,theta}][t],{t,tmin,tmax},opt]

    reuleauxpolygon[n_,a_][t_]:= Which@@Flatten[
        Table[{k/(2*n + 1)<=t<=(k + 1)/(2*n + 1),
            a*{Re[E^(2*Pi*I*k/(2*n + 1)) +
               2*I*Sin[Pi*n/(2*n + 1)]*
               E^(Pi*I*((k + n)/(2*n + 1) + t))],
               Im[E^(2*Pi*I*k/(2*n + 1)) +
               2*I*Sin[Pi*n/(2*n + 1)]*
               E^(Pi*I*((k + n)/(2*n + 1) + t))]}},
              {k,0,2*n}],1]

    reuleauxpolygon[n_,a_][tmin_,tmax_,opt___]:=
        ParametricPlot[Evaluate[reuleauxpolygon[n,a][t]],
        {t,tmin,tmax},opt,Axes->None,AspectRatio->Automatic,
         PlotStyle->{{RGBColor[1,0,0],AbsoluteThickness[2]}}]

    riemannfractal[k_,a_,m_][t_]:=
    Sum[{Cos[2*Pi*n^k*t],Sin[2*Pi*n^k*t]}/n^a,{n,1,m}]

    riemannfractal[k_,a_,m_][tmin_,tmax_,opt___]:=
        pplot[riemannfractal[k,a,m][t],{t,tmin,tmax},opt,
        PlotPoints->1000,PlotStyle->{}]

    rose[n_,a_][t_]:= {a*Cos[n*t]*Cos[t],a*Cos[n*t]*Sin[t]}

    rose[n_,a_][tmin_,tmax_,opt___]:=
        pplot[rose[n,a][t],{t,tmin,tmax},opt]

    rose[{n_,a_,theta_}][t_]:=
    {a*Cos[n*t]*Cos[t + theta],a*Cos[n*t]*Sin[t + theta]}

    rose[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[rose[{n,a,theta}][t],{t,tmin,tmax},opt]

    sc[a_][t_]:= {(a*t)^2,(a*t)^3}

    sc[a_][tmin_,tmax_,opt___]:=
        pplot[sc[a][t],{t,tmin,tmax},opt]

    sc[{a_,theta_}][t_]:=
    {(a*t)^2*Cos[theta] - (a*t)^3*Sin[theta],
     (a*t)^2*Sin[theta] + (a*t)^3*Cos[theta]}

    sc[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[sc[{a,theta}][t],{t,tmin,tmax},opt]

    scarab[a_,b_][t_]:= {(a*Cos[2*t] - b*Cos[t])*Cos[t],
                         (a*Cos[2*t] - b*Cos[t])*Sin[t]}

    scarab[a_,b_][tmin_,tmax_,opt___]:=
        pplot[scarab[a,b][t],{t,tmin,tmax},opt]

    scarab[{a_,b_,theta_}][t_]:=
    {(-b*Cos[t] + a*Cos[2*t])*Cos[t + theta],
     (-b*Cos[t] + a*Cos[2*t])*Sin[t + theta]}

    scarab[{a_,b_theta_}][tmin_,tmax_,opt___]:=
        pplot[scarab[{a,b,theta}][t],{t,tmin,tmax},opt]

    scunitspeed[s_]:= {(s + 8/27)^(2/3) - 4/9,
                           ((s + 8/27)^(2/3) - 4/9)^(3/2)}

    scunitspeed[smin_,smax_,opt___]:=
        pplot[scunitspeed[s],{s,smin,smax},opt]

    semicubic[a_,b_,c_,d_][z_]:=
     {WeierstrassP[z,{b^2/12 - a*c/4,-b^3/216 + a*b*c/48 - a^2*d/16}],
      WeierstrassPPrime[z,{b^2/12 - a*c/4,-b^3/216 + a*b*c/48 - a^2*d/16}]}

    semicubic[a_,b_,c_,d_][tmin_,tmax_,opt___]:=
        pplot[semicubic[a,b,c,d][t],{t,tmin,tmax},opt]

     semicubic[{a_,b_,c_,d_,theta_}][z_]:=
     {Cos[theta]*WeierstrassP[z,{(b^2 - 3*a*c)/12,
      (-2*b^3 + 9*a*b*c - 27*a^2*d)/432}] - 
      Sin[theta]*WeierstrassPPrime[z,{(b^2 - 3*a*c)/12,
       (-2*b^3 + 9*a*b*c - 27*a^2*d)/432}],
      Sin[theta]*WeierstrassP[z,{(b^2 - 3*a*c)/12,
       (-2*b^3 + 9*a*b*c - 27*a^2*d)/432}] + 
      Cos[theta]*WeierstrassPPrime[z,{(b^2 - 3*a*c)/12,
       (-2*b^3 + 9*a*b*c - 27*a^2*d)/432}]}

    semicubic[{a_,b_,c_,d_,theta_}][tmin_,tmax_,opt___]:=
        pplot[semicubic[{a,b,c,d,theta}][t],{t,tmin,tmax},opt]

    serpentine[a_,b_][t_]:= {t,a*t/(1 + b*t^2)}

    serpentine[a_,b_][tmin_,tmax_,opt___]:=
        pplot[serpentine[a,b][t],{t,tmin,tmax},opt]

    serpentine[{a_,b_,theta_}][t_]:=
    {t*(Cos[theta] - a*Sin[theta]/(1 + b*t^2)),
     t*(a*Cos[theta]/(1 + b*t^2) + Sin[theta])}

    serpentine[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[serpentine[{a,b,theta}][t],{t,tmin,tmax},opt]

    sinhspiral[n_,a_][t_]:=
            {a*Cos[t]/Sinh[n*t],a*Sin[t]/Sinh[n*t]}

    sinhspiral[n_,a_][tmin_,tmax_,opt___]:=
        pplot[sinhspiral[n,a][t],
            {t,tmin,tmax},opt,PlotRange->{{-1,1},Automatic}]

    sinhspiral[{n_,a_,theta_}][t_]:=
    {a*Cos[t + theta]*Csch[n*t],a*Sin[t + theta]*Csch[n*t]}

    sinhspiral[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[sinhspiral[{n,a,theta}][t],
            {t,tmin,tmax},opt,PlotRange->{{-1,1},Automatic}]

    sinoval[n_,a_][t_]:= {a*Cos[t],a*Nest[Sin,t,n]}

    sinoval[n_,a_][tmin_,tmax_,opt___]:=
        pplot[sinoval[n,a][t],{t,tmin,tmax},opt]

    spring[a_,b_,v_][t_]:=
        {a*Cos[t],a*Cos[v]*Sin[t] + b*t*Sin[v]}

    spring[a_,b_,v_][tmin_,tmax_,opt___]:=
        pplot[spring[a,b,v][t],{t,tmin,tmax},opt]

    spring[{a_,b_,v_}][t_]:=
    {a*Cos[t]*Cos[theta] - Sin[theta]*(a*Cos[v]*Sin[t] + b*t*Sin[v]),
     a*Cos[t]*Sin[theta] + Cos[theta]*(a*Cos[v]*Sin[t] + b*t*Sin[v])}

    spring[{a_,b_,v_}][tmin_,tmax_,opt___]:=
        pplot[spring[{a,b,v}][t],{t,tmin,tmax},opt]

    strophoid[a_][t_]:= {a*(t^2 - 1)/(t^2 + 1),
                         t*a*(t^2 - 1)/(t^2 + 1)}

    strophoid[a_][tmin_,tmax_,opt___]:=
        pplot[strophoid[a][t],{t,tmin,tmax},opt]

    strophoid[{a_,theta_}][t_]:=
    {a*(t^2 - 1)*(Cos[theta] - t*Sin[theta])/(t^2 + 1),
     a*(t^2 - 1)*(t*Cos[theta] + Sin[theta])/(t^2 + 1)}

    strophoid[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[strophoid[{a,theta}][t],{t,tmin,tmax},opt]

    strophoid[m_,a_][t_]:= {a*(t^2 - m)/(t^2 + 1),
                            t*a*(t^2 - m)/(t^2 + 1)}

    strophoid[m_,a_][tmin_,tmax_,opt___]:=
        pplot[strophoid[m,a][t],{t,tmin,tmax},opt]

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

    strophoid[{m_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[strophoid[{m,a,theta}][t],{t,tmin,tmax},opt]

    talbot[a_,b_,f_][t_]:= {(a^2 + f^2*Sin[t]^2)*Cos[t]/a,
                            (a^2 - 2*f^2 + f^2*Sin[t]^2)*Sin[t]/b}

    talbot[a_,b_,f_][tmin_,tmax_,opt___]:=
        pplot[talbot[a,b,f][t],{t,tmin,tmax},opt]

    talbot[{a_,b_,f_,theta_}][t_]:=
    {Cos[t]*Cos[theta]*(a^2 + f^2*Sin[t]^2)/a - 
     Sin[t]*(a^2 - 2*f^2 + f^2*Sin[t]^2)*Sin[theta]/b,
     Cos[theta]*Sin[t]*(a^2 - 2*f^2 + f^2*Sin[t]^2)/b + 
     Cos[t]*(a^2 + f^2*Sin[t]^2)*Sin[theta]/a}

    talbot[{a_,b_,f_,theta_}][tmin_,tmax_,opt___]:=
        pplot[talbot[{a,b,f,theta}][t],{t,tmin,tmax},opt]

    tanhspiral[a_][t_]:= {Sinh[2*t]/(Cos[2*a*t] + Cosh[2*t]),
                          Sin[2*a*t]/(Cos[2*a*t] + Cosh[2*t])}

    tanhspiral[a_][tmin_,tmax_,opt___]:=
        pplot[tanhspiral[a][t],{t,tmin,tmax},opt]

    tanhspiral[{a_,theta_}][t_]:=
    {-(Sin[2*a*t]*Sin[theta] + Cos[theta]*Sinh[2*t])/
     (Cos[2*a*t] +  Cosh[2*t]),
     (Cos[theta]*Sin[2*a*t] + Sin[theta]*Sinh[2*t])/
     (Cos[2*a*t] + Cosh[2*t])}

    tanhspiral[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[tanhspiral[{a,theta}][t],{t,tmin,tmax},opt]

    teardrop[a_,b_,c_,d_,n_][t_]:=
          {a*Sin[t],d*(b*Sin[t] + c*Sin[2*t])^n}

    teardrop[a_,b_,c_,d_,n_][tmin_,tmax_,opt___]:=
        pplot[teardrop[a,b,c,d,n][t],{t,tmin,tmax},opt]

    teardrop[{a_,b_,c_,d_,n_,theta_}][t_]:=
    {a*Cos[theta]*Sin[t] - d*((b + 2*c*Cos[t])*Sin[t])^n*Sin[theta],
     d*Cos[theta]*((b + 2*c*Cos[t])*Sin[t])^n + a*Sin[t]*Sin[theta]}

    teardrop[{a_,b_,c_,d_,n_,theta_}][tmin_,tmax_,opt___]:=
        pplot[teardrop[{a,b,c,d,n,theta}][t],{t,tmin,tmax},opt]

    tractrix[a_][t_]:= {a*Sin[t],a*(Cos[t] + Log[Tan[t/2]])}

    tractrix[a_][tmin_,tmax_,opt___]:=
        pplot[tractrix[a][t],{t,tmin,tmax},opt,
        PlotRange->{Automatic,{-1.5,1.5}},
        Ticks->{None,Automatic}]

    tractrix[{a_,theta_}][t_]:=
    {a*Cos[theta]*Sin[t] - a*Sin[theta]*(Cos[t] + Log[Tan[t/2]]),
     a*(Cos[t - theta] + Cos[theta]*Log[Tan[t/2]])}

    tractrix[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[tractrix[{a,theta}][t],{t,tmin,tmax},opt,
        PlotRange->{Automatic,{-1.5,1.5}},
        Ticks->{None,Automatic}]

    tractrix[a_,b_][t_]:= {a*Sin[t],b*(Cos[t] + Log[Tan[t/2]])}

    tractrix[a_,b_][tmin_,tmax_,opt___]:=
        pplot[tractrix[a,b][t],{t,tmin,tmax},opt,
        PlotRange->{Automatic,{-1.5,1.5}},
        Ticks->{None,Automatic}]

    tractrix[{a_,b_,theta_}][t_]:=
    {a*Cos[theta]*Sin[t] - b*(Cos[t] + Log[Tan[t/2]])*Sin[theta],
     b*Cos[theta]*(Cos[t] + Log[Tan[t/2]]) + a*Sin[t]*Sin[theta]}

    tractrix[{a_,b_,theta_}][tmin_,tmax_,opt___]:=
        pplot[tractrix[{a,b,theta}][t],{t,tmin,tmax},opt]

    tractrixminus[a_][t_]:= {t,Sqrt[a^2 - t^2] - a*Log[2/(a*t) +
         2*Sqrt[a^2 - t^2]/(a^2*t)] - a*Log[2/t^2]}

    tractrixminus[a_][tmin_,tmax_,opt___]:=
        pplot[tractrixminus[a][t],{t,tmin,tmax},opt]

    tractrixminus[{a_,theta_}][t_]:=
    {t*Cos[theta] + (-Sqrt[a^2 - t^2] + a*Log[2/t^2] + 
       a*Log[(2*(a + Sqrt[a^2 - t^2]))/(a^2*t)])*Sin[theta],
     Cos[theta]*(Sqrt[a^2 - t^2] - a*Log[2/t^2] - 
       a*Log[(2*(a + Sqrt[a^2 - t^2]))/(a^2*t)]) + t*Sin[theta]}

    tractrixminus[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[tractrixminus[{a,theta}][t],{t,tmin,tmax},opt]

    tractrixplus[a_][t_]:= {t,-Sqrt[a^2 - t^2] + a*Log[2/(a*t) +
         2*Sqrt[a^2 - t^2]/(a^2*t)] + a*Log[2/t^2]}

    tractrixplus[a_][tmin_,tmax_,opt___]:=
        pplot[tractrixplus[a][t],{t,tmin,tmax},opt]

    tractrixplus[{a_,theta_}][t_]:=
    {t*Cos[theta] + (Sqrt[a^2 - t^2] - a*Log[2/t^2] -
     a*Log[(2*(a + Sqrt[a^2 - t^2]))/(a^2*t)])*Sin[theta],
     Cos[theta]* (-Sqrt[a^2 - t^2] + a*Log[2/t^2] + 
     a*Log[(2*(a + Sqrt[a^2 - t^2]))/(a^2*t)]) + t*Sin[theta]}

    tractrixplus[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[tractrixplus[{a,theta}][t],{t,tmin,tmax},opt]

    tractrixunitspeed[a_][s_]:= {a*E^(-Sqrt[s^2]/a),
        a*ArcTanh[Sqrt[1 - E^(-2*Sqrt[s^2]/a)]] -
        a*Sqrt[1 - E^(-2*Sqrt[s^2]/a)]}

    tractrixunitspeed[a_][smin_,smax_,opt___]:=
        pplot[tractrixunitspeed[a][s],{s,smin,smax},opt]

    tractrixunitspeed[{a_,theta_}][s_]:=
    {a*Cos[theta]/E^(Sqrt[s^2]/a) + a*(Sqrt[1 - E^((-2*Sqrt[s^2])/a)] - 
      ArcTanh[Sqrt[1 - E^(-2*Sqrt[s^2]/a)]])*Sin[theta],
      -a*Sqrt[1 - E^(-2*Sqrt[s^2]/a)]*Cos[theta] +
        a*ArcTanh[Sqrt[1 - E^(-2*Sqrt[s^2]/a)]]*Cos[theta] +
       (a*Sin[theta])/E^(Sqrt[s^2]/a)}

    tractrixunitspeed[{a_,theta_}][smin_,smax_,opt___]:=
        pplot[tractrixunitspeed[{a,theta}][s],{s,smin,smax},opt]

    triangle[t_]:= If[t<1,t*{1,0},
         If[t<2,{1,0} + (t - 1)*{-1/2,Sqrt[3]/2},
         {1/2,Sqrt[3]/2} + (t - 2)*{-1/2,-Sqrt[3]/2}]]

    triangle[tmin_,tmax_,opt___]:=
        ParametricPlot[Evaluate[triangle[t]],{t,tmin,tmax},opt,
        Axes->None,AspectRatio->Automatic,
        PlotStyle->{{RGBColor[1,0,0],AbsoluteThickness[2]}}]

    trident[a_,b_,c_,d_][t_]:= {t,a*t^2 + b*t + c + d/t}

    trident[a_,b_,c_,d_][tmin_,tmax_,opt___]:=
        pplot[trident[a,b,c,d][t],
            {t,tmin,tmax},opt,PlotRange->{Automatic,{-4,4}}]

    trident[{a_,b_,c_,d_,theta_}][t_]:=
    {(t^2*Cos[theta] - (d + t*(c + b*t + a*t^2))*Sin[theta])/t,
     (c + d/t + t*(b + a*t))*Cos[theta] + t*Sin[theta]}

    trident[{a_,b_,c_,d_,theta_}][tmin_,tmax_,opt___]:=
        pplot[trident[{a,b,c,d,theta}][t],
            {t,tmin,tmax},opt,PlotRange->{Automatic,{-4,4}}]

    trisectrix[a_][t_]:=
        {a*(1 - 4*Cos[t]^2),a*(1 - 4*Cos[t]^2)*Tan[t]}

    trisectrix[a_][tmin_,tmax_,opt___]:=
        pplot[trisectrix[a][t],{t,tmin,tmax},opt]

    trisectrix[{a_,theta_}][t_]:=
    {-a*Cos[t + theta]*(1 + 2*Cos[2*t])*Sec[t],
     -a*Sin[t + theta]*(1 + 2*Cos[2*t])*Sec[t]}

    trisectrix[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[trisectrix[{a,theta}][t],{t,tmin,tmax},opt]

    tschirnhausen[n_,a_][t_]:= {a*Cos[t]/Cos[t/3]^n,
                                a*Sin[t]/Cos[t/3]^n}

    tschirnhausen[n_,a_][tmin_,tmax_,opt___]:=
        pplot[tschirnhausen[n,a][t],{t,tmin,tmax},opt]

    tschirnhausen[{n_,a_,theta_}][t_]:=
    {a*Cos[t + theta]/Cos[t/3]^n,a*Sin[t + theta]/Cos[t/3]^n}

    tschirnhausen[{n_,a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[tschirnhausen[{n,a,theta}][t],{t,tmin,tmax},opt]

    watt[a_,b_,c_][t_]:=
        {Sqrt[b^2 - (a*JacobiCN[t,a^2/c^2] +
              c*JacobiDN[t,a^2/c^2])^2]*Cos[t],
         Sqrt[b^2 - (a*JacobiCN[t,a^2/c^2] +
              c*JacobiDN[t,a^2/c^2])^2]*Sin[t]}

    watt[a_,b_,c_][tmin_,tmax_,opt___]:=
        pplot[watt[a,b,c][t],{t,tmin,tmax},opt]

    watt[{a_,b_,c_,theta_}][t_]:=
    {Cos[t + theta]*Sqrt[b^2 - (a*JacobiCN[t,a^2/c^2] +
     c*JacobiDN[t,a^2/c^2])^2],
     Sqrt[b^2 - (a*JacobiCN[t,a^2/c^2] + 
     c*JacobiDN[t,a^2/c^2])^2]*Sin[t + theta]}

    watt[{a_,b_,c_,theta_}][tmin_,tmax_,opt___]:=
        pplot[watt[{a,b,c,theta}][t],{t,tmin,tmax},opt]

    zcprofile[pm_,a_,c_][t_]:=
       {t,pm*Sqrt[a^2*t^2 - c^2] + c*ArcSin[c/(a*t)]}

    zcprofile[pm_,a_,c_][tmin_,tmax_,opt___]:=
        pplot[zcprofile[pm,a,c][t],{t,tmin,tmax},opt]
        
    zcprofile[{pm_,a_,c_,theta_}][t_]:=
    {t*Cos[theta] - pm*Sqrt[-c^2 + a^2*t^2]*Sin[theta] -
     c*ArcSin[c/(a*t)]*Sin[theta],
     pm*Sqrt[-c^2 + a^2*t^2]*Cos[theta] + 
     c*ArcSin[c/(a*t)]*Cos[theta] + t*Sin[theta]}

    zcprofile[{pm_,a_,c_,theta_}][tmin_,tmax_,opt___]:=
        pplot[zcprofile[{pm,a,c,theta}][t],{t,tmin,tmax},opt]
        
    zetaspiral[a_][t_]:= {a*Re[Zeta[1/2 + I*t]],a*Im[Zeta[1/2 + I*t]]}

    zetaspiral[a_][tmin_,tmax_,opt___]:=
        pplot[zetaspiral[a][t],{t,tmin,tmax},opt]

    zetaspiral[{a_,theta_}][t_]:=
    {a*(Cos[theta]*Re[Zeta[1/2 + I*t]] -  Im[Zeta[1/2 + I*t]]*Sin[theta]),
     a*(Cos[theta]*Im[Zeta[1/2 + I*t]] +  Re[Zeta[1/2 + I*t]]*Sin[theta])}

    zetaspiral[{a_,theta_}][tmin_,tmax_,opt___]:=
        pplot[zetaspiral[{a,theta}][t],{t,tmin,tmax},opt]

(* Implicitly Defined Plane Curves *)

    agnesiimplicit[a_][x_,y_]:= x^2*y - 4*a^2*(2*a - y)

    agnesiimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[agnesiimplicit[a][x,y],{x,xmin,xmax},opt]

    astroidimplicit[n_,a_,b_][x_,y_]:=
              (x/a)^(2/n) + (y/b)^(2/n) - 1

    bulletnoseimplicit[a_,b_][x_,y_]:= (a/x)^2 - (b/y)^2 - 1

    bulletnoseimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[bulletnoseimplicit[a,b][x,y],{x,xmin,xmax},opt]

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

    cardioidimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[cardioidimplicit[a][x,y],{x,xmin,xmax},opt]

    cartesianimplicit[a_,c_,m_][x_,y_]:=
            ((1 - m^2)*(x^2 + y^2) + 2*m^2*c*x + a^2 - m^2*c^2)^2 -
              4*a^2*(x^2 + y^2)

    cartesianimplicit[a_,c_,m_][{xmin_,xmax_,opt___}]:=
        impplot[cartesianimplicit[a,c,m][x,y],{x,xmin,xmax},opt]

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

    cassiniimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[cassiniimplicit[a,b][x,y],{x,xmin,xmax},opt,PlotRange->All]

    circleimplicit[a_][x_,y_]:= x^2 + y^2 - a^2

    circleimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[circleimplicit[a][x,y],{x,xmin,xmax},opt]

    cissoidimplicit[a_][x_,y_]:= x^3 + x*y^2 - 2*a*y^2

    cissoidimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[cissoidimplicit[a][x,y],{x,xmin,xmax},opt]

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

    conchoidimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[conchoidimplicit[a,b][x,y],{x,xmin,xmax},opt]

    crosscurveimplicit[a_,b_][x_,y_]:= (a/x)^2 + (b/y)^2 - 1

    crosscurveimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[crosscurveimplicit[a,b][x,y],{x,xmin,xmax},opt,
        PlotPoints->150]

    cubicparabolaimplicit[a_,b_,c_,d_][x_,y_]:=
          a*x^3 + b*x^2 + c*x + d - y

    cubicparabolaimplicit[a_,b_,c_,d_][{xmin_,xmax_,opt___}]:=
        impplot[cubicparabolaimplicit[a,b,c,d][x,y],{x,xmin,xmax},opt]

    deltoidimplicit[a_][x_,y_]:= (x^2 + y^2)^2 -
         8*a*x*(x^2 - 3*y^2) + 18*a^2*(x^2 + y^2) - 27*a^4

    deltoidimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[deltoidimplicit[a][x,y],{x,xmin,xmax},opt]

    devilimplicit[a_,b_][x_,y_]:= y^2*(y^2 - b^2) - x^2*(x^2 - a^2)

    devilimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[devilimplicit[a,b][x,y],{x,xmin,xmax},opt]

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

    duererimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[duererimplicit[a,b][x,y],{x,xmin,xmax},opt]

    ellipseimplicit[a_,b_][x_,y_]:= x^2/a^2 + y^2/b^2 - 1

    ellipseimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[ellipseimplicit[a,b][x,y],{x,xmin,xmax},opt]

    foliumimplicit[x_,y_]:= x^3 + y^3 - 3*x*y

    foliumimplicit[n_,a_,b_,c_][x_,y_]:=
             a*x^(2*n + 1) + b*y^(2*n + 1) - c*x^n*y^n

    foliumimplicit[n_,a_,b_,c_][{xmin_,xmax_,opt___}]:=
        impplot[foliumimplicit[n,a,b,c][x,y],{x,xmin,xmax},opt]

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

    hippopedeimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[hippopedeimplicit[a,b][x,y],{x,xmin,xmax},opt]

    hyperbolaimplicit[a_,b_][x_,y_]:= x^2/a^2 - y^2/b^2 - 1

    hyperbolaimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[hyperbolaimplicit[a,b][x,y],{x,xmin,xmax},opt]

    kappacurveimplicit[a_][x_,y_]:= (x^2 + y^2)*y^2 - a^2*x^2

    kappacurveimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[kappacurveimplicit[a][x,y],{x,xmin,xmax},opt]

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

    keplerimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[keplerimplicit[a,b][x,y],{x,xmin,xmax},opt]

    lemniscateimplicit[a_][x_,y_]:=
               (x^2 + y^2)^2 - a^2*(x^2 - y^2)

    lemniscateimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[lemniscateimplicit[a][x,y],{x,xmin,xmax},opt]

    lemniscateimplicit[xp_,yp_][x_,y_]:=
	Product[(x - xp[[j]])^2 + (y - yp[[j]])^2,{j,Length[xp]}]

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

    limaconimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[limaconimplicit[a,b][x,y],{x,xmin,xmax},opt]

    nephroidimplicit[a_][x_,y_]:=
           (x^2 + y^2 - 4*a^2)^3 - 108*a^4*y^2

    nephroidimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[nephroidimplicit[a][x,y],{x,xmin,xmax},opt]

    newtonphillipsimplicit[m_,n_,a_,b_][x_,y_]:=
        Sin[y]*Sin[m*y] - a*Sin[x]*Sin[n*x] - b

    parabolaimplicit[a_][x_,y_]:= x^2 - 4*a*y

    parabolaimplicit[a_][{xmin_,xmax_,opt___}]:=
        impplot[parabolaimplicit[a][x,y],{x,xmin,xmax},opt]

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

    perseusimplicit[a_,b_,c_][{xmin_,xmax_,opt___}]:=
        impplot[perseusimplicit[a,b,c][x,y],{x,xmin,xmax},opt]

    piriformimplicit[a_,b_][x_,y_]:= a^4*y^2 - b^2*x^3*(2*a - x)

    piriformimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[piriformimplicit[a,b][x,y],{x,xmin,xmax},opt]

    predatorpreyimplicit[a_,b_,p_,q_][x_,y_]:=
             y^a*x^p*E^(-b*y - q*x)

    predatorpreyimplicit[a_,b_,p_,q_][xmin_,xmax_,ymin_,ymax_,opt___]:=
        ContourPlot[Evaluate[predatorpreyimplicit[a,b,p,q][x,y]],
        {x,xmin,xmax},{y,ymin,ymax},opt,PlotPoints->50,
         Frame->False,AspectRatio->Automatic,ColorFunction->Hue]

    scimplicit[x_,y_]:= y^2 - x^3

    scimplicit[{xmin_,xmax_,opt___}]:=
        impplot[scimplicit[x,y],{x,xmin,xmax},opt]

    semicubicimplicit[a_,b_,c_,d_][x_,y_]:=
                      a*x^3 + b*x^2 + c*x + d - y^2

    semicubicimplicit[a_,b_,c_,d_][{xmin_,xmax_,opt___}]:=
        impplot[semicubicimplicit[a,b,c,d][x,y],{x,xmin,xmax},opt]

    serpentineimplicit[a_,b_][x_,y_]:= y - a*x + b*x^2*y

    serpentineimplicit[a_,b_][{xmin_,xmax_,opt___}]:=
        impplot[serpentineimplicit[a,b][x,y],{x,xmin,xmax},opt,
        PlotPoints->150]

    strophoidimplicit[a_][x_,y_]:= y^2*(a - x) - x^2*(a + x)

    strophoidimplicit[m_,a_][x_,y_]:= y^2*(a - x) - x^2*(a*m + x)

    strophoidimplicit[m_,a_][{xmin_,xmax_,opt___}]:=
        impplot[strophoidimplicit[m,a][x,y],{x,xmin,xmax},opt]

    tridentimplicit[a_,b_,c_,d_][x_,y_]:=
                      a*x^3 + b*x^2 + c*x + d - x*y

    tridentimplicit[a_,b_,c_,d_][{xmin_,xmax_,opt___}]:=
        impplot[tridentimplicit[a,b,c,d][x,y],{x,xmin,xmax},opt]

    wattimplicit[a_,b_,c_][x_,y_]:= (x^2 + y^2)^3 -
        2*(a^2 + b^2 - c^2)^2*(x^2 + y^2)^2 +
        ((a^2 + b^2 - c^2)^4 + 4*a^2y^2)*(x^2 + y^2) -
        4*a^2*b^2*y^2

(* Polar Defined Plane Curves *)

    archimedesspiralpolar[n_,a_][theta_]:= a*theta^(1/n)

    archimedesspiralpolar[n_,a_][thetamin_,thetamax_,opt___]:=
        polarpplot[archimedesspiralpolar[n,a][theta],
        {theta,thetamin,thetamax},opt]

    bowpolar[a_][theta_]:= a*(1 - Tan[theta]^2)

    bowpolar[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[bowpolar[a][theta],{theta,thetamin,thetamax}]

    cardioidpolar[a_][theta_]:= 2*a*(1 + Cos[theta])

    cardioidpolar[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[cardioidpolar[a][theta],{theta,thetamin,thetamax},opt]

    cassinipolar[a_,b_,pm1_,pm2_][theta_]:=
           pm1*Sqrt[a^2*Cos[2*theta] +
           pm2*Sqrt[b^4 - a^4*Sin[2*theta]^2]]
       
    cassinipolar[a_,b_,pm1_,pm2_][thetamin_,thetamax_,opt___]:=
        polarpplot[cassinipolar[a,b,pm1,pm2][theta],
        {theta,thetamin,thetamax}]

    cochleoidpolar[n_,a_][theta_]:= a*(Sin[theta]/theta)^n

    cochleoidpolar[n_,a_][thetamin_,thetamax_,opt___]:=
        polarpplot[cochleoidpolar[n,a][theta],
        {theta,thetamin,thetamax}]

    conchoidpolar[a_,b_][theta_]:= a + b*Sec[theta]

    conchoidpolar[a_,b_][thetamin_,thetamax_,opt___]:=
        polarpplot[conchoidpolar[a,b][theta],
        {theta,thetamin,thetamax}]

    epispiralpolar[n_,a_][theta_]:= a/Sin[n*theta]

    epispiralpolar[n_,a_][thetamin_,thetamax_,opt___]:=
        polarpplot[epispiralpolar[n,a][theta],
        {theta,thetamin,thetamax}]

    fermatspiralpolar[a_][theta_]:= a*Sqrt[theta]

    fermatspiralpolar[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[fermatspiralpolar[a][theta],
        {theta,thetamin,thetamax},opt]

    foliumpolar[theta_]:=
         3*Cos[theta]*Sin[theta]/(Cos[theta]^3 + Sin[theta]^3)

    foliumpolar[n_,a_,b_,c_][theta_]:=
         c*Cos[theta]^n*Sin[theta]^n/
         (a*Cos[theta]^(2*n + 1) + b*Sin[theta]^(2*n + 1))

    foliumpolar[n_,a_,b_,c_][thetamin_,thetamax_,opt___]:=
        polarpplot[foliumpolar[n,a,b,c][theta],
        {theta,thetamin,thetamax},opt]

    hyperbolicspiralpolar[a_][theta_]:= a/theta

    hyperbolicspiralpolar[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[hyperbolicspiralpolar[a][theta],
        {theta,thetamin,thetamax},opt]

    kampylepolar[a_][theta_]:= a*Sec[theta]^2

    kampylepolar[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[kampylepolar[a][theta],
        {theta,thetamin,thetamax},opt]

    keplerorbitpolar[a_,e_][theta_]:= a/(1 - e*Cos[theta])

    keplerorbitpolar[a_,e_][thetamin_,thetamax_,opt___]:=
        polarpplot[keplerorbitpolar[a,e][theta],
        {theta,thetamin,thetamax},opt]

    lemniscatepolar[a_][theta_]:= a*Sqrt[Cos[2*theta]]

    lemniscatepolar[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[lemniscatepolar[a][theta],
        {theta,thetamin,thetamax},opt]

    limaconpolar[a_,b_][theta_]:= 2*a*Cos[theta] + b

    limaconpolar[n_,a_,b_][theta_]:= 2*a*Cos[n*theta] + b

    limaconpolar[n_,a_,b_][thetamin_,thetamax_,opt___]:=
        polarpplot[limaconpolar[n,a,b][theta],
        {theta,thetamin,thetamax},opt]

    lituuspolar[a_][theta_]:= a/Sqrt[theta]

    lituuspolar[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[lituuspolar[a][theta],
        {theta,thetamin,thetamax},opt]

    logspiralpolar[a_,b_][theta_]:= a*E^(b*theta)

    logspiralpolar[a_,b_][thetamin_,thetamax_,opt___]:=
        polarpplot[logspiralpolar[a,b][theta],{theta,thetamin,thetamax},opt]

    nephroidpolar[a_][theta_]:=
        2*a*(Sin[theta/2]^(2/3) + Cos[theta/2]^(2/3))^(3/2)
 
    nephroidpolar[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[nephroidpolar[a][theta],{theta,thetamin,thetamax},opt]

    pacman[n_][theta_]:= 1 + Cos[theta]^n

    pacman[a_][thetamin_,thetamax_,opt___]:=
        polarpplot[pacman[a][theta],{theta,thetamin,thetamax},opt]

    rosepolar[n_,a_][theta_]:= a*Cos[n*theta]

    rosepolar[n_,a_][thetamin_,thetamax_,opt___]:=
        polarpplot[rosepolar[n,a][theta],
        {theta,thetamin,thetamax},opt]

    tschirnhausenpolar[n_,a_][theta_]:= a*Cos[theta/3]^-n

    tschirnhausenpolar[n_,a_][thetamin_,thetamax_,opt___]:=
        polarpplot[tschirnhausenpolar[n,a][theta],
        {theta,thetamin,thetamax},opt]

(* Parametrically Defined Space Curves *)

    ast3d[n_,a_,b_][t_]:= {a*Cos[t]^n,b*Sin[t]^n,Cos[2*t]}

    ast3d[n_,a_,b_][tmin_,tmax_,opt___]:=
        scplot[ast3d[n,a,b][t],{t,tmin,tmax},opt]

    baseballseam[a_,c_][t_]:=
        {a*Sin[Pi/2 - (Pi/2 - c)*Cos[t]]*Cos[t/2 + c*Sin[2*t]],
         a*Sin[Pi/2 - (Pi/2 - c)*Cos[t]]*Sin[t/2 + c*Sin[2*t]],
         a*Cos[Pi/2 - (Pi/2 - c)*Cos[t]]}

    baseballseam[a_,c_][tmin_,tmax_,opt___]:=
        scplot[baseballseam[a,c][t],{t,tmin,tmax},opt]

    besselcurve[a_,b_,c_][t_]:=
         {BesselJ[a,t],BesselJ[b,t],BesselJ[c,t]}

    besselcurve[a_,b_,c_][tmin_,tmax_,opt___]:=
        scplot[besselcurve[a,b,c][t],{t,tmin,tmax},opt]

    bicylinder[a_,b_,pm_][t_]:= {a*Cos[t],a*Sin[t],
                                 pm*Sqrt[b^2 - a^2*Sin[t]^2]}

    bicylinder[a_,b_,pm_][tmin_,tmax_,opt___]:=
        scplot[bicylinder[a,b,pm][t],{t,tmin,tmax},opt]

    biquadratic[a_,k_][t_]:=
           {a*JacobiSN[t,k],a*JacobiCN[t,k],a*JacobiDN[t,k]}

    biquadratic[a_,k_][tmin_,tmax_,opt___]:=
        scplot[biquadratic[a,k][t],{t,tmin,tmax},opt]

    clelia[a_][m_,n_][t_]:=
           {a*m*Cos[t]*Sin[n*t],a*m*Sin[t]*Sin[n*t],
            a*Sqrt[1 - m^2*Sin[n*t]^2]}

    clelia[a_][m_,n_][tmin_,tmax_,opt___]:=
        scplot[clelia[a][m,n][t],{t,tmin,tmax},opt]

    conicalhelix[a_,b_,c_,d_][t_]:=
        {a*t*Cos[d*t],b*t*Sin[d*t],c*t}

    conicalhelix[a_,b_,c_,d_][tmin_,tmax_,opt___]:=
        scplot[conicalhelix[a,b,c,d][t],
            {t,tmin,tmax},opt,PlotPoints->1000]

    cubicalellipse[a_,b_,c_,d_,e_,f_,pm_][t_]:=
           {a*t/(1 + pm*t^2),
           (b*t + c*t^2)/(1 + pm*t^2),
           (d*t + e*t^2 + f*t^3)/(1 + pm*t^2)}

    cubicalellipse[a_,b_,c_,d_,e_,f_,pm_][tmin_,tmax_,opt___]:=
        scplot[cubicalellipse[a,b,c,d,e,f,pm][t],
            {t,tmin,tmax},opt]

    eightknot[t_]:=
        {10*(Cos[t] + Cos[3*t]) + Cos[2*t] + Cos[4*t],
         6*Sin[t] + 10*Sin[3*t],
         4*Sin[3*t]*Sin[5*t/2] + 4*Sin[4*t] - 2*Sin[6*t]}

    eightknot[tmin_,tmax_,opt___]:=
        scplot[eightknot[t],{t,tmin,tmax},opt]

    epitrochoid3d[a_,b_,h_,omega_][t_]:=
        {Cos[t]*(a + b*Cos[omega]) - 
         h*(Cos[t]*Cos[a*t/b]*Cos[omega] - Sin[t]*Sin[a*t/b]),
         Sin[t]*(a + b*Cos[omega]) - 
         h*(Sin[t]*Cos[a*t/b]*Cos[omega] + Cos[t]*Sin[a*t/b]),
        (b - h*Cos[a*t/b])*Sin[omega]}

    epitrochoid3d[a_,b_,h_,omega_][tmin_,tmax_,opt___]:=
        scplot[epitrochoid3d[a,b,h,omega][t],
            {t,tmin,tmax},opt]

    genhelix[a_,b_][f_][t_]:= {a*Cos[t],b*Sin[t],f[t]}

    genhelix[a_,b_][f_][tmin_,tmax_,opt___]:=
        scplot[genhelix[a,b][f][t],{t,tmin,tmax},opt]

    helix[a_,c_][t_]:= {a*Cos[t],a*Sin[t],c*t}

    helix[a_,b_,c_][t_]:= {a*Cos[t],b*Sin[t],c*t}

    helix[a_,b_,c_][tmin_,tmax_,opt___]:=
        scplot[helix[a,b,c][t],{t,tmin,tmax},opt]

    helixunitspeed[a_,b_][s_]:=
        {a*Cos[s/Sqrt[a^2 + b^2]],a*Sin[s/Sqrt[a^2 + b^2]],
         b*s/Sqrt[a^2 + b^2]}

    helixunitspeed[a_,b_][smin_,smax_,opt___]:=
        scplot[helixunitspeed[a,b][s],{s,smin,smax},opt]

    horopter[a_,b_][t_]:=
        {2*a/(1 + b^2*t^2),2*a*b*t/(1 + b^2 t^2),t}

    horopter[a_,b_][tmin_,tmax_,opt___]:=
        scplot[horopter[a,b][t],{t,tmin,tmax},opt]

    line3d[a1_,b1_,c1_,a2_,b2_,c2_][t_]:=
         {a1*t + a2,b1*t + b2,c1*t + c2}

    lissajous3d[n_,d_,a_,b_,c_][t_]:=
         {a*Sin[n*t + d],b*Sin[t],c*Cos[t]} 

    lissajous3d[n_,d_,a_,b_,c_][tmin_,tmax_,opt___]:=
        scplot[lissajous3d[n,d,a,b,c][t],{t,tmin,tmax},opt]

    oneparametersubgroup[mat_,vec_][t_]:=
        ComplexExpand[Normal[MatrixExp[Rationalize[t*mat]]].vec]

    oneparametersubgroup[mat_,vec_][tmin_,tmax_,opt___]:=
        scplot[oneparametersubgroup[mat,vec][t],
            {t,tmin,tmax},opt,PlotRange->All]

    pseudospheregeodesic[a_,c_,u0_,pm_][t_]:=
        {a*Sin[t]*Cos[u0 + pm*Sqrt[a^2 - 2*c^2 - a^2*Cos[2*t]]/
         (Sqrt[2]*c*Sin[t])],
         a*Sin[t]*Sin[u0 + pm*Sqrt[a^2 - 2*c^2 - a^2*Cos[2*t]]/
         (Sqrt[2]*c*Sin[t])],a*(Cos[t] + Log[Tan[t/2]])}

    pseudospheregeodesic[a_,c_,u0_,pm_][tmin_,tmax_,opt___]:=
        scplot[pseudospheregeodesic[a,c,u0,pm][t],
            {t,tmin,tmax},opt]

    pseudosphericalloxodrome[a_,b_][t_]:=
         {a*Cos[b*Cosh[t]]/Cosh[t],
          a*Sin[b*Cosh[t]]/Cosh[t],a*(t - Tanh[t])}

    pseudosphericalloxodrome[a_,b_][tmin_,tmax_,opt___]:=
        scplot[pseudosphericalloxodrome[a,b][t],
            {t,tmin,tmax},opt]

    seiffertspiral[a_,k_][t_]:=
        {a*Cos[t]*JacobiSN[t,k],a*Sin[t]*JacobiSN[t,k],
         a*JacobiCN[t,k]}

    seiffertspiral[a_,k_][tmin_,tmax_,opt___]:=
        scplot[seiffertspiral[a,k][t],
            {t,tmin,tmax},opt]

    sphericalcardioid[a_][t_]:=
        {2*a*Cos[t] - a*Cos[2*t],2*a*Sin[t] - a*Sin[2*t],
         Sqrt[8]*a*Cos[t/2]}

    sphericalcardioid[a_][tmin_,tmax_,opt___]:=
        scplot[sphericalcardioid[a][t],{t,tmin,tmax},opt]

    sphericalellipse[a_,b_,c_,pm_][t_]:=
        {b*JacobiCN[t,(c^2 - b^2)/(a^2 - b^2)],
         c*JacobiSN[t,(c^2 - b^2)/(a^2 - b^2)],
         pm*Sqrt[a^2 - b^2]*JacobiDN[t,(c^2 - b^2)/(a^2 - b^2)]}

    sphericalellipse[a_,b_,c_,pm_][tmin_,tmax_,opt___]:=
        scplot[sphericalellipse[a,b,c,pm][t],
            {t,tmin,tmax},opt]

    sphericalhelix[a_,b_][t_]:=
        {(a + b)*Cos[t] - b*Cos[(a + b)*t/b],
         (a + b)*Sin[t] - b*Sin[(a + b)*t/b],
          2*Sqrt[a*b + b^2]*Cos[a*t/(2*b)]}

    sphericalhelix[a_,b_][tmin_,tmax_,opt___]:=
        scplot[sphericalhelix[a,b][t],{t,tmin,tmax},opt]

    sphericalloxodrome[a_,b_,c_][t_]:=
        {2*a*b*E^(c*t)*Cos[t]/(1 + b^2*E^(2*c*t)),
         2*a*b*E^(c*t)*Sin[t]/(1 + b^2*E^(2*c*t)),
         a*(-1 + b^2*E^(2*c*t))/(1 + b^2*E^(2*c*t))}

    sphericalloxodrome[a_,b_,c_][tmin_,tmax_,opt___]:=
        scplot[sphericalloxodrome[a,b,c][t],{t,tmin,tmax},opt]

    sphericalloxodromeunitspeed[a_,d_,c_][s_]:=
    {a*Cos[c*s/(a*Sqrt[1 + c^2])]*
       Cos[(d - 2*ArcTanh[Tan[c*s/(2*a*Sqrt[1 + c^2])]])/c], 
     -a*Cos[c*s/(a*Sqrt[1 + c^2])]*
       Sin[(d - 2*ArcTanh[Tan[c*s/(2*a*Sqrt[1 + c^2])]])/c], 
      a*Sin[c*s/(a*Sqrt[1 + c^2])]}

    sphericalloxodromeunitspeed[a_,d_,c_][smin_,smax_,opt___]:=
        scplot[sphericalloxodromeunitspeed[a,d,c][s],
        {s,smin,smax},opt]

    sphericalnephroid[a_][t_]:=
        {3*a*Cos[t] - a*Cos[3*t],3*a*Sin[t] - a*Sin[3*t],
         Sqrt[12]*a*Cos[t]}

    sphericalnephroid[a_][tmin_,tmax_,opt___]:=
        scplot[sphericalnephroid[a][t],
            {t,tmin,tmax},opt]

    sphericalspiral[a_][m_,n_][t_]:=
          {a*Cos[m*t]*Cos[n*t],a*Sin[m*t]*Cos[n*t],a*Sin[n*t]}

    sphericalspiral[a_][m_,n_][tmin_,tmax_,opt___]:=
        scplot[sphericalspiral[a][m,n][t],
            {t,tmin,tmax},opt,PlotPoints->1000]

    tennisballseam[n_,a_,b_,c_][t_]:=
          {-a*(Cos[t]*(b^2 + c^2*n^2*Cos[n*t]^2) + 
           c^2*n^3*Sin[t]*Sin[2*n*t]/2)/
           (Sqrt[b^2*Cos[t]^2 + c^2*n^2*Cos[n*t]^2 + 
            a^2*Sin[t]^2]*Sqrt[a^2*b^2*Cos[t]^4 + 
            b^2*c^2*n^2*Cos[n*t]^2*Sin[t]^2 + 
            a^2*b^2*Sin[t]^4 + a^2*c^2*n^4*Sin[t]^2*Sin[n*t]^2 + 
            Cos[t]^2*(a^2*c^2*n^2*Cos[n*t]^2 + 
            b^2*(2*a^2*Sin[t]^2 + c^2*n^4*Sin[n*t]^2)) + 
            a^2*c^2*n^3*Sin[2*t]*Sin[2*n*t]/2 - 
            b^2*c^2*n^3*Sin[2*t]*Sin[2*n*t]/2]),
            b*(-2*(a^2 + c^2*n^2*Cos[n*t]^2)*Sin[t] + 
            c^2*n^3*Cos[t]*Sin[2*n*t])/
            (2*Sqrt[b^2*Cos[t]^2 + c^2*n^2*Cos[n*t]^2 + 
            a^2*Sin[t]^2]*Sqrt[a^2*b^2*Cos[t]^4 + 
            b^2*c^2*n^2*Cos[n*t]^2*Sin[t]^2 + 
            a^2*b^2*Sin[t]^4 + a^2*c^2*n^4*Sin[t]^2*Sin[n*t]^2 + 
            Cos[t]^2*(a^2*c^2*n^2*Cos[n*t]^2 + 
            b^2*(2*a^2*Sin[t]^2 + c^2*n^4*Sin[n*t]^2)) + 
            a^2*c^2*n^3*Sin[2*t]*Sin[2*n*t]/2 - 
            b^2*c^2*n^3*Sin[2*t]*Sin[2*n*t]/2]),
            -c*n*((a^2 - b^2)*Cos[t]*Cos[n*t]*Sin[t] + 
            b^2*n*Cos[t]^2*Sin[n*t] + a^2*n*Sin[t]^2*Sin[n*t])/
            (Sqrt[b^2*Cos[t]^2 + c^2*n^2*Cos[n*t]^2 + 
            a^2*Sin[t]^2]*Sqrt[a^2*b^2*Cos[t]^4 + 
            b^2*c^2*n^2*Cos[n*t]^2*Sin[t]^2 + 
            a^2*b^2*Sin[t]^4 + a^2*c^2*n^4*Sin[t]^2*Sin[n*t]^2 + 
            Cos[t]^2*(a^2*c^2*n^2*Cos[n*t]^2 + 
            b^2*(2*a^2*Sin[t]^2 + c^2*n^4*Sin[n*t]^2)) + 
            a^2*c^2*n^3*Sin[2*t]*Sin[2*n*t]/2 - 
            b^2*c^2*n^3*Sin[2*t]*Sin[2*n*t]/2])}

    tennisballseam[n_,a_,b_,c_][tmin_,tmax_,opt___]:=
        scplot[tennisballseam[n,a,b,c][t],{t,tmin,tmax},opt]

    torusknot[a_,b_,c_][p_,q_][t_]:=
         {(a + b*Cos[q*t])*Cos[p*t],
          (a + b*Cos[q*t])*Sin[p*t],c*Sin[q*t]}

    torusknot[a_,b_,c_][p_,q_][tmin_,tmax_,opt___]:=
        scplot[torusknot[a,b,c][p,q][t],
            {t,tmin,tmax},opt,PlotPoints->1000]

(*
     tprime[a_,b_][v_]:=1/Sqrt[-Cos[v]*(a/b + Cos[v])]


    torusasym[c_,d_,a_,b_][s_]:=
         Module[{uu,tmp1,tmp2,ss},
             tmp1=NDSolve[{uu'[ss]==tprime[a,b][ss],
               uu[c]==d},uu,{ss,0.5001Pi,1.4999Pi},
               MaxSteps->2000];
               tmp2=First[uu[s] /. tmp1];
		{(a + b*Cos[s])*Cos[tmp2],(a + b*Cos[s])*Sin[tmp2],
		b*Sin[s]}]
*)

    twicubic[t_]:= {t,t^2,t^3}

    twicubic[tmin_,tmax_,opt___]:=
        scplot[twicubic[t],{t,tmin,tmax},opt]

    veryflat[1][t_]:= {t,0,5*E^(-1/t^2)}

    veryflat[1][tmin_,tmax_,opt___]:=
        scplot[veryflat[1][t],{t,tmin,tmax},opt]

    veryflat[2][t_]:= If[t<0,{t,5*E^(-1/t^2),0},
                     If[t==0,{0,0,0},{t,0,5*E^(-1/t^2)}]]

    veryflat[2][tmin_,tmax_,opt___]:=
        scplot[veryflat[2][t],{t,tmin,tmax},opt]

    viviani[a_][t_]:= {a*(1 + Cos[t]),a*Sin[t],2*a*Sin[t/2]}

    viviani[a_][tmin_,tmax_,opt___]:=
        scplot[viviani[a][t],{t,tmin,tmax},opt]

    viviani[a_,b_][t_]:=
    {b + (2*a - b)*Cos[t],(2*a - b)*Sin[t],2*Sqrt[b*(2*a - b)]*Sin[t/2]}

    viviani[a_,b_][tmin_,tmax_,opt___]:=
        scplot[viviani[a,b][t],{t,tmin,tmax},opt]

    wigglyellipse[n_,a_,b_,c_][t_]:= {a*Cos[t],b*Sin[t],c*Sin[n*t]}

    wigglyellipse[n_,a_,b_,c_][tmin_,tmax_,opt___]:=
        scplot[wigglyellipse[n,a,b,c][t],{t,tmin,tmax},opt]

(* Parametrically Defined Curves in R^n*)

    besseln[n_,a_][t_]:= Table[a*BesselJ[k,t],{k,0,n - 1}]

    bernstein[p_,q_][t_]:=Binomial[p,q]*t^q*(1 - t)^(p - q)
    
    beziercurve[pts_][t_]:= Module[{n=Length[pts]},
         Sum[bernstein[n - 1,k][t]*pts[[k + 1]],
         {k,0,n - 1}]]

         cosn[n_,a_][t_]:= Table[a*Sin[k*t],{k,1,n}]

    hermitecurve[p0_,m0_,m1_,p1_][t_]:=
        (bernstein[3,0][t] + bernstein[3,1][t])*p0 +
        bernstein[3,1][t]*m0 - bernstein[3,2][t]*m1 +
        (bernstein[3,2][t] + bernstein[3,3][t])*p1

    ell[ts_,i_,t_]:=Product[(ts[[j]]-t)/(ts[[j]]-ts[[i]]),{j,1,i-1}]*
        Product[(ts[[j]]-t)/(ts[[j]]-ts[[i]]),{j,i+1,Length[ts]}]
		
    lagrangecurve[pts_,ts_][t_]:= Module[{n=Length[pts]},
         Sum[ell[ts,k+1,t]*pts[[k+1]],{k,0,n - 1}]]

    sinn[n_,a_][t_]:= Table[a*Sin[k*t],{k,1,n}]

    twistedn[n_,a_][t_]:= Table[a*t^k,{k,1,n}]

(* Approximations *)

    regularpolygon[n_,a_]:=
        approxcurve[0,2*Pi,n][circle[a]]

    teeth[m_,t_][shortlines_]:=
	    Nest[Flatten[(tooth[t]/@#)]&,shortlines,m]

(*
    tooth[t_][Line[{{a_,b_},{c_,d_}}]]:=
        {Line[{{a,b},(2/3)*{a,b} + (1/3)*{c,d}}],
         Line[{(2/3)*{a,b} + (1/3)*{c,d},
         (t/2)*{a - b + c + d,a + b - c + d} +
         ((1 - t)/2)*{a + b + c - d,-a + b + c + d}}],
         Line[{(t/2)*{a - b + c + d,a + b - c + d} +
         ((1 - t)/2)*{a + b + c - d,-a + b + c + d},
         (1/3)*{a,b} + (2/3)*{c,d}}],
        Line[{(1/3)*{a,b} + (2/3)*{c,d},{c,d}}]}
*)

tooth[t_][Line[{{a_,b_},{c_,d_}}]]:=
        {Line[{{a,b},(2/3)*{a,b} + (1/3)*{c,d}}],
         Line[{(2/3)*{a,b} + (1/3)*{c,d},
         {(a + c)/2 + t*(b - d),(b + d)/2 +t*(c - a)}}],
         Line[{{(a + c)/2 + t*(b - d),(b + d)/2 +t*(c - a)},
         (1/3)*{a,b} + (2/3)*{c,d}}],
        Line[{(1/3)*{a,b} + (2/3)*{c,d},{c,d}}]}

   kochsnowflake[m_]:=
    teeth[m,-Sqrt[3]/6][linetolines[regularpolygon[3,1]]]

   kochsnowflake[m_,t_]:=
    teeth[m,t][linetolines[regularpolygon[3,1]]]
    
End[]

EndPackage[]