(***************************************************************************

            REMEMBER: To run programs involving
            curves, you need to input
            radecas.m

**************************************************************************)


(*  Control polygons for poly and rational curves
    and examples of calls to the various programs to draw them            
             11/15/2002                                                  

All points need a weight. 
If poly curve, then weight = 1


                                                                           *)


(**************************************************************************)
(*  A rational version of the de Casteljau algorithm.                     *)
(*  Control polygon     cpoly = (b0,...,bm)                               *)
(*  Each control point is either a weighed point {x1,...xn,w}, w <> 0,
*  or a control vector  {x1,...xn,0}.
*  Affine frame on real line: (r, s)
*  pdecas[cpoly, r, s, t]: Computation of the current point F(t) 
*  To display pdecas, use showpic2D, or showpic3D
*  gsubdiv[cpoly, r, s, l, m, n]: perform n steps  of subdivision.
*  starting with t = l, and ending with t = m. First, recomputes new control
*  polygon, and then use interval [l, m], and subdivide with $t = (l + m)/2.
*  To display gsubdiv, use subdiv2D, or subdiv3D
*  subdiv[cpoly, r, s, n]: perform n steps  of subdivision (using t = (r + 2)/2)
*  To display subdiv, use showsub2D, or showsub3D
*  To display subdiv with flagpol, use showsub2, or showsub3
*  ggdecas[cpoly, tt, r, s]: computes the polar value f(t1,...,tm),
*  where tt = {t1,...,tm}, using the control polygon cpoly with m edges
*  To display ggdecas, use showc2D, or showc3D
*  showclo2D[cpoly, r, s, lx, mx, ly, my, n] will render
*  a closed rational curve, performing n subdivision steps
*  This program first computes a new control polygon wrt [-1, 1],
*  and then computes the opposite control polygon (In the hat space)
*  where oppoly[[i]] = (-1)^(m - i) poly[[i]]    
*  showclo3D[cpoly, r, s, lx, mx, ly, my, lz, mz, n] is the 3D version
*  rendclo2D[cpoly, r, s, lx, mx, ly, my, flagpol, n] will render
*  a closed rational curve, performing n subdivision steps
*  If flagpol =!= 1, do not display control polygon       
* If flagpol === -1, show original arc in black, and dual arc in red     
*  rendlo3D[cpoly, r, s, lx, mx, ly, my, lz, mz, flagpol, n] is the 3D version        
***************************************************************************)

(* Display["rose3.ps", Out[13], "EPS"] *)

(*******************************************************)
Quarter of a circle
inpol = {{1,0,1},{1,1,1},{0,1,2}};
showpic2D[inpol,0,1,-2,2,-2,2,0.5];
showsub2D[inpol,0,1,-2,2,-2,2,4];
showclo2D[inpol,0,1,-2,2,-2,2,6];
rendclo2D[inpol,0,1,-2,2,-2,2,1,6];
rendcloz2D[inpol,0,1,-2,2,-2,2,1,6];

Quarter of ellipse
inpol = {{2,0,1},{2,1,1},{0,1,2}}
showsub2D[inpol,0,1,-2,2,-2,2,4];
subdiv2D[inpol,0,1,-4,4,-2,2,-2,2,6]
showclo2D[inpol,0,1,-2,2,-2,2,6]
rendclo2D[inpol,0,1,-2,2,-2,2,1,6]
rendcloz2D[inpol,0,1,-2,2,-2,2,1,6]

Another example (a cubic with asymptote)
inpol := {{0,0,1},{2,6,1},{6,8,2},{10,0,1}};
showc2D[inpol, tt, 0, 1, 0, 10, 0, 8]
showsub2D[inpol,0,1,0,10,0,8,4]
subdiv2D[inpol,0,1,-2,3/2,-2,20,-10,20,4]
rendclo2D[inpol,0,1,-2,15,-6,16,1,6]
rendclo2D[inpol,0,1,-2,15,-6,16,-1,6]
rendcloz2D[inpol,0,1,-2,15,-6,16,1,6]

And another
inpol := {{0,0,4},{2,6,1},{6,8,0.1},{10,0,1}};
showsub2D[inpol,0,1,-2,20,-10,20,4]
rendclo2D[inpol,0,1,-8,20,-10,20,1,6]
rendcloz2D[inpol,0,1,-8,20,-10,20,1,6]

Folium of Descartes
inpol := {{0,0,1},{1,0,1},{2,1,1},{1.5,1.5,2}}
showsub2D[inpol,0,1,-4,4,-4,4,4];
rendclo2D[inpol,0,1,-4,4,-4,4,0,4] (* interesting *)
rendclo2D[inpol,0,1,-10,10,-10,10,1,4] (* interesting *)

Lemniscate of Bernoulli
inpol := {{0,0,1},{0.25,0.25,1},{0.5,0.5,1},{1,0.5,1},{1,0,2}};
inpol := {{0,0,1},{1/4,1/4,1},{1/2,1/2,1},{1,1/2,1},{1,0,2}};
showsub2D[inpol,0,1,-2,2,-1,1,4]
subdiv2D[inpol,0,1,-2,2,-2,2,-1,1,6]
subdiv2D[inpol,0,1,-1,1,-2,2,-1,2.1,6]
subdiv2D[inpol,0,1,0,1,-1,1.1,-0.5,0.5,6]
subdiv2D[inpol,0,1,-0.5,1,-1,1.1,-0.5,1,6]
showcur2D[inpol,0,1,-2,2,-1,1,-8,8,160]
showcur2D[inpol,0,1,-2,2,-1,1,-8,8,50]
showclo2D[inpol,0,1,-2,2,-1,1,6]
rendclo2D[inpol,0,1,-1.1,1.1,-0.38,0.38,0,6]
rendcloz2D[inpol,0,1,-1.1,1.1,-0.5,0.5,1,6]

(* closed quartic, standard example  *)
inpol := {{0,0,1},{2,6,1},{6,8,2},{10,4,1},{10,0,1}};
(* complementary arc  *)
inpolc := {{0,0,1},{2,6,-1},{6,8,2},{10,4,-1},{10,0,1}};
showcur2D[inpol,0,1,-2,15,-3,12,-20,20,400]
subdiv2D[inpol,0,1,-20,20,-2,15,-3,12,8]
subdiv2D[inpol,0,1,-2,2,-2,15,-3,12,7]

(* arc over [0, 1] *)
subdiv2D[inpol,0,1,0,1,-2,13,-3,12,7]
(* complementary arc over [0, 1] *)
subdiv2D[inpolc,0,1,0,1,-2,13,-3,12,7]

(* arc treated over [-1, 1] *)
subdiv2Db[inpol,-1,1,-1,1,-2,13,-3,12,7,1]
(* complementary arc treated over [-1, 1] *)
subdiv2Db[inpolc,-1,1,-1,1,-2,13,-3,12,7,2]

showclo2D[inpol,0,1,-2,13,-3,12,6]
rendclo2D[inpol,0,1,-2,13,-3,12,0,6]
rendclo2D[inpol,0,1,-2,13,-3,12,1,6]
rendcloz2D[inpol,0,1,-2,13,-3,12,1,7]


(* Another closed quartic *)
inpol =  {{0, 0, 2}, {2, 6, 1}, {6, 8, 2}, {10, 4, 1}, {10, 0, 2}}
showclo2D[inpol,-1,1,-1,11,-4,9,6]

(* Another closed quartic *)
inpol =  {{0, 0, 1}, {2, 6, 0}, {6, 8, 0.333}, {10, 4, 0}, {10, 0, 1}}
showclo2D[inpol,-1,1,-10,22,-8,15,6]

Half Circle (control point at infinity):
inpol := {{1,0,1},{0,1,0},{-1,0,1}};
showsub2D[inpol,0,1,-4,4,-4,4,4]

Full circle (control points at infinity):
inpol := {{-1,0,1},{0,1,0},{3,0,1/3},{0,-1,0},{-1,0,1}};
showsub2D[inpol,0,1,-3,3,-3,3,5]

Full circle:
inpol := {{-1,0,1},{-1,1,1},{1/3,4/3,1},{-1,1,-1},{-1,0,1}};
showsub2D[inpol,0,1,-3,3,-3,3,5]

Segment of Hyperbola (* control point at infinity *):
inpol := {{0,1,0},{0,2,0.5},{1,1,1}};
showsub2D[inpol,0,1,-4,12,-4,12,4] (* infinite for t = 0 *)
showsub2D[inpol,0.1,1,-4,12,-4,12,4]
rendclo2D[inpol,0.1,1,-12,12,-12,12,1,4] (* interesting *)

(* Segment of cuspidal cubic (control points at infinity): *)
inpol := {{0,0,1},{0,0,0},{0.333,0,0},{0,1,0}};
showsub2D[inpol,0,1,0,12,0,12,4] (* infinite for t = 1 *)
showsub2D[inpol,0,0.8,0,12,0,12,4] 
showsub2D[inpol,-10,0.8,0,12,-8,12,4] 


cpoly = {{0, -4, 1}, {10, 15, 20}, {5, -12, 40}, {0, 15, 20}, {10, -4, 1}}


(*******************************************************)
3D curves:

(* Viviani  window with control vectors *)
inpol := {{0,0,-1,1},{1/2,0,0,0},{0,4,1,1/7},{-1/2,0,0,0},{0,-16/3,5/3,3/35},
          {1/2,0,0,0},{0,4,1,1/7},{-1/2,0,0,0},{0,0,-1,1}};
showsub3D[inpol,0,1,-1,3,-6,4,-1,2,4] 
subdiv3D[inpol,0,1,0,1,-1,3,-6,4,-1,2,5]


(* Viviani Window    *)

cpoly =   {{0, 0, 1, 1}, {1/2, 0, 1, 1}, {3/4, 1/2, 3/4, 4/3}, {1/2, 1, 1/2, 2},
   {0, 1, 0, 4}}

showclo3D[cpoly,0,1,-1,3,-6,4,-1,2,5]  (* fine  *)
rendclo3D[cpoly,0,1,-0.55,0.55,0,1,-1,1,1,5]  (* great *)
rendclo3D[cpoly,0,1,-0.55,0.55,0,1,-1,1,-1,6]  (* great *)
rendclo3D[cpoly,0,1,-0.55,0.55,0,1,-1,1,-1,2] 
rendcloz3D[cpoly,0,1,-0.6,0.6,0,1.1,-1,1,-1,6]  (* great *)
rendcloz3D[cpoly,0,1,-0.6,0.8,0,1.1,-1,1,1,6]  (* great *)
rendcloz3D[cpoly,0,1,-0.6,0.75,0,1.2,-1,1,1,3] 


Twisted cubic:
inpol := {{0,0,0,1},{1/3,0,0,1},{2/3,1/3,0,1},{1,1,1,1}};
showsub3D[inpol,-1,2,-1,1.5,-1,1.5,-1,1.5,4] 
showsub3D[inpol,0,1,-1,1,-1,1,-1,1,4] 
subdiv3D[inpol,0,1,-1,1,-1,1,-1,1,-1,1,4] 

Twisted cubic (with control vectors):
inpol := {{0,0,1,0},{0,1/3,0,0},{4/15,0,-1/5,0},{0,-1/2,0,2/5},
          {-4/15,0,1/5,0},{0,1/3,0,0},{0,0,-1,0}};
showsub3D[inpol,0.2,0.8,-3,3,-3,3,-3,3,4] 
showsub3D[inpol,0.2,0.8,-10,10,-10,10,-10,10,4] 

(****************************************************************)
(* This example is a quarter circle as a cubic arc, using
the magic number magic = 4(sqrt(2) - 1)/3
****************************************************************)

magic := 4 * (Sqrt[2] - 1)/3;

cpoly = 
{{4, 0, 1}, {4, N[4 * magic], 1}, {N[4 * magic], 4, 1}, {0, 4, 1}};
 
showpic2D[cpoly, 0, 1, -0.5, 4.5, -0.5, 4.5, 0.4]

showsub2D[cpoly, 0, 1, -0.5, 4.5, -0.5, 4.5, 4]

(*      some polynomial curves    *)

cpoly := {{0, 0, 1}, {3, 8, 1}, {6, 6, 1}, {10, 0, 1}};

showpic2D[cpoly, 0, 1, 0, 10, 0, 10, 0.5]

cpoly := {{0, 0, 1}, {3, 8, 1}, {6, 6, 1}, {10, 0, 1}};

showpic2D[cpoly, 0, 1, -5, 15, -25, 15, -0.5]

showcur2D[cpoly, 0, 1, 0, 10, 0, 10, 0, 1, 20]

showsub2D[cpoly, 0, 1, 0, 10, 0, 8, 1]
showsub2D[cpoly, 0, 1, 0, 10, 0, 8, 2]
showsub2D[cpoly, 0, 1, 0, 10, 0, 8, 3]
showsub2D[cpoly, 0, 1, 0, 10, 0, 8, 4]
showsub2D[cpoly, 0, 1, 0, 10, 0, 8, 5]

(***)

cpoly := {{0, 0, 1}, {1, 8, 1}, {2, 2, 1}, {6, -3, 1}, 
{10, 2, 1}, {11, 8, 1}, {12, 0, 1}};

showcur2D[cpoly, 0, 1, 0, 12, -3, 8, 0, 1, 40]
showcur2D[cpoly, 0, 1, -2, 15, -10, 12, 0, 1, 40]

showpic2D[cpoly, 0, 1, 0, 12, -3, 8, 0.5]
showpic2D[cpoly, 0, 1, -2, 15, -15, 15, -0.05]
showsub2D[cpoly, 0, 1, 0, 12, -3, 8, 1]
showsub2D[cpoly, 0, 1, 0, 12, -3, 8, 5]


(***)

cpoly := {{2, -3, 1}, {0, 0, 1}, {1, 8, 1}, {2, 2, 1}, 
{6, -3, 1}, {10, 2, 1}, {11, 8, 1}, {12, 0, 1}, {14, -4, 1}};

showpic2D[cpoly, 0, 1, 0, 14, -4, 8, 0.5]
showsub2D[cpoly, 0, 1, 0, 14, -4, 8, 1]
showsub2D[cpoly, 0, 1, 0, 14, -4, 8, 4]
showsub2D[cpoly, 0, 1, 0, 14, -4, 8, 5]

(***)

cpoly = {{8, 6, 1}, {4, 6, 1}, {0, 0, 1}, {12, -4, 1}, {12, 4, 1}, {4, -4, 1}}

cpoly = {{8, 6, 1}, {4, 6, 1}, {0, 0, 1}, {12, -4, 1}, {12, 4, 1}, {4, 2, 1}}

cpoly = {{8, 6, 1}, {4, 6, 1}, {0, 0, 1}, {12, -4, 1}, {12, 4, 1}, {0, 2, 1}}

cpoly = {{8, 6, 1}, {4, 6, 1}, {0, 0, 1}, {12, -4, 1}, {12, 4, 1}, {0, 5, 1}}

cpoly = {{8, 6, 1}, {4, 6, 1}, {0, 0, 1}, {12, -4, 1}, {12, 4, 1}, {0, -4, 1}}

showpic2D[cpoly, 0, 1, 0, 12, -4, 7, 0.4]
showsub2D[cpoly, 0, 1, 0, 12, -4, 7, 5]

(***)


cpoly = {{0, 0, 1}, {6, 8, 1}, {0, 8, 1}, {6, 0, 1}, {9, 4, 1}, {-3, 4, 1}}

cpoly = {{0, 0, 1}, {6, 8, 1}, {0, 8, 1}, {6, 0, 1}, {9, 4, 1}, {0, 0, 1}}


cpoly = {{0, 0, 1}, {6, 8, 1}, {0, 8, 1}, {6, 0, 1}, {9, 4, 1}, {0, 8, 1}}

cpoly = {{0, 0, 1}, {6, 8, 1}, {0, 8, 1}, {6, 0, 1}, {9, 4, 1}, {3, 4, 1}}

cpoly = {{0, 0, 1}, {6, 8, 1}, {0, 8, 1}, {6, 0, 1}, {9, 4, 1}, {5, 3, 1}}

showpic2D[cpoly, 0, 1, -4, 9, 0, 8, 0.4]

showsub2D[cpoly, 0, 1, -4, 9, 0, 8, 4]




(*** Quartics for homework ***)

cpoly = {{6, 0, 1}, {0, 0, 1}, {6, 6, 1}, {0, 6, 1}, {6, -1, 1}}
cpoly = {{6, 0, 1}, {0, 0, 1}, {6, 10, 1}, {0, 10, 1}, {6, -1, 1}}

(*** funny ones (two double nodes) (ampersand looking) ***)
cpoly = {{6, 0, 1}, {3, 0, 1}, {8, 4, 1}, {0, 4, 1}, {6, -1, 1}}

cpoly = {{6, 0, 1}, {3, -1, 1}, {6, 6, 1}, {0, 6, 1}, {6.5, -1, 1}}

cpoly = {{6, 0}, {3, -1}, {6, 6}, {0, 6}, {7, -1}}

(*** really funny ones (one double node)  ***)

cpoly = {{0, 0, 1}, {9, 8, 1}, {5, -4, 1}, {1, 8, 1}, {10, 0, 1}}
showsub2D[cpoly, 0, 1, -1, 11, -5, 9, 4]

cpoly = {{0, -4, 1}, {10, 8, 1}, {5, -4, 1}, {0, 8, 1}, {10, -4, 1}}

(** excellent (three double points) **)
cpoly = {{0, -4, 1}, {10, 30, 1}, {5, -20, 1}, {0, 30, 1}, {10, -4, 1}}
showsub2D[cpoly, 0, 1, -1, 11, -21, 31, 4]
showsub2D[cpoly, 0, 1, 3.5, 6.5, 6, 9, 5]
showsub2D[cpoly, 0, 1, 4, 6, 6, 9, 5]

showsub2[cpoly, 0, 1, -1, 11, 4, 13, 0, 1] 
showsub2[cpoly, 0, 1, 3.5, 6.5, 6, 9, 0, 2] 
showsub2[cpoly, 0, 1, 3.5, 6.5, 6, 9, 0, 3] 
showsub2[cpoly, 0, 1, 3.5, 6.5, 6, 9, 0, 4] 
showsub2[cpoly, 0, 1, 3.5, 6.5, 6, 9, 0, 5] 
showsub2[cpoly, 0, 1, 3.5, 6.5, 6, 9, 0, 6] 


cpoly = {{0, -4, 1}, {10, 15, 1}, {5, -10, 1}, {0, 15, 1}, {10, -4, 1}}

cpoly = {{0, -4, 1}, {10, 15, 1}, {5, -12, 1}, {0, 15, 1}, {10, -4, 1}}
showsub2D[cpoly, 0, 1, -1, 11, -16, 16, 4]

(*** also very cute ***)
cpoly =  {{-1,0, 1}, {4,18, 1}, {0,-20, 1}, {-4,18, 1}, {1,0, 1}}
showsub2D[cpoly, 0, 1, -5, 5, -21, 19, 4]

inpol = {{-1,0, 1}, {6,4,1}, {-6,-12,1}, {0,20.33333,1}, 
{6,-12,1}, {-6,4,1}, {1,0,1}}
subdiv2D[inpol,0,1,0,1,-7,7,-13,21,6]
subdiv2D[inpol,0,1,0,1,-2,2,-3,3,6]

inpol = {{-1,-2,1}, {5,8,1}, {-8,-6,1}, {8,-6,1}, {-5,8,1}, {1,-2,1}}
subdiv2D[inpol,0,1,0,1,-2,2,-3,3,6]
subdiv2D[inpol,0,1,0,1,-9,9,-7,9,6]

(* parametric quartics  *)

showpic2D[-3,4,-10,10,-10,10]

(*** Examples of double nodes in rational curves  ***)
(*** degree 3 ***)
dd = {5, 10, 10, 5};
ww = {1, 1, 1, 1};
cpoly = conspoly[dd,ww,3];
showsub2D[cpoly, 0, 1, -10, 10, -10, 10, 4]

(*** degree 4 ***)
(*** next one nice ***)
dd = {5, 10, 10, 10, 5};
ww = {1, 4, 3, 4, 1};
cpoly = conspoly[dd,ww,4];
showsub2D[cpoly, 0, 1, -10, 10, -10, 10, 4]

(*** next one very nice ***)
dd = {10, 10, 10, 10, 10};
ww = {1, 4, 3, 4, 1};
cpoly = conspoly[dd,ww,4];
showsub2D[cpoly, 0, 1, -10, 10, -10, 10, 4]



(*** degree 5 ***)
dd = {10, 10, 10, 10, 10, 10};
ww = {1, 9, 15, 15, 9, 1};

(*** next one very nice ***)
dd = {10, 10, 10, 10, 10, 10};
ww = {1, 10, 17, 17, 10, 1};
cpoly = conspoly[dd,ww,5];
showsub2D[cpoly, 0, 1, -10, 10, -10, 10, 4]

cpoly =  {{-1, -2, 1}, {4, 10, 1}, {-9, -9, 1}, {9, -9, 1}, {-4, 10, 1}, {1, -2, 1}}
showsub2D[cpoly, 0, 1, -10, 10, -10, 10, 4]

(*** degree 6 ***)
dd = {10, 10, 10, 10, 10, 10, 10};
ww = {1, 9, 15, 20, 15, 9, 1};
ww = {1, 10, 40, 25, 40, 10, 1};
ww = {1, 5, 15, 20, 15, 5, 1};
ww = {1, 5, 20, 50, 20, 5, 1};
ww = {1, 9, 25, 60, 25, 9, 1};
ww = {1, 5, 15, 15, 15, 5, 1};
dd = {2, 10, 10, 10, 10, 10, 2};
ww = {1, 5, 15, 25, 15, 5, 1};

(*** next one very nice ***)
dd = {2, 10, 10, 10, 10, 10, 2};
ww = {1, 3, 15, 25, 15, 3, 1};
cpoly = conspoly[dd,ww,6];
showsub2D[cpoly, 0, 1, -4, 4, -4, 4, 5]
showsub2[cpoly, 0, 1, -3, 3, -2, 2, 0, 6]      (* great *)
showsub2D[cpoly, 0, 1, -3, 3, -2, 2, 6]     
showsub2D[cpoly, 0, 1, -11, 11, -11, 11, 4]

(*** degree 7 ***)
dd = {2, 10, 10, 10, 10, 10, 10, 2};
ww = {1, 3, 15, 25, 25, 15, 3, 1};
cpoly = conspoly[dd,ww,7];
showsub2D[cpoly, 0, 1, -11, 11, -11, 11, 4]

dd = {4, 10, 10, 10, 10, 10, 10, 4};
ww = {1, 3, 15, 20, 20, 15, 3, 1};
cpoly = conspoly[dd,ww,7];
showsub2D[cpoly, 0, 1, -11, 11, -11, 11, 4]

(*** cute ***)
dd = {3, 10, 10, 10, 10, 10, 10, 3};
ww = {1, 4, 10, 20, 20, 10, 4, 1};
cpoly = conspoly[dd,ww,7];
showsub2D[cpoly, 0, 1, -5, 5, -5, 5, 4]
showsub2D[cpoly, 0, 1, -11, 11, -11, 11, 4]

dd = {3, 10, 10, 10, 10, 10, 10, 3};
ww = {1, 3, 9, 16, 16, 9, 3, 1};
cpoly = conspoly[dd,ww,7];
showsub2D[cpoly, 0, 1, -5, 5, -5, 5, 4]

dd = {3, 10, 10, 10, 10, 10, 10, 3};
ww = {1, 4, 9, 16, 16, 9, 4, 1};
cpoly = conspoly[dd,ww,7];
showsub2D[cpoly, 0, 1, -5, 5, -5, 5, 4]




(* three leaved rose    *)

cpoly =  {{0, 0, 1}, {3/4, 0, 1}, {9/8, 3/8, 4/3}, {1, 3/4, 2}, {1/2, 1/2, 4}}

showclo2D[cpoly,0,1,-1,1,-1.2,1,6] (* great *)
rendclo2D[cpoly,0,1,-1,1,-1.2,1,1,6] (* great *)
rendclo2D[cpoly,0,1,-1,1,-1.2,1,-1,7] (* great *)


(* four leaved rose WRONG  *)
(* funny looking curve!    *)
cpoly =   {{0, 0, 1}, {2/3, 0, 1}, {10/9, 4/9, 6/5}, {1, 1, 8/5}, {4/9, 8/9, 12/5},
   {0, 0, 4}, {0, 0, 8}}

showclo2D[cpoly,0,1,-2,2,-1,9,5];
rendclo2D[cpoly,0,1,-2,2,-1,9,1,5];
rendcloz2D[cpoly,0,1,-2,2,-1,9,1,5];

(* four leaved rose *)

cpoly =    {{0, 0, 1}, {2/3, 0, 1}, {10/9, 4/9, 6/5}, {1, 1, 8/5}, {4/9, 10/9, 12/5},
   {0, 2/3, 4}, {0, 0, 8}}

showclo2D[cpoly,0,1,-1,1,-1,1,5];  (* great *)
rendclo2D[cpoly,0,1,-1,1,-1,1,1,5]; 
rendclo2D[cpoly,0,1,-1,1,-1,1,-1,7]; 

(* five leaved rose *)

cpoly =   {{0, 0, 1}, {5/6, 0, 1}, {25/18, 5/18, 6/5}, {5/4, 5/8, 8/5},
   {5/9, 5/9, 12/5}, {-1/6, 0, 4}, {-1/2, -1/2, 8}}

showclo2D[cpoly,0,1,-1,1,-1,1.1,6];  (* great *)
rendclo2D[cpoly,0,1,-1,1,-1,1.1,-1,6];  
rendclo2D[cpoly,0,1,-1,1,-1,1.1,-1,7];  

(* seven leaved rose *)


cpoly =  {{0, 0, 1}, {7/8, 0, 1}, {49/32, 7/32, 8/7}, {7/5, 21/40, 10/7},
   {35/68, 35/68, 68/35}, {-21/40, 0, 20/7}, {-35/32, -21/32, 32/7},
   {-1, -7/8, 8}, {-1/2, -1/2, 16}}


showclo2D[cpoly,0,1,-1,1,-1.1,1,6]  (* great  *) 
rendclo2D[cpoly,0,1,-1,1,-1.1,1,0,6]  (* great  *) 
rendclo2D[cpoly,0,1,-1,1,-1.1,1,-1,6]  (* great, arcs in 2 colors   *) 
rendclo2D[cpoly,0,1,-1.2,1.2,-1.1,1,-1,2]  (* great, arcs in 2 colors   *) 
rendclo2D[cpoly,0,1,-1.6,1.6,-1.1,1.6,-1,1]  (* great, arcs in 2 colors   *) 
rendcloz2D[cpoly,0,1,-1.6,1.6,-1.1,1.6,1,6]  (* great, arcs in 2 colors   *) 

(* frame w.r.t. (-1, 1)  *)
cpoly2 = {{1/2, -1/2, 16}, {8, -6, 0}, {-7/4, 7/2, 16/7}, {-8, 2, 0}, 
   {0, -35/6, 48/35}, {8, 2, 0}, {7/4, 7/2, 16/7}, {-8, -6, 0}, 
   {-1/2, -1/2, 16}}

rendcloz2D[cpoly2,0,1,-1.6,1.6,-1.1,1.6,1,6]  (* great, arcs in 2 colors   *) 
rendclo2D[cpoly2,0,1,-1.6,1.6,-1.1,1.6,-1,6]  (* great, arcs in 2 colors   *) 

(* frame w.r.t. (-2, 1)  *)

cpoly3 = {{-278/625, 556/625, 625}, {-233/1000, 13/10, -125}, 
   {2233/1600, -2933/1600, 400/7}, {28/85, -3031/680, -170/7}, 
   {-1225/548, -805/548, 548/35}, {-77/136, 91/34, -68/7}, 
   {133/64, 175/64, 64/7}, {1, 5/8, -8}, {-1/2, -1/2, 16}}


rendcloz2D[cpoly3,0,1,-1.6,1.6,-1.1,1.6,1,6]  (* great, arcs in 2 colors   *) 
rendclo2D[cpoly3,0,1,-1.6,1.6,-1.1,1.6,-1,6]  (* great, arcs in 2 colors   *) 



(* Viviani Window    *)

cpoly =   {{0, 0, 1, 1}, {1/2, 0, 1, 1}, {3/4, 1/2, 3/4, 4/3}, {1/2, 1, 1/2, 2},
   {0, 1, 0, 4}}

showclo3D[cpoly,0,1,-1,3,-6,4,-1,2,5]  (* fine  *)
rendclo3D[cpoly,0,1,-0.55,0.55,0,1,-1,1,1,5]  (* great *)
rendclo3D[cpoly,0,1,-0.55,0.55,0,1,-1,1,0,6]  (* great *)
rendclo3D[cpoly,0,1,-0.55,0.55,0,1,-1,1,-1,6]  (* great *)
rendcloz3D[cpoly,0,1,-0.56,0.56,-0.2,1.1,-1,1,-1,6]  (* great *)



(* Lissajoux type curve     *)

p1 = (1 - t^2)*(1 - 14*t^2 + t^4)
p2 = 4*t*(1 - t^2)*(1 + t^2)
p3 = (1 + t^2)^3


cpoly =   {{1, 0, 1}, {1, 2/3, 1}, {0, 10/9, 6/5}, {-5/4, 5/4, 8/5},
   {-5/3, 10/9, 12/5}, {-1, 2/3, 4}, {0, 0, 8}}

rendclo2D[cpoly,0,1,-1.1,1.1,-1.1,1.1,0,6] (* great *)
rendclo2D[cpoly,0,1,-1.1,1.1,-1.1,1.1,1,6] 
rendclo2D[cpoly,0,1,-1.1,1.1,-1.1,1.1,-1,6] 
rendcloz2D[cpoly,0,1,-1.1,1.1,-1.1,1.1,-1,6] 


(*  Limacons de Pascal  *)

poly = {{{aa - ll, 4}, {-2*aa, 2}, {aa + ll, 0}},
        {{2*(aa + ll), 1}, {2*(ll - aa), 3}},
        {{1, 4}, {2, 2}, {1, 0}}}

cpoly = gcontpoly[poly, 4, 0, 1]

glimpoly = {{aa + ll, 0, 1}, {aa + ll, (aa + ll)/2, 1},
   {(3*((2*aa)/3 + ll))/4, (3*(aa + ll))/4, 4/3},
   {ll/2, ((-aa + ll)/2 + (3*(aa + ll))/2)/2, 2},
   {0, (2*(-aa + ll) + 2*(aa + ll))/4, 4}}

(* Case aa = 1, ll = 1, cardioid *)

aa = 1; ll = 1;

limpoly =   {{2, 0, 1}, {2, 1, 1}, {5/4, 3/2, 4/3}, {1/2, 3/2, 2}, {0, 1, 4}}

rendclo2D[limpoly,0,1,-1.1,2.1,-2.1,2.1,-1,6] 


aa = 1; ll = 3/2; (* no loop *)

limpoly2 =  {{5/2, 0, 1}, {5/2, 5/4, 1}, {13/8, 15/8, 4/3}, {3/4, 2, 2}, {0, 3/2, 4}}

rendclo2D[limpoly2,0,1,-1.1,2.5,-2.5,2.5,-1,6] 
rendcloz2D[limpoly2,0,1,-1.1,2.5,-2.5,2.5,1,6] 

aa = 1; ll = 1/2;  (* loop *)

limpoly3 =  {{3/2, 0, 1}, {3/2, 3/4, 1}, {7/8, 9/8, 4/3}, {1/4, 1, 2}, {0, 1/2, 4}}

rendclo2D[limpoly3,0,1,-1.1,2,-2,2,-1,6] 
rendcloz2D[limpoly3,0,1,-1.1,2,-2,2,1,6] 

aa = 1; ll = 5/2; (* no loop *)

limpoly4 = {{7/2, 0, 1}, {7/2, 7/4, 1}, {19/8, 21/8, 4/3}, 
{5/4, 3, 2}, {0, 5/2, 4}}

rendclo2D[limpoly4,0,1,-1.5,3.5,-3.5,3.5,-1,6] 


(* WRONG, but nice anyway  *)
(* the rose r = sin 2*theta/3  *)

p1 = 4*aa*t*(1 - t^2)^2*(4*(1 - t^2)^2 - 3*(1 + t^2)^2)
p2 = 8*aa*t^2*(1 - t^2)*(3*(1 + t^2)^2 - 8*t^2)   

p1 = 4*aa*t*(1 - t^2)^2*(1 - 14*t^2 + t^4);
p2 = 8*aa*t^2*(1 - t^2)*(3 - 2*t^2 + 3*t^4);    (* -2 * t^2 should be   -10 * t^2 *)
p3 = (1 + t^2)^5;

po1 = Expand[p1]; pol1 = InputForm[po1]
po2 = Expand[p2]; pol2 = InputForm[po2]
po3 = Expand[p3]; pol3 = InputForm[po3]

pol1 = 4*aa*t - 64*aa*t^3 + 120*aa*t^5 - 64*aa*t^7 + 4*aa*t^9
pol2 = 24*aa*t^2 - 40*aa*t^4 + 40*aa*t^6 - 24*aa*t^8
pol3 = 1 + 5*t^2 + 10*t^4 + 10*t^6 + 5*t^8 + t^10

poly =  {{{4*aa, 1}, {-64*aa, 3}, {120*aa, 5}, {-64*aa, 7}, {4*aa, 9}},
         {{24*aa, 2}, {-40*aa, 4}, {40*aa, 6}, {-24*aa, 8}},
         {{1, 0}, {5, 2}, {10, 4}, {10, 6}, {5, 8}, {1, 10}}};



cpoly = gcontpoly[poly, 10, 0, 1]


cpoly =   {{0, 0, 1}, {(2*aa)/5, 0, 1}, {(18*aa)/25, (12*aa)/25, 10/9},
   {aa/2, (6*aa)/5, 4/3}, {(-14*aa)/45, (79*aa)/45, 12/7},
   {(-45*aa)/37, (69*aa)/37, 148/63}, {(-71*aa)/45, (14*aa)/9, 24/7},
   {(-6*aa)/5, (11*aa)/10, 16/3}, {(-12*aa)/25, (18*aa)/25, 80/9},
   {0, (2*aa)/5, 16}, {0, 0, 32}}


(* Case aa = 1 *)
 

cpoly =  {{0, 0, 1}, {2/5, 0, 1}, {18/25, 12/25, 10/9}, {1/2, 6/5, 4/3},
   {-14/45, 79/45, 12/7}, {-45/37, 69/37, 148/63}, {-71/45, 14/9, 24/7},
   {-6/5, 11/10, 16/3}, {-12/25, 18/25, 80/9}, {0, 2/5, 16}, {0, 0, 32}}



rendclo2D[cpoly,0,1,-2.1,2.1,-2.1,2.1,-1,5]  (* cute *)
rendcloz2D[cpoly,0,1,-2.1,2.1,-2.1,2.1,1,5]  (* cute *)



(* The right rose  r = sin 2*theta/3  *)

p1 = 4*aa*t*(1 - t^2)^2*(1 - 14*t^2 + t^4);
p2 = 8*aa*t^2*(1 - t^2)*(3 - 10*t^2 + 3*t^4);  
p3 = (1 + t^2)^5;

po1 = Expand[p1]; pol1 = InputForm[po1]
po2 = Expand[p2]; pol2 = InputForm[po2]
po3 = Expand[p3]; pol3 = InputForm[po3]


pol1 = 4*aa*t - 64*aa*t^3 + 120*aa*t^5 - 64*aa*t^7 + 4*aa*t^9
pol2 = 24*aa*t^2 - 104*aa*t^4 + 104*aa*t^6 - 24*aa*t^8
pol3 = 1 + 5*t^2 + 10*t^4 + 10*t^6 + 5*t^8 + t^10

poly =  {{{4*aa, 1}, {-64*aa, 3}, {120*aa, 5}, {-64*aa, 7}, {4*aa, 9}},
         {{24*aa, 2}, {-104*aa, 4}, {104*aa, 6}, {-24*aa, 8}},
         {{1, 0}, {5, 2}, {10, 4}, {10, 6}, {5, 8}, {1, 10}}};

cpoly = gcontpoly[poly, 10, 0, 1]

cpoly =  {{0, 0, 1}, {(2*aa)/5, 0, 1}, {(18*aa)/25, (12*aa)/25, 10/9},
   {aa/2, (6*aa)/5, 4/3}, {(-14*aa)/45, (71*aa)/45, 12/7},
   {(-45*aa)/37, (45*aa)/37, 148/63}, {(-71*aa)/45, (14*aa)/45, 24/7},
   {(-6*aa)/5, -aa/2, 16/3}, {(-12*aa)/25, (-18*aa)/25, 80/9},
   {0, (-2*aa)/5, 16}, {0, 0, 32}};


cpoly = fcontpoly[poly, 10, 0, 1]
cpoly = yfcontpoly[poly, 10, 0, 1]

cpoly = {{0, 0, 1}, {(2*aa)/5, 0, 1}, {(18*aa)/25, (12*aa)/25, 10/9},
   {aa/2, (6*aa)/5, 4/3}, {(-14*aa)/45, (71*aa)/45, 12/7},
   {(-45*aa)/37, (45*aa)/37, 148/63}, {(-71*aa)/45, (14*aa)/45, 24/7},
   {(-6*aa)/5, -aa/2, 16/3}, {(-12*aa)/25, (-18*aa)/25, 80/9},
   {0, (-2*aa)/5, 16}, {0, 0, 32}}


(* Case aa = 1 *)
 

rcpoly =  {{0, 0, 1}, {2/5, 0, 1}, {18/25, 12/25, 10/9}, {1/2, 6/5, 4/3},
   {-14/45, 71/45, 12/7}, {-45/37, 45/37, 148/63}, {-71/45, 14/45, 24/7},
   {-6/5, -1/2, 16/3}, {-12/25, -18/25, 80/9}, {0, -2/5, 16}, {0, 0, 32}};


rendclo2D[rcpoly,0,1,-1.1,1.1,-1.1,1.1,-1,6]  (* great *)
rendcloz2D[rcpoly,0,1,-1.1,1.1,-1.1,1.1,1,6]  (* great *)

rendclo2D[rcpoly,0,1,-1.1,1.1,-1.1,1.1,-1,1] 

rendclo2D[rcpoly,0,1,-1.1,1.1,-1.1,1.1,-1,2] 

rendclo2D[rcpoly,0,1,-1.1,1.1,-1.1,1.1,-1,3] 

(* another ellipse *)

p1 = 4*t;
p2 = t^2 -3*t + 2;
p3 = t^2 + 1;

poly = {{{4, 1}},
        {{1, 2}, {-3, 1}, {2, 0}},
        {{1, 2}, {1, 0}}};
 
cpoly = gcontpoly[poly, 2, -1, 1]

cpoly = {{-2, 3, 2}, {0, 1, 0}, {2, 0, 2}}

cpoly1 = gcontpoly[poly, 2, 0, 1]

cpoly1 = {{0, 2, 1}, {2, 1/2, 1}, {2, 0, 2}}

rendclo2D[cpoly1,0,1,-2,2.1,-1.1,3.5,0,6] 
rendclo2D[cpoly1,0,1,-2,2.1,-1.1,3.5,-1,6] 


(* Tricuspoid of Steiner *)
(* x = a(2 Cos t + Cos 2t), y = a(2 Sin t - Sin 2t)  *)

x = a(3 - 6t^2 - t^4)
y = 8at^3
z = 1 + 2t^2 + t^4

trpoly = {{{3*aa, 0}, {-6*aa, 2}, {-aa, 4}},
          {{8*aa, 3}},
          {{1, 0}, {2, 2}, {1, 4}}};

cspoly = fcontpoly[trpoly, 4, 0, 1];

cspoly =  {{3*aa, 0, 1}, {3*aa, 0, 1}, {(3*aa)/2, 0, 4/3}, {0, aa, 2}, {-aa, 2*aa, 4}};

aa = 2;

cspoly =  {{6, 0, 1}, {6, 0, 1}, {3, 0, 4/3}, {0, 2, 2}, {-2, 4, 4}};

rendclo2D[cspoly,0,1,-3,6,-6,6,-1,1] 
rendclo2D[cspoly,0,1,-3,6,-6,6,-1,2] 
rendclo2D[cspoly,0,1,-3,6,-6,6,-1,3] 

rendclo2D[cspoly,0,1,-3,6,-6,6,-1,6] 
rendcloz2D[cspoly,0,1,-3,6,-6,6,1,6] 
rendclo2D[cspoly,0,1,-3,6,-6,6,1,3] 


(* Astroid *)
(* x = a Cos^3 t, y = a Sin^3 t  *)

x = a(1 - 3t^2 + 3t^4 - t^6)
y = 8at^3
z = 1 + 3t^2 + 3t^4 + t^6

aspoly = {{{aa, 0}, {-3*aa, 2}, {3*aa, 4}, {-aa, 6}},
          {{8*aa, 3}},
          {{1, 0}, {3, 2}, {3, 4}, {1, 6}}};

ascpoly = fcontpoly[aspoly, 6, 0, 1];

ascpoly =  {{aa, 0, 1}, {aa, 0, 1}, {(2*aa)/3, 0, 6/5}, {aa/4, aa/4, 8/5},
   {0, (2*aa)/3, 12/5}, {0, aa, 4}, {0, aa, 8}};

aa = 2;

ascpoly = {{2, 0, 1}, {2, 0, 1}, {4/3, 0, 6/5}, {1/2, 1/2, 8/5}, {0, 4/3, 12/5},
   {0, 2, 4}, {0, 2, 8}};

rendclo2D[ascpoly,0,1,-2,2,-2.1,2.1,-1,1]
rendclo2D[ascpoly,0,1,-2,2,-2.1,2.1,-1,2] 
rendclo2D[ascpoly,0,1,-2,2,-2.1,2.1,-1,4] 
rendclo2D[ascpoly,0,1,-2,2,-2.1,2.1,1,4]


(* Bicorne *)
(* y^2(a^2 - x^2) =   (x^2 + 2ay - a)^2 *)
(* x = a Cos t, y = -a Sin^2 t/(Sin t -2) *)

x = a(-t^4 + t^3 - t + 1)
y = 2at^2
z = t^4 - t^3 + 2t^2 - t + 1


bipoly = {{{aa, 0}, {-aa, 1}, {aa, 3}, {-aa, 4}},
          {{2*aa, 2}},
          {{1, 0}, {-1, 1}, {2, 2}, {-1, 3}, {1, 4}}};

bicpoly = fcontpoly[bipoly, 4, 0, 1];


bicpoly =  {{aa, 0, 1}, {aa, 0, 3/4}, {(3*aa)/5, (2*aa)/5, 5/6}, {aa/2, aa, 1},
   {0, aa, 2}}

aa = 2;

bicpoly = {{2, 0, 1}, {2, 0, 3/4}, {6/5, 4/5, 5/6}, {1, 2, 1}, {0, 2, 2}};

rendclo2D[bicpoly,0,1,-2.1,2.1,-0.1,2,-1,4]
rendclo2D[bicpoly,0,1,-2.1,2.1,-0.1,2,-1,1]
rendclo2D[bicpoly,0,1,-2.1,2.1,-0.1,2,1,1]
rendclo2D[bicpoly,0,1,-2.1,2.1,-0.1,2,1,2]
rendclo2D[bicpoly,0,1,-2.1,2.1,-0.1,2,1,4]
rendcloz2D[bicpoly,0,1,-2.1,2.1,-0.1,2,1,5]


(* Freeth's Nephroid *)
(* r = a(1 + 2 Sin theta/2) *)


exp1 = aa*(t^2 + 4*t + 1)*(t^4 - 6*t^2 + 1);
exp2 = 4*aa*t*(t^2 + 4*t + 1)*(1 - t^2);
exp3 = t^6 + 3*t^4 + 3*t^2 + 1;

e1 = Expand[exp1];  InputForm[e1];
e2 = Expand[exp2];  InputForm[e2];
e3 = exp3;

e1 = aa + 4*aa*t - 5*aa*t^2 - 24*aa*t^3 - 5*aa*t^4 + 4*aa*t^5 + aa*t^6
e2 = 4*aa*t + 16*aa*t^2 - 16*aa*t^4 - 4*aa*t^5
e3 = exp3

npoly = {{{aa, 0}, {4*aa, 1}, {-5*aa, 2}, {-24*aa, 3}, {-5*aa, 4}, {4*aa, 5}, {aa, 6}},
          {{4*aa, 1}, {16*aa, 2}, {-16*aa, 4}, {-4*aa, 5}},
          {{1, 0}, {3, 2}, {3, 4}, {1, 6}}};

ncpoly = fcontpoly[npoly, 6, 0, 1];

ncpoly =   {{aa, 0, 1}, {(5*aa)/3, (2*aa)/3, 1}, {(5*aa)/3, 2*aa, 6/5},
   {aa/2, (13*aa)/4, 8/5}, {(-13*aa)/9, (10*aa)/3, 12/5}, {-3*aa, 2*aa, 4},
   {-3*aa, 0, 8}};

aa = 2;

ncpoly =   {{2, 0, 1}, {10/3, 4/3, 1}, {10/3, 4, 6/5}, {1, 13/2, 8/5},
   {-26/9, 20/3, 12/5}, {-6, 4, 4}, {-6, 0, 8}};

rendclo2D[ncpoly,0,1,-6.1,3.1,-6.1,6.1,-1,4]
rendclo2D[ncpoly,0,1,-6.1,3.1,-6.1,6.1,1,5]
rendclo2D[ncpoly,0,1,-6.1,3.1,-6.1,6.1,-1,6]


(* Trisectrix of Maclaurin *)
(* y^2(a + x) = x^2(3a - x) *)
(* x = a(3 - m^2)/(m^2  + 1), y = at(3 - m^2)/(m^2 + 1) *)

mcpoly = {{{3*aa, 0}, {-aa, 2}},
          {{3*aa, 1}, {-aa, 3}},
          {{1, 0}, {1, 2}}};

mlcpoly = fcontpoly[mcpoly, 3, 0, 1];

mlcpoly = {{3*aa, 0, 1}, {3*aa, aa, 1}, {2*aa, (3*aa)/2, 4/3}, {aa, aa, 2}};

aa = 2;

mlcpoly = {{6, 0, 1}, {6, 2, 1}, {4, 3, 4/3}, {2, 2, 2}};

rendclo2D[mlcpoly,0,1,-2.1,6.1,-8.1,8.1,1,5]


ClearAll["Global`*"]