(* ::Package:: *)

(* ------------------------------------------------------------------------ *)
(* MinimalSurfaces.m Header                                                 *)
(* ------------------------------------------------------------------------ *)

(* :Title: Minimal Surfaces *)

(* :Context: MinimalSurfaces` *)

(* :Author: O. Michael Melko *)

(* :Copyright: Copyright 2005, O. Micheal Melko *)

(* :Package Version: 0.1 *)

(* :Mathematica Version: 7.0 *)

(* :Summary:
	This package provides tools for rendering minimal surfaces graphically,
	or by means of stereolythography.
*)

(* :History:
	Revised 2 December 2005
	Revised 5 January 2010
*)

(* :Keywords:
	Differential geometry; Riemann surfaces;
	minimal surfaces; Costa's surface; Delaunay triangulation;
	stereolythography; 3D printing.
 *)

(* :Sources: *)

(* ------------------------------------------------------------------------ *)
(* Begin Package                                                            *)
(* ------------------------------------------------------------------------ *)

BeginPackage["MinimalSurfaces`"]
Needs["ComputationalGeometry`"]

(* ------------------------------------------------------------------------ *)
(* Function Usage / Global Context Definition                               *)
(* ------------------------------------------------------------------------ *)

CounterClockwise::usage = ""
Clockwise::usage = ""

Arc::usage =
   "Arc[{x, y}, r, {th1, th2}, orientation] is a two-dimensional graphics
	primitive specifying a circular arc with center {x, y}, radius r, and 
	initial and terminal angles {th1, th2} which are assumed to satisfy
	-Pi <= th1, th2 <= Pi. The orientation may be either Clockwise or 
	CounterClockwise."

CostaSurface::usage =
   "CostaSurface[u, v] returns a list of three complex numbers. The real
	part of this list specifies a point on Costa's minimal surface. The real
	numbers u and v are assumed not to be half-integers."

ParallelSurface::usage =
   "ParallelSurface[X_List, u_List, d_] specifes the parameterization of a
	surface obtained from the parameterization X via a normal displacement
	of distance |d|. Orientation of the normal is assgined to X by means of
	the right-hand rule. When d > 0, displacement is in direction of this
	normal. Otherwise, it's in the opposite direction."

CreateVertexData::usage =
   "CreateVertexData[c_List, ops___Rule]"

	MeshSize::usage =
   "MeshSize -> {n1, n2} is an option of CreateVertexData that specifies the
	number of sample points to use in the x and y directions. The defaults are
	n1 = 10, and n2 =10."

	Options[CreateVertexData] = {MeshSize -> {10, 10}}

GlueComponents::usage =
   "GlueComponents[top_List, bot_List] glues together a pair polyhedral 
	surfaces specified by top and bot. Each list has the form
	{vertices, faces, boundary}, where vertices is a list of points in R^3.
	The faces list give incidence relations in vertices, which define polygons, 
	and boudary is a list of sublists, each pointing to vertices that lie on a 
	component of the boundary of a surface. The boundary components of top and 
	bottom are assumed to be of the same shape. The output specifies the
	polyhedral boundary of a solid in R^3. It is in the form of
	{vertices, faces}. Here, vertices is a concatenation of the corresponding
	lists in top and bottom, and faces = {f1,f2, f3}. The lists f1 and f2
	are the incidence relations for the two surfaces, and f3 is a list whose
	sublists specify connections between the boundary components."

Triangulate::usage =
   "Triangulate[v_List, c_List, ops___Rule] creates a triangulation of a
	Riemann surface."

	Identifications::usage =
   "Specifies boundary indentifications of planar region."
	
	Options[Triangulate] = {Identifications -> {}}

CreatePolyhedron::usage =
   "CreatePolyhedron"

ExportGraphics::usage =
   "ExportGraphics[file, data, format] exports data to file using format"

Begin["`Private`"]

(* ------------------------------------------------------------------------ *)
(* Error messages for exported objects                                      *)
(* ------------------------------------------------------------------------ *)

MinimalSurfaces::badarg = "bad argument"

(* ------------------------------------------------------------------------ *)
(* Constants & Data Types                                                   *)
(* ------------------------------------------------------------------------ *)

Clockwise := 1

CounterClockwise := -1

Arc[c, r, a, o]

(* ------------------------------------------------------------------------ *)
(* Surface Generating Functions                                             *)
(* ------------------------------------------------------------------------ *)

(*	Func:	CostaSurface

	Pre:	u + I v is not a pole of the "WeierstrassP" function.

	Post:	Returns complex parameterization for Costa Surface.
*)

CostaSurface[u_, v_]:=
Module[{x, y, a, b, e, g, p, z},

	g = Chop[N[WeierstrassInvariants[{1/2, I/2}]]];

	z[x_, y_] := WeierstrassZeta[x + I y, g];
	p[x_, y_] := WeierstrassP[x + I y, g];

	e1 = N[p[1/2, 0]];	(* The "first root" *)
	a  = Pi^2/(4 e1);
	b  = Pi/(2 e1);

	phi = {
		  -z[u, v] + Pi u + a +   b (z[u - 1/2, v] - z[u, v - 1/2]),
		-I z[u, v] + Pi v + a - I b (z[u - 1/2, v] - z[u, v - 1/2]),
		   Sqrt[Pi/2] Log[(p[u, v] - e1)/(p[u, v] + e1)]
	};
	
	phi/2
]

(*	Func:	ParallelSurface

	Desc:	Create Parallel Surface of distance |d| to parameterized surface
		given by X. Assign orientation of normal by right-hand rule.

	Pre:	X	Real or complex parametrization.

			u	Point in the domain of X.

			d	Real number.

	Post:	Returns real parameterization for surface parallel to X.
*)

ParallelSurface[X_List, u_List, d_]:=
Module[{Xu, Xv, N},

	Xu = Re[D[X,u[[1]]]];
	Xv = Re[D[X,u[[2]]]];

	N = {	Xu[[2]] Xv[[3]] - Xu[[3]] Xv[[2]],
			Xu[[3]] Xv[[1]] - Xu[[1]] Xv[[3]],
			Xu[[1]] Xv[[2]] - Xu[[2]] Xv[[1]]
		};
	N = N/norm[N];

	Re[X] + d N
]

(* ------------------------------------------------------------------------ *)
(* Vertex Creation and Triangulation                                        *)
(* ------------------------------------------------------------------------ *)

(*	Func:	CreateVertexData

	Pre:	c	specifies the boundary of a planar region. It consists of a 
				list of line segments and circular arcs traversed
				counterclockwise.

	Post:	Returns a list {interior, boundary} of two sublists of vertices
			within the boundary specified by b:

				interior	List of interior vertices.

				boundary	List of lists. Each sublist contains vertices
							for a line segment or arc from the boundary "b" 
							of the region to be triangulated.

	Rem:	ops may specify method of generating points on boundary
			elements "c", and in interior of the domain. Defaults might be:

				Distribution -> Uniform
				MeshSize -> {n1, n2}

			where n1, n2 specifies the number of mesh points in x, y	
			directions. An option could be provided to exclude subregions 
			from the domain bounded by "c1". It might read as follows:

				Exclude -> c1

			where "c1" is a list similar to "c".
*)

CreateVertexData[c_List, ops___Rule]:=
Module[
	{	k = 0,
		dX, dY, n = {}, q1 = 0, q2 = 0,
		m = {}, boundary = {}, interior = {}
	},
	
	n = MeshSize/.{ops}/.Options[CreateVertexData];

	m = boundingRectangle[c];

	dX = (m[[2, 1]] - m[[1, 1]]) / n[[1]];  (* horizontal stepsize *)
	dY = (m[[2, 2]] - m[[1, 2]]) / n[[2]];  (* vertical stepsize *)

	v = Flatten[
			Table[m[[1]] + {dX * i, dY * j}, {i, 0, n[[1]]}, {j, 0, n[[2]]}],
			1
		];

	interior = Select[v, insideQ[c, #]&];


	While[k++ < Length[c],
		Switch[Head[ c[[k]] ],
			Line,
			boundary = Append[
							boundary,
							partitionLine[c[[k]], {q1, q2}, {dX, dY}]
						],

			Arc,
			boundary = Append[
							boundary,
							partitionArc[c[[k]], {q1, q2}, {dX, dY}]
						]
		];
	];

	{interior, boundary}

]

(*	Func:	GlueComponents

	Desc:	Glues together a pair polyhedral surfaces obtained by normal
			displacement of a given parameterization. The result specifies
			the polyhedral boundary of a solid in R^3.

	Pre:	The lists "top" and "bot" are in the same form as output from
			the Triangulate function

			top = {verts1, faces1, boundary1}
			bot = {verts2, faces2, boundary2}

	Post:	Returns {w, f}, where w = Join[verts1, verts2]. The list f is of
			the form f = {faces1, faces2, f3}, where f3 is a list whose
			sublists consist of incidence relations between the corresponding
			boundary components of "top" and "bot". Hence f3 "glues" "top and
		   "bot " together. Face orientations are chosen in a consistent
			manner.
*)

GlueComponents[top_List, bot_List]:=
Module[{n = Length[top[[1]]], w = {}, f={}, b = {}},

	swap[l_] := l /. {l[[2]] -> l[[3]], l[[3]] -> l[[2]]};

(*	the tr function creates a "regular triangulation" of two flat lists of the
	same length:
*)
	tr[l1_List, l2_List] :=
		Module[{n = Length[l1], t1 = {}, t2 = {}},
    	t1 = Table[{l1[[k]], l2[[k]], l1[[ Mod[k + 1, n, 1] ]]}, {k, 1, n}];
    	t2 = Table[
				{l2[[k]], l2[[Mod[k + 1, n, 1]]], l1[[ Mod[k + 1, n, 1] ]]},
				{k, 1, n}
			];
    	Join[t1, t2]
    ];

(*	w is a flat list of vertices *)
	w = Join[ top[[1]], bot[[1]] ];

(*	adjust incidence references in b[[3]] by n; triangulate corresponding
	components
*)
	b = MapThread[tr, {top[[3]], bot[[3]] + n}];

(*	reverse orientation of triangles in bot *)
	f = {top[[2]], Map[ swap, bot[[2]] ] + n, b};

	{w, f}
]

(*	Func:	Triangulate

	Desc:	Produces triangulation of a Riemann surface with disks about
			marked points removed.

	Pre:	v	List of two sublists in the same form as output from
				"CreateVertexData". Specifically:

					v[[1]]	List of interior vertices.

					v[[2]]	List of lists. Each sublist contains vertices for
							a line segment or arc from the boundary "c" of the
							region to be triangulated.

			c	specifies the boundary of a planar region; it consists of a
				list of line segments and circular arcs traversed
				counterclockwise.

			ops Specifies boundary identifications of sublists in v[[2]].
				The syntax is as follows:

					Identifications -> {}	(default, no identifications)

					Identifications -> {l1, ,l2}

				The sublist l1 specifies edges in v[[2]] that are to be
				identified, and the sublist l2 specifies the closed curves
				comprising the new boundary, For example:

					l1 = {{1, 11}, {3, 9}, {5, 15}, {7, 13}}
					l2 = {{2, 10}, {4, 16, 12, 8}, {6, 14}}

	Post:	Returns vertex data and incidence relations of a triangulation for
			a Riemann surface with "marked disks" removed. Returns incidence
			list of boundary vertices.

			Ouput:

			List "w" and list "f" of incidence relations, i.e., a list
			whose elements are lists of three numbers, such as
				f[[k]]  = {a, b, c}
			where a, b, c "point" to vertices
				v[[a]], v[[b]], v[[c]]

			List cp = {cp1, cp2, ... } is a list of lists. Each element cpj
			gives the position of vertices for a connected component of the
			boundary of the triangulated region.

	Rem:	The adjacency sublists "a" returned by "DelauneyTriangulation are
			not cyclic. Hence, if a[[1]] is an interior point, it is necessary
			to add the last face incident a[[1]] before deleting vertices from
			the triangulation.
*)

Triangulate[v_List, c_List, ops___Rule]:=
Module[
	{	d, i, m, n, eps = 10^(-4), r, x, y, z, gamma,
		catenate, common, cyclicPerm, pos, shift,
		b = {}, b1 = {}, b2 = {}, b3 = {}, cp = {}, l = {}, t = {}, tc = {},
		ids = {}, w = {}, wR = {}, f = {}
	},

	ids = Identifications/.{ops}/.Options[Triangulate];
	If[ ids != {}, ids[[1]] = Transpose[ ids[[1]] ] ];

(*	if ids != {}:
		Reverse direction of sublists of vertices to be discarded; put them
		at the end; return flattened list. The "Select" statement defining
		"b3" removes any duplicate endpoints that are not discarded later.
	else:
		remove duplicate elements; return flattened list
*)

	catenate[z_] := Apply[Join, Map[Drop[#, -1]&, z]];
	shift[x_, y_, z_] :=
		If[	MemberQ[z, x],
			y[[ Part[Flatten[Position[z, x]], 1] ]],
			x
		];

	If[ids != {},
		b1 = Flatten[ Part[ v[[2]], ids[[1, 1]] ], 1];
		b2 = Part[ v[[2]], ids[[1, 2]] ];
		b2 = Flatten[Map[Reverse, b2], 1];
		b3 = Select[
				Flatten[v[[2]], 1],
				(MemberQ[Join[b1, b2], #] == False)&
			 ];
		b  = Join[b3, b1, b2];

		gamma = Map[ catenate, Map[ Extract[v[[2]], #]&, ids[[2]], {2} ] ];
		gamma = Map[ shift[#, b1, b2]&, gamma, {2} ],

		b  = catenate[ v[[2]] ]
	];

(* 	extract marked points from c *)

	p = Map[Part[#, 1]&, Select[c, Head[#] == Arc&]];

(*	augment vertex data with p, randomize  and triangulate *)

	SeedRandom[Prime[100]];
	r[z_] := z + eps*(Random[]-1/2);

	w = Join[v[[1]], b, p];
	wR = Map[r, w];

	t = DelaunayTriangulation[wR];

(*	insert cyclic adjacency for t[[i, 1]] if it is interior *)

	For[i = 1, i <= Length[t], i++,
    	If[ MemberQ[t[[t[[i, 2, 1]], 2]], t[[i, 2, -1]]],
        	t[[i, 2]] = Append[t[[i, 2]], t[[i, 2, 1]]]
		];
    ];

(*	build list of faces f from t; remove redundancies *)

	cyclicPerm[x_] := NestList[RotateLeft, x, 2];
	common[x_, y_] := Length[Intersection[x, y]];

	While[Length[t] > 0,
		tc = Map[Prepend[#, t[[1, 1]]] &, Partition[t[[1, 2]], 2, 1]];
		f = Join[f, Select[tc, (common[f, cyclicPerm[#]] == 0)&]];
		t = Drop[t, 1];
	];

(* remove triangles incident to marked points *)

	w = Drop[w, -Length[p]];
	f = Select[f, (Max[#] <= Length[w])&];

(* remove degenerate triangles resulting from vertex randomization *)

	d[z_] := Det[ { z[[2]] - z[[1]], z[[3]] - z[[1]] } ];
	f = Select[f, d[Part[w, #]] != 0&];

(*	Vertices to be discarded are at end of "w". Let Vp, resp. Vq be a
	vertex with index p, resp. q. If Vp~Vq, replace all occurences of q in
	f with p.
*)

	pos[y_, z_] := Flatten[Map[Position[y, #]&, z]];

	If[ids != {},
		m = Length[b2];
		n = Length[w];
		l = Range[n - m + 1, n];

		f = ReplaceAll[f, MapThread[Rule, {l, l - m}]];
		w = Drop[w, -m];
		cp = Map[pos[w, #]&, gamma],

		cp = Map[pos[w, #]&, {Take[w, -Length[b]]}]
	];

	{w, f, cp}
]

(* ------------------------------------------------------------------------ *)
(* Graphics Functions                                                       *)
(* ------------------------------------------------------------------------ *)

(*	Func:	CreatePolyhedron

	Pre:	p	List of the form p = {w, x}, where x = {x1, x2, ... }, and
				xj = {gj, fj}. Each fj is a list of faces in the polyhedron,
				defined by incidence relations to w, and each gj is a list of
				graphics directives to be applied to fj.

	Post:	Returns a list of 3D graphics directives and primitives.

	Rem:	This function could easily be modified to apply directives to
			edges and vertices.
*)

CreatePolyhedron[p_List]:=
Module[{f, g = {}, l, q = {}, t},

(*	q = {{g1, g2, ... }, {f1, f2, ...}} *)

	If[ Length[p[[2]]] > 1,
		q = Transpose[p[[2]]], q = {{p[[2, 1, 1]]}, {p[[2, 1, 2]]}}
	];

(*	replace vertex list positions with vertex coordinate data *)

	q[[2]] = Map[Part[p[[1]], #]&, q[[2]], {2}];

(*	{fj1, fj2, ... } -> {Polygon[fj1], Polygon[fj2], ... } *)

	q[[2]] = Map[Polygon, q[[2]], {2}];

	If[ Length[p[[2]]] > 1, Transpose[q], {Flatten[q, 1]} ]
]

(* ------------------------------------------------------------------------ *)
(* Import/Export Functions                                                  *)
(* ------------------------------------------------------------------------ *)

(*	Func:	ExportGraphics

	Pre:	file	Path to output file.

			data	Geometric data describing the surface of
					a solid object.
						A.	Include incidence list from modified  Delauney.
						B.	Include output from "CreatePolyhedron"?

			style	Export format.

	Post:	Export data to one of following formats:
				PLY
				POVRAY
*)

ExportGraphics[file_String, data_List, style_String]:=
Module[{},

	Switch[style,
		"PLY", createPLY[file, data],
		"POVRAY", createPOVRAY[file, data],
		_, createPLY[file, data]
	]
]

(* ------------------------------------------------------------------------ *)
(* Utility Functions                                                        *)
(* ------------------------------------------------------------------------ *)

(*	Func:	boundingRectangle

	Pre:	b	Specifies the boundary of a planar region. It consists of a
				list of line segments and circular arcs, traversed
				counterclockwise.

	Post:	Returns smallest rectangle containing "b".
*)

boundingRectangle[b_List]:=
Module[{l,l1, l2, xMin, xMax, yMin, yMax},

(*	extract vertex data from b *)

	l1 = Map[Part[#,1]&, Select[b,(Head[#] == Arc)&]];
	l2 = Flatten[Map[Part[#,1]&, Select[b,(Head[#] == Line)&]], 1];

	l = Join[l1, l2];

(*	find coordinates of bounding rectangle *)

	l = Sort[l, ( #1[[1]] < #2[[1]] )&];

	xMin = l[[ 1, 1]];
	xMax = l[[-1, 1]];

	l = Sort[l, ( #1[[2]] < #2[[2]] )&];

	yMin = l[[ 1, 2]];
	yMax = l[[-1, 2]];

(*	return bounding rectangle *)

	{{xMin, yMin}, {xMax, yMax}}
]


(*	Func:	createPLY

	Pre:	p	List of the form p = {w, x}, where x = {x1, x2, ... }, and
				xj = {gj, fj}. Each fj is a list of faces in the polyhedron,
				defined by incidence relations to w, and each gj is a list of
				graphics directives to be applied to fj.

	Post:	Creates a file in PLY format:

			header text
			v (1,1) v (1,2) v (1,3)
			...
			v (r,1) v (r,2) v (r,3)
			r1 g1 b1 n1 f (1,1) ... f (1,n1)
			...
			rs gs bs ns f (s,1) ... f (s,ns)
			
			where v (i,j) are components of vertices, and each line thereafter
			specificies a face fj with nj vertices f (j,k), 1 <= k <= nj, and
			color given by rj, bj, bj.

	Rem:	This function could easily be modified to apply color attributes
			to edges and vertices. The file format has been tested on ZCorp's
			Model 402Z 3D printer.
*)

createPLY[file_String, p_List]:=
Module[{c, i, j, k, l = {}, m, n, f = {}, v = {}, q, s, x},

(*	q puts each sublist in p[[2]] in proper form for export. Note that
	in PLY format, indexing starts at 0.
*)
	q[l_] := Map[ Join[ {Length[#]}, # - 1, l[[1]] ] &, l[[2]] ];

	f = Map[q, p[[2]]];
	v = Chop[p[[1]]];

	m = Length[v];
	n = Apply[Plus, Map[Length, f]];
	s = OpenWrite[file];

(*	File Header	*)

	WriteString[s,
		"ply\n",
		"format ascii 1.0\n",
		"comment created by MinimalSurfaces.m, Version 0.1\n",
		"comment Copyright 2005, O. Michael Melko.\n",
		"comment solid model coordinates for a surface.\n",

		"element vertex " <> ToString[m] <> "\n",
		"property float x\n",
		"property float y\n",
		"property float z\n",

		"element face " <> ToString[n] <> "\n",
		"property list uchar int vertex_indices\n",
		"property uchar red\n",
		"property uchar green\n",
		"property uchar blue\n",

		"end_header\n"
	];

(*	write vertex data to file *)

	For[j = 1, j <= m, j++,
		c = " ";
		For[k = 1, k <= Length[v[[j]]], k++,
			If[k == Length[v[[j]]], c = "\n"];
			WriteString[s, v[[j, k]], c];
		];
	];

(*	write face data to file *)

	For[i = 1, i <= Length[f], i++,
		For[j = 1, j <= Length[f[[i]]], j++,
			c = " ";
			For[k = 1, k <= Length[f[[i, j]]], k++,
				If[k == Length[f[[i, j]]], c = "\n"];
				WriteString[s, f[[i, j, k]], c];
			];
		];
	];

	Close[s];
]

(*	Func:	createPOVRAY

	Pre:

	Post:
*)

createPOVRAY[file_String, data_List]:=
Module[{},

	s = OpenWrite[file];

	(* someday, put code here *)

	Close[s]
]

(*	Func:	insideQ

	Desc:	Tests if a point is inside a region of the plane.

	Pre:	b	Specifies the boundary of a planar region. It consists of a
				list of line segments and circular arcs, traversed
				counterclockwise.

			p	Point in the plane.

	Post:	Returns "True" if "p" is in the region bounded by "b",
			otherwise returns "False".

	Rem:	Currently the "Arc" part only excludes points from a disk.
			This could be generalized by testing the orientation of the arc.
	
			There seems to be a Mathematica bug here. The last condition
			in the Switch statement be should read "_Arc", not "_". When this
			change is made, "insideQ" fails to evaluate that condition. It
			works, however, if "insideQ" with the change is copied into a
			notebook and evaluated there.
*)

insideQ[b_List, p_List]:=
Module[{k = 0, test = 1},

	While[(test == 1) && (k++ < Length[b]),
		Switch[b[[k]],
			_Line,
			If[ Det[{ b[[k, 1, 2]] - b[[k, 1, 1]],
					p - b[[k, 1, 1]] }] > 0,
				test = 1, test = -1
			],

			_,
			If[ norm[ p - b[[k, 1]] ] > b[[k, 2]],
				test = 1, test = -1
			]
		]
	];

	If[test == 1, True, False]
]

(*	Func:	partitionArc

	Pre:	a	Specifies a circular arc. It is assumed to be an expression of
				the from "Arc[{c1, c2}, r, {th1, th2}, or]", where
					c1		coordinates of center of a circular arc.
					c2
					r		radius,
					th1		initial angle with respect to the
							horizontal axis,
					th2		terminal angle with respect to the
							horizontal axis.
								-Pi <= th1, th2 <= Pi
					or		orientation.

			q	List containing the horizontal and vertical offset of a mesh.

			s	List containing the horizontal and vertical stepsize of a
				mesh.

	Post:
*)

partitionArc[a_Arc, q_List, s_List]:=
Module[
	{	c = a[[1]], r = a[[2]], th = a[[3]], or = a[[4]],
		x1, x2, y1, y2,
		l = {}, l1 = {}, l2 = {},
		lX = {}, lPX = {}, lMX = {}, lY = {}, lPY = {}, lMY = {},
		z = {}
	},

(*	compute Arg of z - w in range (-Pi, Pi] *)

	ar[z_, w_] := Arg[ (z[[1]] - w[[1]]) + I (z[[2]] - w[[2]]) ];

(*	partition horizontal projection of circle *)

	x1 = c[[1]] - r;
	x2 = c[[1]] + r;

	lX = partitionInterval[{x1, x2}, q[[1]], s[[1]]];

(*	drop endpoints if they are not on the mesh *)

	If[ FractionalPart[ (lX[[ 1]] - q[[1]])/ s[[1]] ] != 0,
		lX = Drop[lX, 1]
	];
	If[ FractionalPart[ (lX[[-1]] - q[[1]])/ s[[1]] ] != 0,
		lX = Drop[lX,-1]
	];

	lPX = Table[
			{lX[[k]], c[[2]] + Sqrt[r^2 - (lX[[k]] - c[[1]])^2]},
			{k, 1, Length[lX]}
		];

	lMX = Table[
			{lX[[k]], c[[2]] - Sqrt[r^2 - (lX[[k]] - c[[1]])^2]},
			{k, 1, Length[lX]}
		];

(*	points on the circle (c, r) that meet vertical mesh lines *)

	lX = Union[lPX, lMX];

(*	partition vertical projection of circle *)

	y1 = c[[2]] - r;
	y2 = c[[2]] + r;

	lY = partitionInterval[{y1, y2}, q[[2]], s[[2]]];

(*	drop endpoints if they are not on the mesh *)

	If[ FractionalPart[ (lY[[ 1]] - q[[2]])/ s[[2]] ] != 0,
		lY = Drop[lY, 1]
	];
	If[ FractionalPart[ (lY[[-1]] - q[[2]])/ s[[2]] ] != 0,
		lY = Drop[lY,-1]
	];

	lPY = Table[
			{c[[1]] + Sqrt[r^2 - (lY[[k]] - c[[2]])^2], lY[[k]]},
			{k, 1, Length[lY]}
		];

	lMY = Table[
			{c[[1]] - Sqrt[r^2 - (lY[[k]] - c[[2]])^2], lY[[k]]},
			{k, 1, Length[lY]}
		];

(*	points on the circle (c, r) that meet horizontal mesh lines *)

	lY = Union[lPY, lMY];

(*	take union, and add endpoints of the arc *)

	l = Union[lX, lY];

(*	Select points that lie on the specified arc; split arc in two when
	crossing over the branch cut of Arg; sort from th[[1]] to th[[2]]
	following specified orientation.
*)
	If[ or == CounterClockwise,
		If[ th[[1]] < th[[2]],
			l = Select[l,((th[[1]] <= ar[#, c]) && (ar[#, c] <= th[[2]]))&];
			l = Sort[l, (ar[#1, c] < ar[#2, c])&],

			l1 = Select[l,((th[[1]] <= ar[#, c]) && (ar[#, c] <= Pi))&];
			l1 = Sort[l1, (ar[#1, c] < ar[#2, c])&];

			l2 = Select[l,((-Pi < ar[#, c]) && (ar[#, c] <= th[[2]]))&];
			l2 = Sort[l2, (ar[#1, c] < ar[#2, c])&];

			l = Join[l1, l2]
		],

		If[ th[[1]] < th[[2]],
			l1 = Select[l,((-Pi < ar[#, c]) && (ar[#, c] <= th[[1]]))&];
			l1 = Sort[l1, (ar[#2, c] < ar[#1, c])&];

			l2 = Select[l,((th[[2]] <= ar[#, c]) && (ar[#, c] <= Pi))&];
			l2 = Sort[l2, (ar[#2, c] < ar[#1, c])&];

			l = Join[l1, l2],

			l = Select[l,((th[[2]] <= ar[#, c]) && (ar[#, c] <= th[[1]]))&];
			l = Sort[l, (ar[#2, c] < ar[#1, c])&]
		]
	];

(*	Add endpoints to l, if they're not already there; return l *)

	z = Map[Plus[c, #] &, r Map[{Cos[#], Sin[#]} &, th]];

	If[!MemberQ[l, z[[1]]], l = Prepend[l, z[[1]]]; ];
	If[!MemberQ[l, z[[2]]], l =  Append[l, z[[2]]]; ];

	l
]

(*	Func:	partionInterval

	Pre:	x	List {x1, x2} of two real numbers. (Here, "real" includes
				the "Integer", "Rational", and "Real" number types built in to
				Mathematica).

			q	Offset of a uniform 1D mesh.

			s	Stepsize of the mesh.

	Post:	Returns list containing the points x1, x2, and any elements of the
			mesh in between. The returned list is oriented, i.e., it's in the
			same order as {x1, x2}. If x1 = x2, x is returned.
*)

partitionInterval[x_List, q_, s_]:=
Module[{l = {}, m1, m2, y = Sort[x]},

	m1 = Ceiling[(y[[1]] - q)/s];
	m2 = Floor[(y[[2]] - q)/s];

	l = Select[
			Table[q + k s, {k, m1, m2}],
			( (y[[1]] < #) && (# < y[[2]]) )&
		];

	Which[
		x[[1]] == x[[2]], x,
		x[[1]] <  x[[2]], l = Sort[Union[x, l]],
		x[[2]] <  x[[1]], l = Sort[Union[x, l], (#2 < #1)&]
	]
]

(*	Func:	partitionLine

	Pre:	l	is a list of two points in the plane defining a line segment.

			q	is a list containing the horizontal and vertical offset of a
		   		mesh.

			s	is a list containing the horizontal and vertical stepsize of
				a mesh.

	Post:	Returns a list of points on the line line segment defined by "x"
			that intersect the mesh grid.
*)

partitionLine[l_, q_List, s_List]:=
Module[{x = l[[1]], lX = {}, lY = {}, m},

	Which[
		(x[[1, 1]] == x[[2, 1]]) && (x[[1, 2]] ==  x[[2, 2]]),
		x,

		(x[[1, 1]] != x[[2, 1]]) && (x[[1, 2]] ==  x[[2, 2]]),
		lX = partitionInterval[
				{x[[1, 1]], x[[2, 1]]}, q[[1]], s[[1]]
			];
		Table[{lX[[k]], x[[1, 2]]}, {k, 1, Length[lX]}],

		(x[[1, 1]] == x[[2, 1]]) && (x[[1, 2]] !=  x[[2, 2]]),
		lY = partitionInterval[
				{x[[1, 2]], x[[2, 2]]}, q[[2]], s[[2]]
			];
		Table[{x[[1, 1]], lY[[k]]}, {k, 1, Length[lY]}],

		(x[[1, 1]] != x[[2, 1]]) && (x[[1, 2]] !=  x[[2, 2]]),
		m = ( x[[2, 2]] - x[[1, 2]] ) / ( x[[2, 1]] - x[[1, 1]] );

		lX = partitionInterval[
				{x[[1, 1]], x[[2, 1]]}, q[[1]], s[[1]]
			];
		lX = Table[
			{lX[[k]], x[[1, 2]] + m ( lX[[k]] - x[[1, 1]] )},
			{k, 1, Length[lX]}
		];

		lY = partitionInterval[
				{x[[1, 2]], x[[2, 2]]}, q[[2]], s[[2]]
			];
		lY = Table[
			{x[[1, 1]] + (1/m) ( lY[[k]] - x[[1, 2]]), lY[[k]]},
			{k, 1, Length[lY]}
		];

		If[ x[[1, 1]] < x[[2, 1]],
			Sort[Union[lX, lY], (#1[[1]] < #2[[1]])&],
			Sort[Union[lX, lY], (#2[[1]] < #1[[1]])&]
		]
	]

]

(* ------------------------------------------------------------------------ *)
(* Short Auxiliary Functions                                                *)
(* ------------------------------------------------------------------------ *)

(*	Func:	norm

	Pre:	Vector with real components.

	Post:	Returns Euclidean norm.
*)

norm[b_List]:= Sqrt[b.b]

(* ------------------------------------------------------------------------ *)
(* ------------------------------------------------------------------------ *)

End[]

Protect[CostaSurface, ParallelSurface, CreateVertexData, GlueComponents, 
	Triangulate, CreatePolyhedron, ExportObjectData
]

EndPackage[]
