(* 
   Crippled Development Version 0.16
   05/6/2001
   The OpenGL features are turned off to allow use of
   current notebooks with MMA3+
   
   This package must be distributed without charge.
   Copyright is by the author.
   Changes must be notified in this header.
   
   This package is compatibel with MMA 4.0 upwards.
   
   Matthias Weber
   Mathematical Sciences Research Institute
   1000 Centennial Drive
   Berkeley, CA 94720
   USA
   
   email: weber@msri.org

*)
 
BeginPackage[
    "Own`Mesh0`", 
    {
      "Utilities`FilterOptions`"}
      ];

StereographicProjection::usage = 
	"StereographicProjection[z] returns the image of z under stereographic \
		projection.";
		
Mesh2D::usage = 
    "Mesh2D is a 2-dimensional graphics format for triangle meshes.";

Boundary::usage = 
    "Returns a list of vertice indices of edges for certain standard \
\
meshes.";

Corners::usage=
	"Corners[mesh] Returns a list of vertice indices of the corners of the mesh";

SplitEdge::usage = 
    "SplitEdge[mesh,e,pos] splits the boundary edge at position e in the \
\
boundary list of the Mesh2D object mesh into two edges at position pos.";

FindVertex::usage=
    "FindVertex[mesh,v] gives a list of boundary edges and positions of \
points \
of the vertex v in the boundary of mesh.";


Mesh3D::usage = 
    "Mesh3D is a 3-dimensional graphics format for triangle meshes.";

Vertices::usage = "Vertices[list] is a list of coordinates.";

Faces::usage = "Faces[list] is a list of triangles.";

Normals::usage = "Normals[list] is a list of normals.";

Edges::usage="Edges[list] is a list of edges given by pairs of vertice indices";

EdgeList::usage="EdgeList[mesh] will return a list of edges of a mesh.";

RectangularGrid::usage = 
    "RectangularGrid[x,y] creates a Mesh2D object representing a rectangular \
\
grid with vertical and horizontal lines given by the lists x and y.";

NRange::usage="NRange";

Subdivide::usage = 
    "Subdivide is an option for Grid generating functions, specifying wether \
\
the mesh is to subdivided into triangles (True) or not (False).";

PolarGrid::usage = 
    "PolarGrid[r,t] creates a Mesh2D object representing a circular grid with \
\
radial and circular lines given by the lists r and t.";

FirstQuadrantGrid::usage =
    "FirstQuadrantGrid[x1,b,x,y] generates a mesh in the first quadrant with \
\
polar coordinates around x1>0. b is a free positive real parameter \
\
influencing the shape of the grid, and x and y are lists of real numbers. \
The \
y-list should be a range from 0 to Pi.";

FullFirstQuadrantGrid::usage="not yet";

FullUpperHalfPlaneGrid::usage="Mesh in the upper half plane avoiding x1,x2, splitting
at 0, Infinity. Free parameter c.";

WedgeGrid::usage="WedgeGrid";

PolarWedgeGrid::usage=
	"PolarWedgeGrid[x,y] creates a triangle mesh with polar coordinates \
	in a wedge with radii x and angles y.";

QuarterDiskGrid::usage="QuarterDiskGrid[x,y,r] creates a triangle grid in \
  the first quadrant with corners at the x- and y-coordinates within a circle \
  of radius r.";

RectangularPolarGrid::usage=
	"RectangularPolarGrid[x,y,b,a] creates a grid in a rectangle with \
corners at 0,1,i b, (1+i)b using polar coordinates centered at a and (1-a)i+i b. \
The x and y lists  parameterize the gridlines.";

RectangularPolarGrid1::usage=
	"RectangularPolarGrid1[x,y,b,a] creates a grid in a rectangle with \
corners at 0,1,i b, (1+i)b using polar coordinates centered at 0 and 1+i b. \
The x and y lists  parameterize the gridlines.";

RectangularSubGrid::usage="RectangularSubGrid[{m,n},xi,yi] creates an edgelist \
  for a rectangular mesh of size m by n with the edges being specified by their 
  x- and y-coordinates.";

MeshApply::usage=
	"MeshApply[f,mesh] evaluates the complex function f on a 2-dimensional mesh
	 and returns the image mesh.";
	   
RepeatEdges::usage="Given a mesh with vs vertices and an edgelist edge based on \
this mesh, RepeatEdges[edges, vs, n] creates an edgelist based on a mesh which  \
 replicates the original mesh n times.";
 
ColorPlot::usage="ColorPlot[mesh,f[x,y],{x,y}] adds color information to a mesh, \
  based on the function f.";


TestMesh::usage="TestMesh[mesh] returns basic information about a mesh.";
  
AddBoundaryColor::"AddBoundaryColor[mesh] adds color at the boundary vertices.";

Mesh2DToTexturedGLGraphics::usage="Mesh2DToTexturedGLGraphics converts a Mesh2D \
	object in a GLGraphics object with textured polygons.";
	
Mesh2DToGraphics::usage = 
    "Mesh2dToGraphics[mesh] converts a Mesh2D object in a graphics object.";
    
Mesh2DToGLGraphics::usage = 
    "Mesh2dToGraphics[mesh] converts a Mesh2D object in a GLGraphics object.";

MeshPlot3D::usage = 
    "MeshPlot3D[fn,m] evaluates the function fn on the mesh vertices of m and \
\
returns a Mesh3D object.";

SymbolicMeshPlot3D::usage="SymbolicMeshPlot3D [fn,m] evaluates the function fn \
on the mesh vertices of m and returns a Mesh3D object.";

MeshJoin::usage = 
    "MeshJoin[mesh1,mesh2,merge] joins two meshes,thereby mergin the points \
\
in merge into pairs.";

MeshJoin2::usage = 
    "MeshJoin2[mesh1,mesh2,merge] joins two meshes,thereby mergin the points \
\
in merge into pairs.";

MeshJoin3::usage="testing";

Mesh3DToGraphics3D::usage = 
    "Mesh3DToGraphics3D[m] converts a Mesh3D object into a Graphics3D \
\
object.";

Mesh3DToGLGraphics::usage = 
    "Mesh3DToGLGraphics[m] converts a Mesh3D object into a GLGraphics \
\
object.";

BoundaryToGLGraphics::usage="none";

BoundaryToGraphics3D::usage="none";

Edges::usage="Edges";

EdgesToGLGraphics::usage="EdgesToGLGraphics";

EdgesToGraphics3D::usage="EdgesToGLGraphics";


MeshToVrml1::usage = "MeshToVrml1[file, mesh] exports a Mesh to Vrml1.";

MeshWrite::usage = 
    "MeshWrite[file,name,mesh] writes a Mesh3D object to a file with name \
\
name.";

MeshWrite2::usage = 
    "MeshWrite2[file,name,mesh] writes a Mesh3D object to a file with name \
\
name.";

MeshRotate::usage = "MeshRotate[mesh,edge] rotates a mesh about an edge";

MeshReflect::usage = "MeshReflect[mesh,plane] reflects a mesh at a plane";

MeshScrew::usage="MeshScrew[mesh, angle, trans] rotates mesh about the z-axes by angle \
and translates it vertically by trans.";
 
DetectSymmetryLine::usage="DetectSymmetryLine[mesh] finds symmetry boundary \
lines of meshes and returns  \
 the  corresponding normal.";
 
VerifyPlanarBoundaryCurve::usage="VerifyPlanarBoundaryCurve[mesh3D,pos,n] 
  returns the average distance of the boundary segment pos to the plane \
  \
given by normal n, together with the mean quadratic error.";
  
VerifyStraightBoundaryCurve::usage="VerifyStraightBoundaryCurve[mesh3D,pos,n] \

  returns the segment vertices of the boundary segment pos,  \
  together \
with the mean quadratic error with respect to the direction vector n.";
  
DetectBoundarySymmetry::usage="DetectBoundarySymmetry[mesh] returns a list of \
special \
  boundary components.";
 
DetectStraightBoundaryLine::usage="DetectStraightBoundaryLine[mesh2d] finds the \
	 straight boundary segments of a planar domain.";
	 
FindBoundarySymmetry::usage="FindBoundarySymmetry";

 
PlanarBoundaryCurve::usage="Primitive representing a plane.";

StraightBoundaryCurve::usage="Primitive representing a line.";

AdjustVertices::usage=
	"AdjustVertices[mesh,pos,v,n] sets the values of the vertices and normals \
	at pos the v and n.";
	
AdjustBoundaryCurve::usage="AdjustBoundaryCurve[mesh,pos,bdry] forces boundary components pos \
	to lie on a symmetry line or plane.";

GenerateDihedralSymmetry::usage="GenerateDihedralSymmetry[m,b1,b2] returns a \
  Mesh3D object by reflecting m repeatedly about the two symmetry planes b1, b2.";

ReflectAtBoundaryCurve::usage=
	"ReflectAtBoundaryCurve[m,b] reflects mesh m at the boundary curve b \
	and returns the joined mesh.";

FlipOrientation::usage = "FlipOrientation is an option for MeshReflect.";

GLShow::usage = "Fake GLShow function.";



Begin["`Private`"];

Unprotect[
	StereographicProjection,
	Vertices,
	Edges,
	Faces,
	Boundary,
	Corners,
	FindVertex,
	SplitEdge,
	SplitEdgeList,
    Mesh2D,
    Mesh3D,
    MeshApply,
    RectangularGrid,
    PolarGrid,
    PolarWedgeGrid,
    FirstQuadrantGrid,
    WedgeGrid,
    FullFirstQuadrantGrid,
    FullUpperHalfPlaneGrid,
    QuarterDiskGrid,
    RectangularPolarGrid,
    RectangularPolarGrid1,
    Mesh2DToGraphics,
    Mesh3DToGraphics3D,
    MeshPlot3D,
    RectangularSubGrid,
	RepeatEdges,
	ColorPlot,
	MeshPlot3D,
	NRange,
	TestMesh,
	AddBoundaryColor,
	SetBoundaryColor,
	MeshJoin2,
	MeshJoin3,
	MeshJoin,
	MeshRotate,
	MeshScrew,
	MeshReflect,
	GenerateDihedralSymmetry,
	ReflectAtBoundaryCurve,
	DetectSymmetryLine,
	VerifyPlanarBoundaryCurve,
	VerifyStraightBoundaryCurve,
	DetectBoundarySymmetry,
	DetectStraightBoundaryLine,
	FindBoundarySymmetry,
	AdjustVertices,
	AdjustStraightBoundaryCurve,
	AdjustBoundaryCurve,
	AdjustPlanarBoundaryCurve,
	Mesh2DToGraphics,
	Mesh2DToGLGraphics,
	Mesh2DToTexturedGLGraphics,
	Mesh3DToGraphics3D,
	Mesh3DToGLGraphics,
	BoundaryToGLGraphics,
	BoundaryToGraphics3D,
	EdgesToGLGraphics,
	EdgesToGraphics3D,
	EdgeList,
	MeshWrite,
	MeshWrite2,
	MeshToVrml1
];

eTolerance = 10.^-7;

norm[v_]:=v.v;

StereographicProjection[z_] := {(2*Re[z])/(1 + Abs[z]^2), 
   (2*Im[z])/(1 + Abs[z]^2), 1 - 2/(1 + Abs[z]^2)};

Vertices[Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]]] := v;

Vertices[Mesh2D[Vertices[v_], Colors[c_], Faces[f_], Boundary[b_]]] := v;

Faces[Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]]] := f;

Faces[Mesh2D[Vertices[v_], Colors[c_], Faces[f_], Boundary[b_]]] := f;

Boundary[Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]]] := b;

Boundary[Mesh2D[Vertices[v_], Colors[c_], Faces[f_],  Boundary[b_]]] := b;

Corners[mesh_] := First/@Boundary[mesh];

Vertices[Mesh3D[Vertices[v_], Faces[f_], Boundary[b_]]] := v;

Vertices[Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]]] := v;

Vertices[Mesh3D[Vertices[v_], Normals[n_], Colors[c_], Faces[f_], Boundary[b_]]] := v;

Boundary[Mesh3D[Vertices[v_],  Faces[f_], Boundary[b_]]] := b;

Boundary[Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]]] := b;

Boundary[Mesh3D[Vertices[v_], Normals[n_], Colors[c_], Faces[f_], Boundary[b_]]] := b;

(* Mesh2D utilities *)

FindVertex[Mesh2D[Vertices[v_], Colors[c_], Faces[f_],Boundary[b_]],p_,eps_:0.0001]:=
  Position[Map[v[[#]]&,b],_?(norm[#-p]<eps&),{2}];
  
FindVertex[Mesh2D[Vertices[v_],Faces[f_],Boundary[b_]],p_,eps_:0.0001]:=
  Position[Map[v[[#]]&,b],_?(norm[#-p]<eps&),{2}];
  
SplitEdge[Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]], {edge_, pos_}] := 
    Mesh2D[Vertices[v], Faces[f], 
      Boundary[Join[
          b[[Range[1, edge - 1]]], {b[[edge]][[Range[1, pos]]], 
            b[[edge]][[Range[pos, Length[b[[edge]]]]]]}, 
          b[[Range[edge + 1, Length[b]]]]]]];
          
SplitEdge[Mesh2D[Vertices[v_], Colors[c_], Faces[f_], Boundary[b_]], {edge_, pos_}] := 
    Mesh2D[Vertices[v], Colors[c],Faces[f], 
      Boundary[Join[
          b[[Range[1, edge - 1]]], {b[[edge]][[Range[1, pos]]], 
            b[[edge]][[Range[pos, Length[b[[edge]]]]]]}, 
          b[[Range[edge + 1, Length[b]]]]]]];
          
SplitEdgeList[Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]], {edge_, posl_}] := 
    Mesh2D[Vertices[v], Faces[f], 
      Boundary[Join[
          b[[Range[1, edge - 1]]], 
          b[[edge]][[#]]& /@ Apply[Range,Partition[Join[{1},posl,{Length[b[[edge]]]}],2,1],{1}], 
          b[[Range[edge + 1, Length[b]]]]]]];
          
SplitEdgeList[Mesh2D[Vertices[v_], Colors[c_], Faces[f_], Boundary[b_]], {edge_, posl_}] := 
    Mesh2D[Vertices[v], Colors[c], Faces[f], 
      Boundary[Join[
          b[[Range[1, edge - 1]]], 
          b[[edge]][[#]]& /@ Apply[Range,Partition[Join[{1},posl,{Length[b[[edge]]]}],2,1],{1}], 
          b[[Range[edge + 1, Length[b]]]]]]];
          
stitch[l1_,l2_]:=Module[{i},
    If[Length[l1]ēLength[l2],
      Join[
        Flatten[Table[{{l1[[i]],l1[[i+1]],l2[[i]]},{l1[[i+1]],l2[[i+1]],
                l2[[i]]}},{i,1,Length[l1]-1}],1],
        Table[{Last[l1],l2[[i+1]],l2[[i]]},{i,Length[l1],Length[l2]-1}]],
      Join[
        Flatten[Table[{{l1[[i]],l1[[i+1]],l2[[i]]},{l1[[i+1]],l2[[i+1]],
                l2[[i]]}},{i,1,Length[l2]-1}],1],
        Table[{l1[[i]],l1[[i+1]],Last[l2]},{i,Length[l2],Length[l1]-1}]]]];

(* mesh generating functions *)

RectangularTriangleGrid[x_, y_] :=
    Mesh2D[
      Vertices[Flatten[Transpose[Outer[List, x, y]], 1]], 
      Faces[Flatten[
          Map[{{#, # + 1, # + Length[x]}, {# + 1, # + Length[x] + 1, # + 
                    Length[x]}} &, 
            Flatten[Table[
                i Length[x] + j, {i, 0, Length[y] - 2}, {j, 1, 
                  Length[x] - 1}], 1]], 1]],
      Boundary[
        {
          Range[Length[x]],
          Length[x] Range[Length[y]],
          Range[Length[y]Length[x], (Length[y] - 1)Length[x] + 1, -1],
          Length[x]Range[Length[y] - 1, 0, -1] + 1
          }
        ]
      ];

RectangularQuadGrid[x_, y_] :=
    Mesh2D[
      Vertices[Flatten[Transpose[Outer[List, x, y]], 1]], 
      Faces[Flatten[
          Map[{{#, # + 1, # + Length[x] + 1, # + Length[x]}} &, 
            Flatten[Table[
                i Length[x] + j, {i, 0, Length[y] - 2}, {j, 1, 
                  Length[x] - 1}], 1]], 1]],
      Boundary[
        {
          Range[Length[x]],
          Length[x] Range[Length[y]],
          Range[Length[y]Length[x], (Length[y] - 1)Length[x] + 1, -1],
          Length[x]Range[Length[y] - 1, 0, -1] + 1
          }
        ]
      ];


RectangularGrid[x_, y_, opts___] := 
    With[{subdiv = Subdivide /. {opts} /. Options[RectangularGrid]}, 
      If[subdiv, RectangularTriangleGrid[x, y], RectangularQuadGrid[x, y]]];

Options[RectangularGrid] = {Subdivide -> False};

PolarTriangleGrid[x_, y_] := 
    Mesh2D[Vertices[
        Flatten[Transpose[Outer[{#1Cos[#2], #1Sin[#2]} &, x, y]], 1]], 
      Faces[Flatten[
          Map[{{#, # + 1, # + Length[x]}, {# + 1, # + Length[x] + 1, # + 
                    Length[x]}} &, 
            Flatten[Table[
                i Length[x] + j, {i, 0, Length[y] - 2}, {j, 1, 
                  Length[x] - 1}], 1]], 1]],
      Boundary[
        {
          Range[Length[x]],
          Length[x] Range[Length[y]],
          Range[Length[y]Length[x], (Length[y] - 1)Length[x] + 1, -1],
          Length[x]Range[Length[y] - 1, 0, -1] + 1
          }
        ]
      ];
      
PolarWedgeGrid[x_,y_]:=
    Mesh2D[
      Vertices[
        Prepend[
          Flatten[
            Transpose[
              Outer[
                {#1Cos[#2],#1Sin[#2]}&,
                Drop[x,1],y]
              ],1
            ],{0,0}
          ]
        ],
      Faces[
        Join[
          Table[
            {1,2+i(Length[x]-1),2+(i+1)(Length[x]-1)},
            {i,0,Length[y]-2}
            ],
          Flatten[
            Map[
              {{#,#+1,#+Length[x]-1},{#+1,#+Length[x],#+Length[x]-1}}&,
              Flatten[
                Table[
                  i (Length[x]-1)+j,{i,0,Length[y]-2},{j,2,Length[x]-1}],1
                ]
              ],1
            ]
          ]
        ],
      Boundary[
        {Range[Length[x]],
          (Length[x]-1) Range[Length[y]]+1,Append[
            Range[Length[y](Length[x]-1)+1,(Length[y]-1)(Length[x]-1)+2,-1],1
            ]
          }
        ]
      ];
      
 

PolarQuadGrid[x_, y_] := 
    Mesh2D[Vertices[
        Flatten[Transpose[Outer[{#1Cos[#2], #1Sin[#2]} &, x, y]], 1]], 
      Faces[Flatten[
          Map[{{#, # + 1, # + Length[x] + 1, # + Length[x]}} &, 
            Flatten[Table[
                i Length[x] + j, {i, 0, Length[y] - 2}, {j, 1, 
                  Length[x] - 1}], 1]], 1]],
      Boundary[{Range[Length[x]],
          Length[x] Range[Length[y]],
          Range[Length[y]Length[x], (Length[y] - 1)Length[x] + 1, -1],
          Length[x]Range[Length[y] - 1, 0, -1] + 1
          }
        ]
      ];

PolarGrid[x_, y_, opts___] := 
    With[{subdiv = Subdivide /. {opts} /. Options[PolarGrid]}, 
      If[subdiv, PolarTriangleGrid[x, y], PolarQuadGrid[x, y]]];


Options[PolarGrid] = {Subdivide -> False};

reim[z_] := {Re[z], Im[z]};

ComplexForm[fn_][{x_, y_}] := reim[fn[x + I y]];

MeshApply[fn_, Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]]] := 
    Mesh2D[Vertices[Map[ComplexForm[fn], v]], Faces[f], Boundary[b]];
    
MeshApply[fn_, Mesh2D[Vertices[v_], Colors[c_], Faces[f_], Boundary[b_]]] := 
    Mesh2D[Vertices[Map[ComplexForm[fn], v]], Colors[c], Faces[f], Boundary[b]];

roughinsert[l_, v_] := 
    Join[Select[l, # < v - eTolerance &], {v}, 
      Select[l, # > v + eTolerance &]];
      

roughposition[l_, v_] := Length[Select[l, # < v - eTolerance &]] + 1;


FirstQuadrantGrid[x1_Real, b_Real, xRng_, yRng_, opts___] := 
    Module[{xnewRng, pos},
      xnewRng = roughinsert[xRng, Log[1/b]];
      pos = Length[xnewRng] + 1 - roughposition[xRng, Log[1/b]];
      SplitEdge[
        MeshApply[
          Chop[x1 Sqrt[1 + b E^#]] &,
          RectangularGrid[
            xnewRng,
            yRng,
            opts]
          ],
        {3,pos}
        ]
      ];

Options[FirstQuadrantGrid] = {Subdivide -> False};

WedgeGrid[x1_Real, k_, xRng_, yRng_, opts___] := 
        MeshApply[
          Chop[x1(1 + 2/(-1 + I*E^#))^(2/k)] &,
          RectangularGrid[
            xRng,
            Pi/2.+Pi/2.*yRng,
            opts]
          ];
  

Options[WedgeGrid] = {Subdivide -> False};


FullFirstQuadrantGrid[x1_Real,  xRng_, yRng_,xsplt_:{}, opts___] := 
	Module[{x2,xnewRng, pos},
		x2=1/2*Log[(# + x1)/(# - x1)]&/@xsplt;
		xnewRng = Fold[roughinsert,xRng, x2];
		pos = Sort[Length[xnewRng] + Range[Length[xsplt]] - Map[roughposition[xRng,#]&,x2]];
        SplitEdgeList[MeshApply[
          			(x1 Tanh[#]) &,
          			RectangularGrid[
           				xnewRng,
            			yRng,
            			opts]],{3,pos}]
     ];
          


Options[FullFirstQuadrantGrid] = {Subdivide -> False};


FullUpperHalfPlaneGrid[x1_Real, x2_Real, c0_, xRng_, yRng_, opts___] := 
	Module[{v1,v2,xnewRng, pos},
		v1=Log[Abs[x1/c0]];
		v2=Log[Abs[x2/c0]];
		xnewRng = Fold[roughinsert,xRng, {v1,v2}];
		m1=RectangularGrid[
           				xnewRng,
            			yRng,
            			opts];
        pos=FindVertex[m1,{v1,Pi}][[1]];
        m1=SplitEdge[m1,pos];
        pos=FindVertex[m1,{v2,0}][[1]];
        m1=SplitEdge[m1,pos];
        MeshApply[(c0 E^#+x1)/(c0 E^#/x2+1) &,
          			m1]   
     ];
          


Options[FullUpperHalfPlaneGrid] = {Subdivide -> False};


QuarterDiskGrid[phi_,xl_,yl_,r_]:=Module[{xs,ys,vs,y,cnt,pat},
	xs=Union[r Cos/@phi,xl,{r}];
	ys=Union[r Sin/@phi,yl,{r}];
    vs=Map[
        Function[{y},{#,y}&/@
            Append[Select[xs,#^2+y^2<r^2+eTolerance&],Sqrt[r^2-y^2]]],
        ys];
    cnt=1;
    pat=Map[cnt++&,vs,{2}];
    
    Mesh2D[
      Vertices[Flatten[vs,1]],
      Faces[
        Flatten[Table[stitch[pat[[i]],pat[[i+1]]],{i,1,Length[pat]-1}],1]],
      Boundary[{pat[[1]],Last/@pat,Reverse[First/@pat]}]
      ]
    ];   
    
RectangularPolarGrid[xl_, yl_, b0_, a0_] := 
  Module[{r0, g2, g3, f1, h, h1, f2, h2, lambda, mo1, x0, mu, tr, ts, z, vs, 
      m1},
    r0 = 2/(1 - Sqrt[ModularLambda[1/b0 I]]) - 1;
    {g2, g3} = WeierstrassInvariants[{1, I b0}];
    f1[z_] := (EllipticF[ArcSin[z], 1/r0^2] + EllipticK[1/r0^2])/(2*
            EllipticK[1/r0^2]);
    h[z_] := -WeierstrassP[z, {g2, g3}] + WeierstrassP[1. + b0 I, {g2, g3}];
    lambda = -((2*h[1.])/(-1 + r0));
    mo1[z_] := (z + r0 lambda)/(-z + lambda);
    h1[z_] := mo1[h[z]];
    x0 = Re[h1[a0]];
    mu = 1/Im[h1[(1 + b0 I)/2.]];
    tr[z_] := (-x0 - r0 mu z)/(-1. + x0 mu  z);
    ts[z_] := (-x0 + z)/((r0 + x0*z)*mu);
    f2[z_] := f1[tr[z]];
    vs = Map[ts, {-r0, -1., 1., r0}];
    m1 = PolarGrid[Union[Exp[xl],Abs[vs]], yl];
    vv = FindVertex[m1, {vs[[1]], 0}][[1]];
    m1 = SplitEdge[m1, vv];
    vv = FindVertex[m1, {vs[[2]], 0}][[1]];
    m1 = SplitEdge[m1, vv];
    vv = FindVertex[m1, {vs[[3]], 0}][[1]];
    m1 = SplitEdge[m1, vv];
    vv = FindVertex[m1, {vs[[4]], 0}][[1]];
    m1 = SplitEdge[m1, vv];
    MeshApply[f2, m1]
    ]; 
    

RectangularPolarGrid1[xl_, yl_, b0_] := 
  Module[{a,f,x,x1,y1,g,xvs,m1,vv},
    a=-1/ModularLambda[I b0]+1;
    f[x_] := 2*I*EllipticF[ArcSin[Sqrt[x]/Sqrt[a]], a];
    x1=f[0]-f[a];
    y1=f[1]-f[0];
    g[x_]:=(-f[x]+y1-x1)/y1;
    xvs=Union[xl,{0,1.,Abs[a]}];
    
    m1 = PolarWedgeGrid[xvs,yl];
    m1=MeshApply[g, m1];
    vv = FindVertex[m1, {1, 0}][[1]];
    m1 = SplitEdge[m1, vv];
    vv = FindVertex[m1, {0,b0}][[1]];
    SplitEdge[m1, vv]
    ]; 

RectangularSubGrid[{m_,n_},xi_,yi_]:=
  Edges[Join[Map[m Range[0,n-1]+#&,xi],Map[Range[1,m]+(#-1)m&,yi]]];  

RepeatEdges[Edges[ed_],nv_,m_]:=Edges[Flatten[Map[ed+nv#&,Range[0,m-1]],1]];

(* Plotting functions *)

ColorPlot[{cn_}, {x_, y_}, m2_] := 
    Module[{c1,v,f,b},
      {v,f,b}={Vertices[m2],Faces[m2],Boundary[m2]};
      c1 = Function[{x, y}, cn];
      Mesh2D[
        Vertices[v],
        Colors[Apply[c1, v, {1}]],
        Faces[f],
        Boundary[b]
        ]
      ];

(* no normals, no colors: *)

MeshPlot3D[{fn_}, {x_, y_}, Mesh2D[Vertices[v_],Faces[f_],Boundary[b_]]] := 
    Module[{f1},
      f1 = Function[{x, y}, fn];
      Mesh3D[
        Vertices[Apply[f1, v, {1}]],
        Faces[f],
        Boundary[b]
        ]
      ];

(* no normals, yes colors: *) 
     
MeshPlot3D[{fn_}, {x_, y_}, Mesh2D[Vertices[v_],Colors[c_],Faces[f_],Boundary[b_]]] := 
    Module[{f1},
      f1 = Function[{x, y}, fn];
      Mesh3D[
        Vertices[Apply[f1, v, {1}]],
        Colors[c],
        Faces[f],
        Boundary[b]
        ]
      ];

(* yes normals, no colors: *)

MeshPlot3D[{fn_, nn_}, {x_, y_}, 
      Mesh2D[Vertices[v_],Faces[f_],Boundary[b_]]] := 
      Module[{f1, n1},
      f1 = Function[{x, y}, fn];
      n1 = Function[{x, y}, nn];
      Mesh3D[
        Vertices[Apply[f1, v, {1}]],
        Normals[Apply[n1, v, {1}]],
        Faces[f],
        Boundary[b]
        ]
      ];
      
(* yes normals, yes colors: *)
     
MeshPlot3D[{fn_, nn_}, {x_, y_},
   Mesh2D[Vertices[v_],Colors[c_],Faces[f_],Boundary[b_]]] := 
     Module[{f1, n1},
      f1 = Function[{x, y}, fn];
      n1 = Function[{x, y}, nn];
      Mesh3D[
        Vertices[Apply[f1, v, {1}]],
        Normals[Apply[n1, v, {1}]],
        Colors[c],
        Faces[f],
        Boundary[b]
        ]
      ];
      
      
SetAttributes[SymbolicMeshPlot3D, HoldFirst];
SetAttributes[MeshPlot3D, HoldFirst]; 
SetAttributes[ColorPlot, HoldFirst];

NRange[a_,b_,n_]:=Range[N[a],N[b],N[b-a]/N[n-1]];

TestMesh[m_] := 
  If[SameQ[Head[m],Mesh3D]||SameQ[Head[m],Mesh2D], Map[{Head[#], Length[#[[1]]]} &, List @@ m], False];

     
AddBoundaryColor[Mesh3D[Vertices[v_], Normals[n_], Faces[f_], 
		Boundary[b_]],cl_:{RGBColor[0,0,0],RGBColor[1,1,1]}]:=
	Module[{bd},
	bd=Union@@b;
	
	Mesh3D[
		Vertices[v], 
		Normals[n], 
		Colors[Map[If[MemberQ[bd,#],cl[[1]],cl[[2]]]&,Range[Length[v]]]],
		Faces[f], 
		Boundary[b]
	  ]
	];
	
SetBoundaryColor[Mesh3D[Vertices[v_], Normals[n_], Colors[c_], Faces[f_], 
		Boundary[b_]],cl_:RGBColor[0,0,0]]:=
	Module[{bd},
	bd=Union@@b;
	
	Mesh3D[
		Vertices[v], 
		Normals[n], 
		Colors[Map[If[MemberQ[bd,#],cl,c[[#]]]&,Range[Length[v]]]],
		Faces[f], 
		Boundary[b]
	  ]
	];

(* Mesh manipulation *)

MeshJoin2[Mesh3D[Vertices[v1_], Faces[f1_], Boundary[b1_]], 
      Mesh3D[Vertices[v2_], Faces[f2_], Boundary[b2_]], 
      merge_:{{},{}}] :=
    Mesh3D[
    	Vertices[Join[v1, v2]],
        Faces[
        	Join[f1, (Length[v1] + f2 /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], 
              	Flatten[b1[[#]]&/@merge[[1]]], List])]],
      	Boundary[
       	 Join[b1, (Length[v1] + b2 /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], 
              	Flatten[b1[[#]]&/@merge[[1]]], List])]]];
              
              
              
MeshJoin2[Mesh3D[Vertices[v1_], Normals[n1_], Faces[f1_], Boundary[b1_]], 
      Mesh3D[Vertices[v2_], Normals[n2_], Faces[f2_], Boundary[b2_]], 
      merge_:{{},{}}] :=
    Mesh3D[Vertices[Join[v1, v2]],
      Normals[Join[n1, n2]],
      Faces[
        Join[f1, ((Length[v1] + f2) /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List])]],
      Boundary[
        Join[b1, ((Length[v1] + b2) /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List])]]];

MeshJoin2[Mesh3D[Vertices[v1_], Normals[n1_], Colors[c1_], Faces[f1_], Boundary[b1_]], 
      Mesh3D[Vertices[v2_], Normals[n2_], Colors[c2_], Faces[f2_], Boundary[b2_]], 
      merge_:{{},{}}] :=
    Mesh3D[Vertices[Join[v1, v2]],
      Normals[Join[n1, n2]],
      Colors[Join[c1, c2]],
      Faces[
        Join[f1, (Length[v1] + f2 /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List])]],
      Boundary[
        Join[b1, (Length[v1] + b2 /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List])]]];

MeshJoin2[Mesh2D[Vertices[v1_], Faces[f1_], Boundary[b1_]], 
      Mesh2D[Vertices[v2_], Faces[f2_], Boundary[b2_]], 
      merge_:{{},{}}] :=
    Mesh2D[Vertices[Join[v1, v2]],
      Faces[
        Join[f1, (Length[v1] + f2 /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List])]],
      Boundary[
        Join[b1, (Length[v1] + b2 /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List])]]];
              
MeshJoin2[Mesh2D[Vertices[v1_],Colors[c1_], Faces[f1_], Boundary[b1_]], 
      Mesh2D[Vertices[v2_],Colors[c2_], Faces[f2_], Boundary[b2_]], 
      merge_:{{},{}}] :=
    Mesh2D[Vertices[Join[v1, v2]],
   		 Colors[Join[c1, c2]],
      Faces[
        Join[f1, (Length[v1] + f2 /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List])]],
      Boundary[
        Join[b1, (Length[v1] + b2 /. 
              Inner[Rule, Length[v1] + Flatten[b2[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List])]]];
                           
              
MeshJoin3[Mesh3D[Vertices[v1_],Normals[n1_],Faces[f1_],Boundary[b1_]],
      Mesh3D[Vertices[v2_],Normals[n2_],Faces[f2_],Boundary[b2_]],
      merge_:{{},{}}]:=
    Mesh3D[Vertices[Join[v1,v2]],Normals[Join[n1,n2]],
      Faces[Join[
          f1,((Length[v1]+f2)//.Inner[Rule,
                Length[v1]+Flatten[b2[[#]]&/@merge[[2]]],
                Flatten[b1[[#]]&/@merge[[1]]],List])]],
      Boundary[Join[
          b1,((Length[v1]+b2)//.Inner[Rule,
                Length[v1]+Flatten[b2[[#]]&/@merge[[2]]],
                Flatten[b1[[#]]&/@merge[[1]]],List])]]];
                
MeshJoin[m1_,m2_,merge_:{}]:=MeshJoin2[m1,m2,{merge,merge}];                      

MeshJoin[Mesh3D[Vertices[v1_], Faces[f1_], Boundary[b1_]],
      merge_:{{},{}}] :=
    Mesh3D[Vertices[v1],
      Faces[
        f1 /.Inner[Rule, Flatten[b1[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List]],
      Boundary[
        b1 /.Inner[Rule, Flatten[b1[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List]]];


                            
MeshJoin[Mesh3D[Vertices[v1_], Normals[n1_], Faces[f1_], Boundary[b1_]],
      merge_:{{},{}}] :=
    Mesh3D[Vertices[v1],
      Normals[n1],
      Faces[
        f1 /.Inner[Rule, Flatten[b1[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List]],
      Boundary[
        b1 /.Inner[Rule, Flatten[b1[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List]]];

MeshJoin[Mesh3D[Vertices[v1_], Normals[n1_], Colors[c1_], Faces[f1_], Boundary[b1_]],
      merge_:{{},{}}] :=
    Mesh3D[Vertices[v1],
      Normals[n1],
      Colors[c1],
      Faces[
        f1 /.Inner[Rule, Flatten[b1[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List]],
      Boundary[
        b1 /.Inner[Rule, Flatten[b1[[#]]&/@merge[[2]]], Flatten[b1[[#]]&/@merge[[1]]], List]]];



(* rotate p about line through v1 and v2 *)

rotate[p_, {v1_, v2_}] := -p + 2v1 + 
      2(p - v1).(v1 - v2)/((v1 - v2).(v1 - v2)) (v1 - v2);
 
(* differential of former map, use for normals *)
      
nrotate[p_, {v1_, v2_}] := -p + 
      2p.(v1 - v2)/((v1 - v2).(v1 - v2)) (v1 - v2);


MeshRotate[Mesh3D[Vertices[v_], Faces[f_], Boundary[b_]], StraightBoundaryCurve[{v1_, v2_}], 
      opts___] := 
    With[{flip = FlipOrientation /. {opts} /. Options[MeshRotate]}, Mesh3D[
        Vertices[Map[rotate[#, {v1, v2}] &, v]],
        Faces[If[flip, Map[Reverse, f], f]],
        Boundary[b]
        ]
      ];

MeshRotate[
      Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]], StraightBoundaryCurve[{v1_, v2_}], \

      opts___] := 
    With[{flip = FlipOrientation /. {opts} /. Options[MeshRotate]}, Mesh3D[
        Vertices[Map[rotate[#, {v1, v2}] &, v]],
        Normals[If[flip, -Map[nrotate[#, {v1, v2}] &, n], Map[nrotate[#, {v1, v2}] &, n]]],
        Faces[If[flip, Map[Reverse, f], f]],
        Boundary[b]
        ]
      ];
      
MeshRotate[
      Mesh3D[Vertices[v_], Normals[n_], Colors[c_], Faces[f_], Boundary[b_]], StraightBoundaryCurve[{v1_, v2_}], \

      opts___] := 
    With[{flip = FlipOrientation /. {opts} /. Options[MeshRotate]}, Mesh3D[
        Vertices[Map[rotate[#, {v1, v2}] &, v]],
        Normals[If[flip, -Map[nrotate[#, {v1, v2}] &, n], Map[nrotate[#, {v1, v2}] &, n]]],
        Colors[c],
        Faces[If[flip, Map[Reverse, f], f]],
        Boundary[b]
        ]
      ];

Options[MeshRotate] = {FlipOrientation -> True};

screw[{x_,y_,z_}, angle_, trans_]:={Cos[angle]x+Sin[angle]y,-Sin[angle]x+Cos[angle]y, z+trans};

MeshScrew[
      Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]], angle_, trans_] := 
    Mesh3D[
        Vertices[Map[screw[#,angle,trans] &, v]],
        Normals[Map[screw[#,angle,0.0] &, n]],
        Faces[f],
        Boundary[b]
        ];
      
      
reflect[p_,{pn_, d_}]:=p-2(p.pn-d)pn;

nreflect[p_,pn_]:=p-2(p.pn)pn;

MeshReflect[
	Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]], PlanarBoundaryCurve[pn_, d_], 
      opts___] := 
      With[{flip = FlipOrientation /. {opts} /. Options[MeshReflect]}, 
      Mesh2D[
        Vertices[Map[reflect[#, {pn, d}] &, v]],
        Faces[If[flip, Map[Reverse, f], f]],
        Boundary[b]
        ]
      ];
      
MeshReflect[
	Mesh2D[Vertices[v_], Colors[c_],Faces[f_], Boundary[b_]], PlanarBoundaryCurve[pn_, d_], 
      opts___] := 
      With[{flip = FlipOrientation /. {opts} /. Options[MeshReflect]}, 
      Mesh2D[
        Vertices[Map[reflect[#, {pn, d}] &, v]],
        Colors[c],
        Faces[If[flip, Map[Reverse, f], f]],
        Boundary[b]
        ]
      ];
           

MeshReflect[
      Mesh3D[Vertices[v_], Faces[f_], Boundary[b_]], PlanarBoundaryCurve[pn_, d_], 
      opts___] := 
    With[{flip = FlipOrientation /. {opts} /. Options[MeshReflect]}, Mesh3D[
        Vertices[Map[reflect[#, {pn, d}] &, v]],
        Faces[If[flip, Map[Reverse, f], f]],
        Boundary[b]
        ]
      ];
      
      
MeshReflect[
      Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]], PlanarBoundaryCurve[pn_, d_], 
      opts___] := 
    With[{flip = FlipOrientation /. {opts} /. Options[MeshReflect]}, Mesh3D[
        Vertices[Map[reflect[#, {pn, d}] &, v]],
        Normals[If[flip, Map[nreflect[#, pn] &, n], -Map[nreflect[#, pn] &, n]]],
        Faces[If[flip, Map[Reverse, f], f]],
        Boundary[b]
        ]
      ];
      
MeshReflect[
      Mesh3D[Vertices[v_], Normals[n_], Colors[c_], Faces[f_], Boundary[b_]], PlanarBoundaryCurve[pn_, d_], 
      opts___] := 
    With[{flip = FlipOrientation /. {opts} /. Options[MeshReflect]}, 
      Mesh3D[
        Vertices[Map[reflect[#, {pn, d}] &, v]],
        Normals[If[flip, Map[nreflect[#, pn] &, n], -Map[nreflect[#, pn] &, n]]],
        Colors[c],
        Faces[If[flip, Map[Reverse, f], f]],
        Boundary[b]
        ]
      ];
      
MeshReflect[
       PlanarBoundaryCurve[pn1_, d1_],PlanarBoundaryCurve[pn_, d_], 
      opts___] := 
    With[{flip = FlipOrientation /. {opts} /. Options[MeshReflect]}, 
      PlanarBoundaryCurve[
        nreflect[pn1, pn],
        d1-2d(pn1.pn)
        ]
      ];
 
Options[MeshReflect] = {FlipOrientation -> True};

equalPBC[PlanarBoundaryCurve[n1_,d1_],PlanarBoundaryCurve[n2_,d2_],eps_:0.001]:=
    norm[Cross[n1,n2]]+(d1-d2)^2<eps;

projection[a_,b_]:=Outer[Times,b,a];


(* detects if normal vectors lie near great circle *)

detectplanarnormal[pi_]:=Module[{ma,es,ld,n1},
    ma=Plus@@(projection[#,#]&/@pi);
    es=Eigensystem[ma];
    ld=Min[es[[1]]];
    n1=es[[2,Position[es[[1]],ld][[1,1]]]];
    {n1,ld}
    ];
    
(* detects if points lie near plane *)    
dtp[pi_]:=Module[{p0,ma,es,ld,n1},
      p0=Plus@@pi/Length[pi];
      ma=Plus@@(projection[#,#-p0]&/@pi);
      es=Eigensystem[ma];
      ld=Min[es[[1]]];
      n1=es[[2,Position[es[[1]],ld][[1,1]]]];
      {n1,p0.n1,ld}];

GenerateDihedralSymmetry[mp_,PlanarBoundaryCurve[n1_,d1_],
    PlanarBoundaryCurve[n2_,d2_],eps_:0.001]:=
  Module[{order,sb,pc,c1,c2,total,refl,ref1,ref2,len,i},
    order=Round[2Pi/ArcCos[Abs[n1.n2]]];
    sb=DetectBoundarySymmetry[mp,eps];
    pc=Union[Flatten[Position[sb,PlanarBoundaryCurve[_,_]],2]];
    c1=Select[pc,equalPBC[sb[[#,2]],PlanarBoundaryCurve[n1,d1],eps]&];
    ref1=sb[[c1[[1]],2]];
    c2=Select[pc,equalPBC[sb[[#,2]],PlanarBoundaryCurve[n2,d2],eps]&];
    ref2=sb[[c2[[1]],2]];
    len=Length[Boundary[mp]];
    total=mp;
    refl=mp;
    ct=c1;
    For[i=1,i<order,i++,
      refl=MeshReflect[refl,ref1];
      total=MeshJoin2[total,refl,{c1+(i-1)len,c1}];
      {c1,c2}={c2,c1};
      {ref1,ref2}={MeshReflect[ref2,ref1],ref1};
      ];
    MeshJoin[total,{c1,c1+(order-1) len}]
    ]; 
    
ReflectAtBoundaryCurve[mp_,PlanarBoundaryCurve[n1_,d1_],eps_:0.001]:=
  Module[{sb,pc,c1,total,refl,ref1},
    sb=DetectBoundarySymmetry[mp,eps];
    pc=Union[Flatten[Position[sb,PlanarBoundaryCurve[_,_]],2]];
    c1=Select[pc,equalPBC[sb[[#,2]],PlanarBoundaryCurve[n1,d1],eps]&];
    ref1=sb[[c1[[1]],2]];
    total=mp;
    ct=c1;
    refl=MeshReflect[mp,ref1];
    total=MeshJoin[total,refl,{c1}]];

DetectSymmetryLine[Mesh3D[Vertices[v_],Normals[n_],Faces[f_],Boundary[b_]]]:=
  Map[detectplanarnormal,Map[n[[#]]&,b]];
  
DetectSymmetryLine[Mesh3D[Vertices[v_],Normals[n_],Colors[c_],Faces[f_],Boundary[b_]]]:=
  Map[detectplanarnormal,Map[n[[#]]&,b]];



VerifyPlanarBoundaryCurve[
    Mesh3D[Vertices[v_],Normals[n_],Faces[f_],Boundary[b_]],pos_,n0_]:=
  Module[{d,e},
    d=Plus@@((n0.#)&/@v[[b[[pos]]]])/Length[b[[pos]]];
    e=Plus@@(((n0.#)-d)^2&/@v[[b[[pos]]]]);
    {d,e}];
    
VerifyPlanarBoundaryCurve[
    Mesh3D[Vertices[v_],Normals[n_],Colors[c_],Faces[f_],Boundary[b_]],pos_,n0_]:=
  Module[{d,e},
    d=Plus@@((n0.#)&/@v[[b[[pos]]]])/Length[b[[pos]]];
    e=Plus@@(((n0.#)-d)^2&/@v[[b[[pos]]]]);
    {d,e}];
  
VerifyStraightBoundaryCurve[
    Mesh3D[Vertices[v_],Normals[n_],Faces[f_],Boundary[b_]],pos_,n0_]:=
  Module[{d,e},
    s={v[[b[[pos,1]]]],v[[b[[pos,Length[b[[pos]]]]]]]};
    e=Sum[
        norm[Cross[n0,(v[[b[[pos,i]]]]-v[[b[[pos,i+1]]]])]],{i,1,
          Length[b[[pos]]]-1}];
    {s,e}];
    
VerifyStraightBoundaryCurve[
    Mesh3D[Vertices[v_],Normals[n_],Colors[c_],Faces[f_],Boundary[b_]],pos_,n0_]:=
  Module[{d,e},
    s={v[[b[[pos,1]]]],v[[b[[pos,Length[b[[pos]]]]]]]};
    e=Sum[
        norm[Cross[n0,(v[[b[[pos,i]]]]-v[[b[[pos,i+1]]]])]],{i,1,
          Length[b[[pos]]]-1}];
    {s,e}];
    
DetectBoundarySymmetry[m_,eps_:0.0001]:=
  Module[{ls,i,ver},
    ls=DetectSymmetryLine[m];
    Table[{i,If[ls[[i,2]]>eps,False,
        ver=
          VerifyPlanarBoundaryCurve[
            m,i,ls[[i,1]]];
        If[ver[[2]]<eps,
          PlanarBoundaryCurve[ls[[i,1]],ver[[1]]],
          ver=VerifyStraightBoundaryCurve[m,i,ls[[i,1]]];
          If[ver[[2]]<eps,
            StraightBoundaryCurve[ver[[1]]],{False}]]]},{i,1,Length[Boundary[m]\
]}]];

(* works for Mesh2D: *)

DetectStraightBoundaryLine[m_,eps_:0.0001]:=
  Module[{ls,i},
    ls=Map[dtp,Map[Vertices[m][[#]]&,Boundary[m]]];
    Table[{i,If[ls[[i,3]]>eps,False,
          PlanarBoundaryCurve[ls[[i,1]],ls[[i,2]]]]},{i,1,Length[Boundary[m]\
]}]];



FindBoundarySymmetry[m_,PlanarBoundaryCurve[n_,d_],eps_:0.0001]:=
  Module[{i,ver},
    Flatten[Position[Table[
        ver=
          VerifyPlanarBoundaryCurve[
            m,i,n];
        ver[[2]]<0.0001 &&Abs[ver[[1]]-d]<eps
          ,{i,1,Length[Boundary[m]]}],_?(#&)]
]];

FindBoundarySymmetry[m_,StraightBoundaryCurve[{p1_,p2_}],eps_:0.0001]:=
  Module[{i,ver,v},
    Flatten[Position[Table[
        ver=
          VerifyStraightBoundaryCurve[
            m,i,p1-p2];v=ver[[1,1]];
        ver[[2]]<eps &&
          norm[v-(((v-p1).(p2-p1)) p2-((v-p2).(p2-p1)) p1)/norm[p2-p1]]<eps
          ,{i,1,Length[Boundary[m]]}
],_?(#&)]]];
  
(* Mesh adjustment *)

projv[v_,StraightBoundaryCurve[{p1_, p2_}]]:=
	(((v-p1).(p2-p1)) p2-((v-p2).(p2-p1)) p1)/norm[p2-p1];

normalize[v_]:=v/Sqrt[norm[v]];
	
projn[n_, StraightBoundaryCurve[{p1_, p2_}]]:=normalize[n-(n.(p2-p1))/((p2-p1).(p2-p1))(p2-p1)];

projv[v_, PlanarBoundaryCurve[n1_,d1_]]:=
	v-((v.n1)-d1)n1;
	
projn[n_, PlanarBoundaryCurve[n1_,d1_]]:=
	normalize[n-(n.n1)/(n1.n1)n1];


AdjustVertices[Mesh3D[Vertices[v_],Normals[n_],Faces[f_],Boundary[b_]],vpos_,
      vval_,nval_]:=Module[{v1=v,n1=n},
      v1[[vpos]]=vval;
      n1[[vpos]]=nval;
      Mesh3D[Vertices[v1],Normals[n1],Faces[f],Boundary[b]]
      ];

AdjustStraightBoundaryCurve[
    	Mesh3D[Vertices[v_],Normals[n_],Faces[f_], Boundary[b_]],
    	pos_,
    	StraightBoundaryCurve[{v1_, v2_}]
	]:=Module[{ind},
    	ind=Union@@(b[[#]]&/@pos); (* adjust these vertices and normals *)
   		Mesh3D[
    		Vertices[Map[If[MemberQ[ind,#],projv[v[[#]],StraightBoundaryCurve[{v1, v2}]],v[[#]]]&,
				Range[Length[v]]]],
    		Normals[Map[If[MemberQ[ind,#],projn[n[[#]],StraightBoundaryCurve[{v1, v2}]],n[[#]]]&,
			Range[Length[n]]]],
    		Faces[f],
    		Boundary[b]
    	]
    ];
    
    
AdjustStraightBoundaryCurve[
    Mesh3D[Vertices[v_],Normals[n_],Colors[c_],Faces[f_], Boundary[b_]],
    pos_,
    StraightBoundaryCurve[{v1_, v2_}]
]:=Module[{ind},
   	ind=Union@@(b[[#]]&/@pos); (* adjust these vertices and normals *)
    Mesh3D[
    	Vertices[Map[If[MemberQ[ind,#],projv[v[[#]],StraightBoundaryCurve[{v1, v2}]],v[[#]]]&,
			Range[Length[v]]]],
    	Normals[Map[If[MemberQ[ind,#],projn[n[[#]],StraightBoundaryCurve[{v1, v2}]],n[[#]]]&,
			Range[Length[n]]]],
    	Colors[c],
    	Faces[f],
    	Boundary[b]
    	]
    ];
    
AdjustPlanarBoundaryCurve[
    Mesh3D[Vertices[v_],Normals[n_],Faces[f_], Boundary[b_]],
    pos_,
    PlanarBoundaryCurve[n1_,d1_]
]:=Module[{ind},
    ind=Union@@(b[[#]]&/@pos); (* adjust these vertices and normals *)
    Mesh3D[
    	Vertices[Map[If[MemberQ[ind,#],projv[v[[#]],PlanarBoundaryCurve[n1,d1]],v[[#]]]&,Range[Length[v]]]],
    	Normals[Map[If[MemberQ[ind,#],projn[n[[#]],PlanarBoundaryCurve[n1,d1]],n[[#]]]&,Range[Length[n]]]],
    	Faces[f],
    	Boundary[b]
    ]
    ];
    
    
AdjustPlanarBoundaryCurve[
    Mesh3D[Vertices[v_],Normals[n_],Colors[c_],Faces[f_], Boundary[b_]],
    pos_,
    PlanarBoundaryCurve[n_,d_]
]:=Module[{ind},
    ind=Union@@(b[[#]]&/@pos); (* adjust these vertices and normals *)
    Mesh3D[
    	Vertices[Map[If[MemberQ[ind,#],projv[v[[#]],PlanarBoundaryCurve[n1,d1]],v[[#]]]&,Range[Length[v]]]],
    	Normals[Map[If[MemberQ[ind,#],projn[n[[#]],PlanarBoundaryCurve[n1,d1]],n[[#]]]&,Range[Length[n]]]],
    	Colors[c],
    	Faces[f],
    	Boundary[b]
    ]
    ];
    
 
AdjustBoundaryCurve[m_, pos_,PlanarBoundaryCurve[n_,d_]]:=
	AdjustPlanarBoundaryCurve[m,pos,PlanarBoundaryCurve[n,d]];
	
AdjustBoundaryCurve[m_, pos_,StraightBoundaryCurve[{v1_, v2_}]]:=
	AdjustStraightBoundaryCurve[m,pos,StraightBoundaryCurve[{v1, v2}]];

(* Conversion to display graphics formats *)

insertValues[f_, v_] := Map[v[[#]] &, f];

conv2D3D[l_]:=Map[Append[#,0]&,l];




Mesh2D /: MakeBoxes[Mesh2D[___], StandardForm] := 
    RowBox[{"\[SkeletonIndicator]", "Mesh2D", "\[SkeletonIndicator]"}];

Mesh3D /: MakeBoxes[Mesh3D[___], StandardForm] := 
    RowBox[{"\[SkeletonIndicator]", "Mesh3D", "\[SkeletonIndicator]"}];

Edges /: MakeBoxes[Edges[___], StandardForm] := 
    RowBox[{"\[SkeletonIndicator]", "Edges", "\[SkeletonIndicator]"}];


Mesh2DToGraphics[Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]]] := 
    Graphics[Map[Line, insertValues[edgelist[f], v]]];
    
Mesh2DToGraphics[Mesh2D[Vertices[v_], Colors[c_], Faces[f_], Boundary[b_]],cfn_] := 
    Graphics[Map[{cfn[(Plus@@c[[#]])/Length[#]],Polygon[v[[#]]]}&,f]];
    
Mesh2DToGLGraphics=Mesh2DToGraphics;

Mesh2DToTexturedGLGraphics[Mesh2D[Vertices[v_], Faces[f_], Boundary[b_]],tex_,dir_:{}] := 
        	If[Length[f[[1]]] == 3, 
     		 GLGraphics[Join[dir,{tex},{
        		Triangles[insertValues[Flatten[f], conv2D3D[v]], TextureCoords->{1,1}]}]], 
      		GLGraphics[Join[dir,{tex},{
        		Quads[insertValues[Flatten[f], conv2D3D[v]], TextureCoords->{1,1}]}]]
        ];
        
        



Mesh3DToGraphics3D[Mesh3D[Vertices[v_], Faces[f_], Boundary[b_]],dir_:{}] := 
    Graphics3D[Append[dir,Map[Polygon, insertValues[f, v]]]];
    
Mesh3DToGraphics3D[Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]],dir_:{}] := 
    Graphics3D[Append[dir,Map[Polygon, insertValues[f, v]]]];

Mesh3DToGLGraphics=Mesh3DToGraphics3D;



rescale[w_, e_] := w^e/(1 + w^e);
        

BoundaryToGraphics3D[
      Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]], bi_, dir_:{}] := 
    Graphics3D[Join[dir,
        Map[Line[insertValues[b[[#]], v]]&,bi]]];
        
BoundaryToGLGraphics=BoundaryToGraphics3D;        

BoundaryToGraphics3D[
      Mesh3D[Vertices[v_], Faces[f_], Boundary[b_]], bi_, dir_:{}] := 
    Graphics3D[Join[dir,
        Map[Line[insertValues[b[[#]], v]]&,bi]]];
         


 
 EdgesToGraphics3D[Mesh3D[Vertices[v_],Normals[n_],Faces[f_],Boundary[b_]],
      Edges[ed_] ,dir_:{}]:=
    Graphics3D[
      Join[dir,Map[
          Line[insertValues[#,v]]&,ed]]];
          
 EdgesToGraphics3D[Mesh3D[Vertices[v_],Normals[n_],Colors[c_], Faces[f_],Boundary[b_]],
      Edges[ed_] ,dir_:{}]:=
    Graphics3D[
      Join[dir,Map[
          Line[insertValues[#,v]]&,ed]]];
          
  EdgesToGraphics3D[Mesh3D[Vertices[v_],Faces[f_],Boundary[b_]],
      Edges[ed_] ,dir_:{}]:=
    Graphics3D[
      Join[dir,Map[
          Line[insertValues[#,v]]&,ed]]];
                           
 EdgesToGLGraphics=EdgesToGraphics3D;
                            
(* Export functions *)

ClosedLine[l_] := Append[l, First[l]];

getEdges[l_] := Thread[{l, RotateLeft[l]}];

SortedEdges[l_] := Map[Sort, getEdges[l]];

edgelist[f_] := Union[Flatten[Map[SortedEdges, f], 1]];

EdgeList[Mesh3D[Vertices[v_], Normals[n_], Colors[c_], Faces[f_], 
        Boundary[b_]]]:=
        Edges[edgelist[f]];
        
EdgeList[Mesh3D[Vertices[v_], Normals[n_], Faces[f_], 
        Boundary[b_]]]:=
        Edges[edgelist[f]];
        
EdgeList[Mesh3D[Vertices[v_], Faces[f_], 
        Boundary[b_]]]:=
        Edges[edgelist[f]];
        
EdgeList[Mesh2D[Vertices[v_], Faces[f_], 
        Boundary[b_]]]:=
        Edges[edgelist[f]];


CompactForm[x_, digits_] := 
    SetPrecision[Chop[x], digits] // CForm // ToString;

PointToString[pts_List, prec_Integer, tchar_:" "] := 
    StringJoin[{" ", CompactForm[#, prec]} & /@ pts, tchar];


MeshWrite[file_, name_, 
      Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]]] :=
    
    Module[
      {strm = OpenWrite[file, FormatType -> FortranForm],
        i},
      
      WriteString[strm, "createTopSegment verts ", "\n"];
      WriteString[strm, "createTopSegment normals ", "\n"];
      WriteString[strm, "createTopSegment faces ", "\n"];
      WriteString[strm, "\n"];
      
      WriteString[strm, "setSegment verts\n"];
      WriteString[strm, "c 1 0 0\n"];
      For[i = 1, i <= Length[v], i++,
        WriteString[strm, "v ", i - 1, " ", PointToString[v[[i]], 4], "\n"]
        ];
      
      WriteString[strm, "\n"];
      WriteString[strm, "setSegment faces\n"];
      WriteString[strm, "setVertSegment verts\n"];
      WriteString[strm, "c 0 1 0\n"];
      For[i = 1, i <= Length[f], i++,
        WriteString[strm, "t ", i - 1, " ", f[[i, 1]] - 1, " ", f[[i, 2]] - \
1,
           " ", f[[i, 3]] - 1, "\n"]];
      
      Close[strm]];

MeshWrite2[file_, name_, 
      Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]]] :=
    
    Module[
      {strm = OpenWrite[file, FormatType -> FortranForm],
        i,
        edges = edgelist[f]},
      
      WriteString[strm, "createTopSegment verts ", "\n"];
      WriteString[strm, "createTopSegment edges ", "\n"];
      WriteString[strm, "createTopSegment normals ", "\n"];
      WriteString[strm, "createTopSegment faces ", "\n"];
      WriteString[strm, "\n"];
      
      WriteString[strm, "setSegment verts\n"];
      WriteString[strm, "c 1 0 0\n"];
      For[i = 1, i <= Length[v], i++,
        WriteString[strm, "v ", i - 1, " ", PointToString[v[[i]], 4], "\n"]
        ];
      
      WriteString[strm, "setSegment edges\n"];
      WriteString[strm, "setVertSegment verts\n"];
      WriteString[strm, "c 0 0 1\n"];
      For[i = 1, i <= Length[edges], i++,
        WriteString[strm, "e ", i - 1, " ", edges[[i, 1]] - 1, " ", 
          edges[[i, 2]] - 1, "\n"]
        ];
      
      WriteString[strm, "\n"];
      WriteString[strm, "setSegment faces\n"];
      WriteString[strm, "setVertSegment verts\n"];
      WriteString[strm, "c 0 1 0\n"];
      For[i = 1, i <= Length[f], i++,
        WriteString[strm, "t ", i - 1, " ", f[[i, 1]] - 1, " ", f[[i, 2]] - \
1,
           " ", f[[i, 3]] - 1, "\n"]];
      
      Close[strm]];


MeshToVrml1[file_, 
      Mesh3D[Vertices[v_], Normals[n_], Faces[f_], Boundary[b_]]] :=
    
    Module[
      {strm = OpenWrite[file, FormatType -> FortranForm],
        lf,
        i},
      
      WriteString[strm, "#VRML V1.0 ascii", "\n"];
      WriteString[strm, "\n"];
      WriteString[strm, "Separator {", "\n"];
      WriteString[strm, "Coordinate3 {point[", "\n"];
      WriteString[strm, "\n"];
      
      For[i = 1, i <= Length[v] - 1, i++,
        WriteString[strm, PointToString[v[[i]], 4], ",\n"]
        ];
      WriteString[strm, PointToString[Last[v], 4], "\n"];
      WriteString[strm, "]}\n"];
      
      WriteString[strm, "\n"];
      WriteString[strm, "Normal {vector[\n"];
      
      For[i = 1, i <= Length[n] - 1, i++,
        WriteString[strm, PointToString[n[[i]], 4], ",\n"]
        ];
      WriteString[strm, PointToString[Last[n], 4], "\n"];
      WriteString[strm, "]}\n"];
      
      WriteString[strm, "\n"];
      WriteString[strm, "IndexedFaceSet {\n"];
      WriteString[strm, "coordIndex[\n"];
      
      lf = Length[f];
      If[Length[f[[1]]] == 3,
        For[i = 1, i <= lf - 1, i++,
          
          WriteString[strm, f[[i, 1]] - 1, ", ", f[[i, 2]] - 1, ", ", 
            f[[i, 3]] - 1, ", -1,\n"]
          ];
        WriteString[strm, f[[lf, 1]] - 1, ", ", f[[lf, 2]] - 1, ", ", 
          f[[lf, 3]] - 1, ", -1\n"],(*else*)
        
        For[i = 1, i <= lf - 1, i++,
          
          WriteString[strm, f[[i, 1]] - 1, ", ", f[[i, 2]] - 1, ", ", 
            f[[i, 3]] - 1, ", ", f[[i, 4]] - 1, ", -1,\n"]];
        WriteString[strm, f[[lf, 1]] - 1, ", ", f[[lf, 2]] - 1, ", ", 
          f[[lf, 3]] - 1, ", ", f[[lf, 4]] - 1, ", -1\n"]
        ];
      
      WriteString[strm, "]}}\n"];
      
      Close[strm]
      ];

GLShow=Show;

Protect[
	StereographicProjection,
	Vertices,
	Edges,
	Faces,
	Boundary,
	Corners,
	FindVertex,
	SplitEdge,
	SplitEdgeList,
    Mesh2D,
    Mesh3D,
    MeshApply,
    RectangularGrid,
    PolarGrid,
    PolarWedgeGrid,
    FirstQuadrantGrid,
    WedgeGrid,
    FullFirstQuadrantGrid,
    FullUpperHalfPlaneGrid,
    QuarterDiskGrid,
    RectangularPolarGrid,
    RectangularPolarGrid1,
    Mesh2DToGraphics,
    Mesh3DToGraphics3D,
    MeshPlot3D,
    RectangularSubGrid,
	RepeatEdges,
	ColorPlot,
	MeshPlot3D,
	NRange,
	TestMesh,
	AddBoundaryColor,
	SetBoundaryColor,
	MeshJoin2,
	MeshJoin3,
	MeshJoin,
	MeshRotate,
	MeshScrew,
	MeshReflect,
	GenerateDihedralSymmetry,
	ReflectAtBoundaryCurve,
	DetectSymmetryLine,
	VerifyPlanarBoundaryCurve,
	VerifyStraightBoundaryCurve,
	DetectBoundarySymmetry,
	DetectStraightBoundaryLine,
	FindBoundarySymmetry,
	AdjustVertices,
	AdjustStraightBoundaryCurve,
	AdjustBoundaryCurve,
	AdjustPlanarBoundaryCurve,
	Mesh2DToGraphics,
	Mesh2DToGLGraphics,
	Mesh2DToTexturedGLGraphics,
	Mesh3DToGraphics3D,
	Mesh3DToGLGraphics,
	BoundaryToGLGraphics,
	BoundaryToGraphics3D,
	EdgesToGLGraphics,
	EdgesToGraphics3D,
	EdgeList,
	MeshWrite,
	MeshWrite2,
	MeshToVrml1
];


End[];

EndPackage[];