(* 2024-04-02
random collection of WolframLang code.
from

Xah Lee WolframLang Packages
http://xahlee.info/M/xah_lee_wolframlang_packages.html

*)

(* ::Package:: *)

BeginPackage["Transform2DPlot`"] ;



Begin["`Private`"];

Clear[graphicsBoundary];

graphicsBoundary[gra_]:=({Min@#,Max@#}&)/@Transpose@(#/.Point->List&)@Flatten@Cases[gra/.Rectangle[{xm_,ym_},{xM_,yM_}]:>Polygon[{{xm,ym},{xM,ym},{xM,yM},{xm,yM},{xm,ym}}],Polygon[pts_]|Line[pts_]:>(Point@@#&)/@pts,Infinity]

Clear[boundaryPoints];

boundaryPoints[li_List]:=Flatten@Outer[Point,Sequence@@li]/.Point->List

Clear[insertMiddlePoint];

insertMiddlePoint[a_,b_,maxDist_,f_,0]=Point@b;

insertMiddlePoint[a_,b_,maxDist_,f_,n_]:=
If[Sqrt[# . #]&[f[a]-f[b]]>maxDist,{insertMiddlePoint[a,(a+b)/2,maxDist,f,n-1],insertMiddlePoint[(a+b)/2,b,maxDist,f,n-1]},Point@b]

Clear[linearInterpolate];

linearInterpolate[pts_List,maxDist_]:=Module[{res={},p1,p2,n},
Do[{p1,p2}=pts[[{i,i+1}]];
If[(n=Ceiling[(Sqrt[#1 . #1]&)[p2-p1]/maxDist])>1,res={res,Table[Point[p1+(p2-p1) j],{j,0,1-1/(2 n),1/n}]},res={res,Point[p1]}],{i,1,Length[pts]-1}];(#1/. Point[t_]->t&)[Flatten[{res,Point[Last[pts]]}]]]

linearInterpolate[pts_List,maxDist_,f_,n_]:=(#1/. Point[t_]->t&)[Flatten[{Point[First[pts]],(insertMiddlePoint[Sequence@@#1,maxDist,f,n]&)/@Partition[pts,2,1]}]]

Clear[breakCircleDiskRectangle];
breakCircleDiskRectangle[gra_,maxDist_/;N[maxDist]!=0]:=With[{md=maxDist},gra/. {(Circle|Disk)[{c1_,c2_},0|0.,___]:>Point[{c1,c2}],Circle[{c1_,c2_},r_?NumberQ]:>With[{iStep=(2 Pi)/Ceiling[(2 r Pi)/md]},Line[Table[r {Cos[t],Sin[t]}+{c1,c2},{t,0,2 Pi,iStep}]]],Circle[{c1_,c2_},{a_,b_}]:>With[{iStep=(2 Pi)/Ceiling[(2 Pi Max[a,b])/md]},Line[Table[{a Cos[t],b Sin[t]}+{c1,c2},{t,0,2 Pi,iStep}]]],Circle[{c1_,c2_},r_,{a_,b_}]:>With[{iStep=(b-a)/Ceiling[(2 r Pi)/md]},Line[Table[r {Cos[t],Sin[t]}+{c1,c2},{t,a,b,iStep}]]],Disk[{c1_,c2_},r_?NumberQ]:>With[{iStep=(2 Pi)/Ceiling[(2 r Pi)/md]},Polygon[Table[r {Cos[t],Sin[t]}+{c1,c2},{t,0,2 Pi,iStep}]]],Disk[{c1_,c2_},{a_,b_}]:>With[{iStep=(2 Pi)/Ceiling[(2 Pi Max[a,b])/md]},Polygon[Table[{a Cos[t],b Sin[t]}+{c1,c2},{t,0,2 Pi,iStep}]]],Disk[{c1_,c2_},r_,{a_,b_}]:>With[{iStep=(b-a)/Ceiling[(2 r Pi)/md]},Polygon[Table[r {Cos[t],Sin[t]}+{c1,c2},{t,a,b,iStep}]]],Rectangle[{xm_,ym_},{xM_,yM_}]:>Polygon[{{xm,ym},{xM,ym},{xM,yM},{xm,yM}}]}]

(* N[gra_,Infinity]:=gra; *)

IncreaseGraphicsResolution[gra_,maxDist_/;(N@maxDist!=0)]:=With[{md=maxDist},(breakCircleDiskRectangle[gra,maxDist]/.{Polygon[pts_]:>Polygon@(Drop[#,-1]&)@(linearInterpolate[#,md]&)@Append[pts,First@pts],Line[pts_]:>Line[linearInterpolate[pts,md]]})]

IncreaseGraphicsResolution[gra_,maxDist_/;(N@maxDist!=0),f_,n_]:=breakCircleDiskRectangle[gra,maxDist]/.{Polygon[points_]:>Polygon@linearInterpolate[Append[points,First@points],maxDist,f,n],Line[points_]:>Line[linearInterpolate[points,maxDist,f,n]]}

PolarGrid[{rMin_,rMax_,rStep_:Automatic},{thetaMin_,thetaMax_,thetaStep_:Automatic}]:=PolarGrid[{rMin,rMax,rStep},{thetaMin,thetaMax,thetaStep},{Hue[0.7,#1,0.8]&,Hue[0,#1,0.8]&}];

PolarGrid[{rMin_,rMax_,rStep_:Automatic},{thetaMin_,thetaMax_,thetaStep_:Automatic},{cf1_,cf2_}]:=
Module[{thetaRange=N[thetaMax-thetaMin],rRange=N[rMax-rMin],thetaStep2,rStep2},
		If[thetaStep===Automatic,thetaStep2 = thetaRange/14,thetaStep2=thetaStep];
		If[rStep===Automatic,rStep2 = rRange/14,rStep2=rStep];
		N@{Table[{cf2[(i - thetaMin)/thetaRange],Line[{rMin {Cos[i],Sin[i]},rMax {Cos[i],Sin[i]}}]},{i,thetaMin,thetaMax,N[thetaStep2]}],Table[{cf1[(i - rMin)/rRange],Circle[{0,0},i,{thetaMin,thetaMax}]},{i,rMin,rMax,N[rStep2]}]}];

RectangularGrid[{xMin_,xMax_,xStep_:Automatic},{yMin_,yMax_,yStep_:Automatic}]:=RectangularGrid[{xMin,xMax,xStep},{yMin,yMax,yStep},{Hue[0,#1,0.8]&,Hue[0.7,#1,0.8]&}];

RectangularGrid[{xMin_,xMax_,xStep_:Automatic},{yMin_,yMax_,yStep_:Automatic},{cf1_,cf2_}]:=
Module[{xRange=N[xMax-xMin],yRange=N[yMax-yMin],xStep2,yStep2},
		If[xStep===Automatic,xStep2=xRange/First@Flatten@{PlotPoints /. Options[Plot3D, N[PlotPoints]]},xStep2=xStep];

		If[yStep===Automatic,yStep2=yRange/Last@Flatten@{PlotPoints /. Options[Plot3D, N[PlotPoints]]},yStep2=yStep];

		N@{Table[{cf1[(i - xMin)/xRange],Line[{{i,yMin},{i,yMax}}]},{i,xMin,xMax,N[xStep2]}],Table[{cf2[(i - yMin)/yRange],Line[{{xMin,i},{xMax,i}}]},{i,yMin,yMax,N[yStep2]}]}
];

TriangularGrid[d_,n_]:=TriangularGrid[d,n,{0,0},{GrayLevel,Hue[0,#1,0.8]&,Hue[0.7,#1,0.8]&}]/;NumberQ[N[d]]&&NumberQ[N[n]];

TriangularGrid[d_,n_,{c1_,c2_}]:=TriangularGrid[d,n,{c1,c2},{GrayLevel,Hue[0,#1,0.8]&,Hue[0.7,#1,0.8]&}]/;And@@NumberQ/@N[{d,n,c1,c2}];

TriangularGrid[d_,n_,{cf1_,cf2_,cf3_}]:=TriangularGrid[d,n,{0,0},{cf1,cf2,cf3}]/;NumberQ[N[d]]&&NumberQ[N[n]];

TriangularGrid[d_,n_,{c1_,c2_},{cf1_,cf2_,cf3_}]:=Module[{colorList,points1,cos=N[Cos[Pi/3]],sin=N[Sin[Pi/3]]},colorList=N[Transpose[({cf1[#1],cf2[#1],cf3[#1]}&)/@N[Range[0,1,1/(n-1)]]]];points1=Line/@N[Table[{{-((2 d-Abs[i])/Sqrt[3.]),i},{(2 d-Abs[i])/Sqrt[3.],i}},{i,N[-d],N[d],N[(2 d)/(n-1)]}]];points1={Transpose[{colorList[[1]],points1}],Transpose[{colorList[[2]],points1/.line:Line[_]:>Map[{{cos,-sin},{sin,cos}} . #1&,line,{2}]}],Transpose[{colorList[[3]],points1/.line:Line[_]:>Map[{{cos,sin},{-sin,cos}} . #1&,line,{2}]}]};N[If[N[{c1,c2}]===N[{0,0}],points1,points1=points1/.line:Line[_]:>Map[N[#1+{c1,c2}]&,line,{2}]]]]/;And@@NumberQ/@N[{d,n,c1,c2}];

RectangularPatchGraphics[{xMin_, xMax_, xNumber_,xGapScale2_:Automatic}, {yMin_,yMax_,yNumber_,yGapScale2_:Automatic},cf2_:Automatic]:=Module[{xStep,xGapWidth,xSideLength,yStep,yGapWidth,ySideLength,gp,maxLength},
If[cf2===Automatic,cf=(Hue[(#1-xMin)/(xMax-xMin),1-(#2-yMin)/(yMax-yMin),.9]&),cf=cf2];
If[xGapScale2===Automatic,xGapScale=1/10,xGapScale=xGapScale2];
If[yGapScale2===Automatic,yGapScale=1/10,yGapScale=yGapScale2];
xStep=(xMax-xMin)/(xNumber);
yStep=(yMax-yMin)/(yNumber);
xGapWidth=xStep*xGapScale;
yGapWidth=yStep*yGapScale;
xSideLength=xStep-xGapWidth;
ySideLength=yStep-yGapWidth;
gp=Table[Rectangle[{x-xSideLength/2,y-ySideLength/2},{x+xSideLength/2,y+ySideLength/2}],{x, xMin + xStep/2, xMax - xStep/2 + xStep/4, xStep},{y,yMin+yStep/2,yMax- yStep/2 + yStep/4,yStep}];
Graphics@gp/.rect:Rectangle[{a_,b_},{c_,d_}]:>{cf[(a+c)/2,(b+d)/2],rect}
]

PolarPatchGraphics[{rMin_,rMax_,rNumber_,rGapScale2_:Automatic},{thetaMin_, thetaMax_, thetaNumber_,thetaGapScale2_:Automatic},resolution2_:Automatic,cf2_:Automatic]:=Module[{rStep,rGapSweep,rSideSweep,thetaStep,thetaGapSweep,thetaSideSweep,radialBeginEndPairs,sweepBeginEndPairs,sectorsArray,gp,center,sector,resolution},

If[cf2===Automatic,cf=(Hue[#2/(thetaMax-thetaMin),1-(#1-rMin)/(rMax-rMin),.9]&),cf=cf2];If[rGapScale2===Automatic,rGapScale=1/10,rGapScale=rGapScale2];If[thetaGapScale2===Automatic,thetaGapScale=1/10,thetaGapScale=thetaGapScale2];
If[resolution2===Automatic,resolution=1/10,resolution=resolution2];

rStep=(rMax-rMin)/rNumber;
rGapSweep=rStep*rGapScale;
rSideSweep=rStep-rGapSweep;
thetaStep=(thetaMax-thetaMin)/thetaNumber;
thetaGapSweep=thetaStep*thetaGapScale;
thetaSideSweep=thetaStep-thetaGapSweep;

radialBeginEndPairs=Table[{r - rSideSweep/2,r + rSideSweep/2},{r,rMin + rStep/2, rMax - rStep/2 + rStep/4, rStep}];
sweepBeginEndPairs=Table[{theta - thetaSideSweep/2,theta + thetaSideSweep/2},{theta,thetaMin + thetaStep/2, thetaMax - thetaStep/2 + thetaStep/4, thetaStep}];

sectorsArray=Function[{s},sector[#,s]&/@radialBeginEndPairs]/@sweepBeginEndPairs;

gp=sectorsArray/.{sector[{rA_,rB_},{tA_,tB_}]:>sector[center[(rA+rB)/2,(tA+tB)/2],

Polygon@Flatten[{({Cos@#,Sin@#}rA&)/@Range[tA,tB,(tB-tA)/(Ceiling[(rA (-tA+tB))/resolution])],

({Cos@#,Sin@#}rB&)/@Range[tB,tA,(tA-tB)/(Ceiling[(rA (-tA+tB))/resolution])]},1]

]};
Graphics@{gp/.{sector->List,center[r_,theta_]:>cf[r,theta]}}

]

vectorFunctionCanonicalForm[fun_]:=Module[{fun2},
Off[First::"normal",Last::"normal",Function::"slotn",Function::"fpct"];
fun2=Which[((VectorQ@#)&&(And@@(NumberQ/@#))&)@N@fun[3.,4.],fun@@#&,
NumberQ@N@fun[3.,4.]&&Not@NumberQ@N@fun[3.],(({Re[#],Im[#]}&)@(fun@@#))&,
NumberQ@N@fun[3.],({Re[#],Im[#]}&)@fun[First@#+I*Last@#]&,
((VectorQ@#)&&(And@@(NumberQ/@#))&)@N@fun[{3.,4.}],fun@#&];
On[First::"normal",Last::"normal",Function::"slotn",Function::"fpct"];
fun2
]

Options[Transform2DPlot]=Union[{ColorFunction->(Hue[0,#1,0.8]&),PlotPoints->{24,24}},
  Options[Graphics]/.{(AspectRatio->x_):>AspectRatio->Automatic,(Axes->False):>Axes->True}];

Transform2DPlot[expr_,{x_,xMin_,xMax_},{y_,yMin_,yMax_},opts___Rule]:=Transform2DPlot[Function[{x,y},expr],{xMin,xMax},{yMin,yMax},opts];
Transform2DPlot[fun_,{xMin_,xMax_},{yMin_,yMax_},opts___Rule]:=Module[{plotPoints,colorFunction,coordArray,coordArray2,colorInfoArray,minColor,maxColor,vLineGP,hLineGP},plotPoints=PlotPoints/. {opts}/. Options[Transform2DPlot];
plotPoints={First[Flatten[{plotPoints}]],Last[Flatten[{plotPoints}]]};
colorFunction=ColorFunction/. {opts}/. Options[Transform2DPlot];
coordArray=Table[{x,y},{x,N[xMin],N[xMax],N[(xMax-xMin)/(First[plotPoints]-1)]},{y,N[yMin],N[yMax],N[(yMax-yMin)/(Last[plotPoints]-1)]}];
coordArray2=Map[vectorFunctionCanonicalForm[fun],coordArray,{-2}];
colorInfoArray=MapThread[Sqrt[Plus@@((#1-#2)^2)]&,{coordArray,coordArray2},2];
minColor=Min[colorInfoArray];
maxColor=Max[colorInfoArray];
colorInfoArray=Map[colorFunction,(colorInfoArray-minColor)/(maxColor-minColor+0.0001),{-1}];
vLine=Map[Line,(Partition[#1,2,1]&)/@coordArray2,{-3}];
hLine=Map[Line,(Partition[#1,2,1]&)/@Transpose[coordArray2],{-3}];
vLine=Transpose[{Transpose[Drop[Transpose[colorInfoArray],-1]],vLine},{3,2,1}];
hLine=Transpose[{Transpose[Drop[colorInfoArray,-1]],hLine},{3,2,1}];
Show[Graphics[{vLine,hLine},Sequence@@FilterRules[Flatten[{opts,Options[Transform2DPlot]}],Options[Graphics]]]]]/;
And@@NumberQ/@N[{xMin,xMax,yMin,yMax}]

Options[Transform2DGraphicsPlot]=Union[{ResolutionLength ->Automatic,ResolutionDefaultDivision ->50,AdaptiveSampling->True,MaxRecursion->10}, Options[Graphics]];

Transform2DGraphicsPlot[gra_,fun_,opts___Rule]:=Module[{gp,maxLength,rdDivision,as,fun2,maxRecurse},maxLength=ResolutionLength/. {opts}/. Options[Transform2DGraphicsPlot];rdDivision=ResolutionDefaultDivision/. {opts}/. Options[Transform2DGraphicsPlot];as=AdaptiveSampling/. {opts}/. Options[Transform2DGraphicsPlot];maxRecurse=MaxRecursion/. {opts}/. Options[Transform2DGraphicsPlot];fun2=vectorFunctionCanonicalForm[fun];If[Head[gra]===Graphics,gp=N[First[gra]],gp=N[gra]];If[maxLength===Automatic,maxLength=1/rdDivision Max[(Last[#1]-First[#1]&)/@graphicsBoundary[gp]]];(Show[Graphics[#1,Sequence@@FilterRules[Flatten[{opts,Options[Transform2DGraphicsPlot]}],Options[Graphics]]]]&)[If[as,IncreaseGraphicsResolution[gp,maxLength,fun2,maxRecurse],IncreaseGraphicsResolution[gp,maxLength]]/. {poly:Polygon[_]:>Map[fun2,poly,{2}],li:Line[_]:>Map[fun2,li,{2}],pt:Point[_]:>Map[fun2,pt,{1}]}]]/;If[(ResolutionLength/. {opts}/. Options[Transform2DGraphicsPlot])=!=Automatic&&N[ResolutionLength/. {opts}/. Options[Transform2DGraphicsPlot]]==0,Message[ResolutionLength::"zero",fun];False,True]

StereographicProject[{p1_,p2_}]:=With[{cR1=1+p1^2+p2^2},{p1/cR1,p2/cR1,1-1/cR1}]

StereographicProject[{p1_,p2_},r_]:=With[{cR1=4 r^2+4 Sqrt[r^4],cR2=p1^2+p2^2+4 r^2},{(cR1 p1)/(2 cR2),(cR1 p2)/(2 cR2),2 r-(cR1 r)/cR2}]

StereographicProject[{p1_,p2_},1/2,1/2]:=StereographicProject[{p1,p2}]

StereographicProject[{p1_, p2_}, r_, c3_] :=
   With[{n3 = r + c3},
     With[{cR1 = Sqrt[4*(c3 - n3)^2*n3^2 -
               4*(n3^2 + p1^2 + p2^2)*(c3^2 -
                    2*c3*n3 + n3^2 - r^2)],
         cR2 = n3^2 + p1^2 + p2^2},
       {((cR1 - 2*c3*n3 + 2*n3^2)*p1)/(2*cR2),
         ((cR1 - 2*c3*n3 + 2*n3^2)*p2)/(2*cR2),
         n3 - (n3*(cR1 - 2*c3*n3 + 2*n3^2))/
             (2*cR2)}]]

Options[StereographicProjectionPlot3D] =
     Union[{ResolutionLength -> 0.07,
         AdaptiveSampling -> True, MaxRecursion ->
           10}, Options[Graphics3D] /.
         {(Axes -> False) -> Axes -> True,
           (Lighting -> Automatic) -> Lighting ->
               "Neutral"}];

StereographicProjectionPlot3D[gra_Graphics,opts___Rule]:=StereographicProjectionPlot3D[gra,1/2,1/2,opts];

StereographicProjectionPlot3D[gra_Graphics,r2_/;NumberQ[r2],opts___Rule]:=StereographicProjectionPlot3D[gra,r2,1/2,opts];

StereographicProjectionPlot3D[gra_Graphics,r_,c_,opts___Rule]:=Module[{gp},maxLength=ResolutionLength/. {opts}/. Options[StereographicProjectionPlot3D];as=AdaptiveSampling/. {opts}/. Options[Transform2DGraphicsPlot];Show[Graphics3D[If[as,IncreaseGraphicsResolution[List@@gra,maxLength,StereographicProject[#1,r,c]&,maxRecurse],IncreaseGraphicsResolution[List@@gra,maxLength]]/. {p_Polygon:>Map[StereographicProject[#1,r,c]&,p,{2}],li_Line:>Map[StereographicProject[#1,r,c]&,li,{2}]},Sequence@@FilterRules[Flatten[{opts,Options[StereographicProjectionPlot3D]}],Options[Graphics3D]]]]]/;If[(ResolutionLength/. {opts}/. Options[Transform2DGraphicsPlot])=!=Automatic&&N[ResolutionLength/. {opts}/. Options[Transform2DGraphicsPlot]]==0,Message[ResolutionLength::"zero",fun];False,True]&&NumberQ[r]&&NumberQ[c]

MobiusTransform[a_Complex,b_Complex,c_Complex,d_Complex]=Function[(a*#+b)/(c*#+d)];

complexPlus[{a_,b_},{c_,d_}]:={a+c,b+d}
complexTimes[{a_,b_},{c_,d_}]:={a*c-b*d,b*c+a*d}
complexInversion[{a_,b_}]:=1/(a^2+b^2)*{a,-b}

MobiusTransform2[{a1_,a2_},{b1_,b2_},{c1_,c2_},{d1_,d2_}][{z1_,z2_}]:=complexTimes[complexPlus[complexTimes[{a1,a2},{z1,z2}],{b1,b2}],complexInversion[complexPlus[complexTimes[{c1,c2},{z1,z2}],{d1,d2}]]]

End[];



EndPackage[];



BeginPackage["PlaneCurvePlot`"];

Begin[ "Private`" ]

Clear[ normalLineGP ];

normalLineGP = Function[{tPts, tangentVecs, lnlen},
Function[{norms},  Line /@ Transpose[{tPts - norms, tPts + norms}] ][Map[ Function[ Function[{-Last[#], First[#]}][ Normalize[ # ] lnlen/2. ] ] , tangentVecs] ]
];

Clear[ tangentLineGP ];

tangentLineGP = Function[{tPts, tangentVecs, lnlen},
Function[{tans},
Line /@ Transpose[{tPts - tans, tPts + tans}]]
[ Map[ Function[ Normalize[ # ] lnlen/2. ] , tangentVecs] ]
];

parallelDotGP = Function[{cuvPts, velociVectors, distList},
With[{xnomals = Map[ Function[ Function[{-Last[#], First[#]}][Normalize[ # ]] ], velociVectors]},
Map[Point, Map[ Function[ cuvPts + xnomals*#1 ] , distList], {2}]
] ];

Clear[ reflDirVector ];

reflDirVector = Function[{rayV, wallV},
If[Norm[ wallV ] == 0 , {0, 0},
Normalize[- rayV + (2*( wallV . rayV * wallV))/(wallV . wallV) ]]
] ;

Clear[ causticLineGP ];

causticLineGP = Function[{tPts, tangentVecs, rayOrigin, lnlen},
With[ { xsign = Sign[lnlen], lnlen2 = Abs[lnlen]},
If[ Norm[ rayOrigin ] == 0,
MapThread[
If[Norm[#1] >= lnlen2,
Line[{{0,0}, Normalize[#1]*lnlen2}],
Line[{{0,0}, #1, #1 + reflDirVector[#1, #2]* (lnlen2 - Norm[#1])* xsign}]
] &, {tPts, tangentVecs}]
,
MapThread[
If[Norm[#1-rayOrigin] >= lnlen2,
Line[{rayOrigin , rayOrigin + Normalize[#1-rayOrigin]*lnlen2}],
Line[{rayOrigin , #1, #1 + reflDirVector[#1-rayOrigin, #2]* (lnlen2 - Norm[#1-rayOrigin])* xsign}]
] &, {tPts, tangentVecs}]
] ] ];

Clear[ secantLineGP ];

secantLineGP = Function[{tPts, m, n}, Map[ Function[(Line[{First[#1], Last[#1]}] )], Partition[tPts, m + 1, n]]];

Clear[ conchoidLineAndDotGP ];

conchoidLineAndDotGP = Function[{pointP, tPts, distList},
With[{ dirVList = Map[ Function[Normalize[#1 - pointP]], N[tPts]] },
{Line /@ Transpose[{
tPts + Min[distList]*dirVList,
tPts + Max[distList]*dirVList}],
Map[Point, Function[tPts + dirVList*#1 ] /@ distList, {2}]}] ];

Clear[ getPedalPoint ];

getPedalPoint = Function[{cT, tp, p},
Function[If[ SameQ[# == {0, 0}, True ] , {}, tp + ((p - tp) . #) * #]][Normalize[cT]]
 ];

Clear[ pedalLineAndDotGP ];

pedalLineAndDotGP = Function[{pointP, tPts, tangentVecs},
Function[{pedalPoints},
{MapThread[If[#2 === {}, {}, Line[{#1, #2, pointP}]] & , {tPts, pedalPoints}],
(If[#1 === {}, {}, Point[#1]] & ) /@ pedalPoints}
][ MapThread[getPedalPoint[#1, #2, pointP] & , {tangentVecs, tPts}]]
];

Clear[ inversionDotGP ];

inversionDotGP = Function[{center, radius, tPts},
If[ Norm[center] == 0 ,
Point /@
Map[ Function[ If[Norm[#1] == 0, Nothing, (radius^2*Normalize[#1])/ Norm[#1]]] , tPts] ,
Point /@
Map[ Function[ If[Norm[#1 - center] == 0, Nothing, (radius^2*Normalize[#1 - center])/ Norm[#1 - center] + center  ]] , tPts]
]
];

Clear[involuteLineAndDotGP];

involuteLineAndDotGP = Function[{fx, fy, tList, curvePoints, velociVectors},
 Module[{a, b, t, nArcLenF, f, arcLengths, involutePoints},
 {a, b} = N[{First[tList], Last[tList]}];
nArcLenF = NDSolve[{Derivative[1][f][t] == Sqrt[Derivative[1][fx][t]^2 + Derivative[1][fy][t]^2], f[a] == 0}, f, {t, a, b}][[1,1,2]];
 arcLengths = nArcLenF /@ tList;
involutePoints = arcLengths*Normalize /@ (-velociVectors) + curvePoints;
{MapThread[Line[{#1, #2}] & , {curvePoints, involutePoints}], Point /@ involutePoints}]
];

Clear[ nCurvatureK ];

nCurvatureK = Function[
MapThread[
Function[ With[{xnorm=Norm[ #1 ]}, If[ xnorm == 0 , 0, (#2 . {-Last[#1], First[#1]})/ xnorm^3, 0] ] ] ,
{#1, #2}] ];

Clear[ curvatureGPs ];

curvatureGPs = Function[{tPts, tangentVs, acceVs, circleEvoluteRadial },
Function[{pts2, tanV2, accV2},
With[ {
normalUnitVList = Map[ Function[ Function[{-Last[#], First[#]}][Normalize[ # ]] ], tanV2],
radiusList = 1/nCurvatureK[ tanV2, accV2]
},
Block[ { centerVectors, centerPoints },
centerVectors = normalUnitVList*radiusList;
centerPoints = centerVectors + pts2;
{
If[circleEvoluteRadial[[1]], MapThread[Circle[#1, #2] & , {centerPoints, Abs[radiusList]}], {}],
If[circleEvoluteRadial[[2]], {MapThread[Line[{#1, #2}] & , {pts2, centerPoints}], Point /@ centerPoints}, {{}, {}}],
If[circleEvoluteRadial[[3]], {(Line[{{0, 0}, #1}] & ) /@ centerVectors, Point /@ centerVectors} , {{}, {}}]}
 ] ] ] @@
(Function[ (Delete@Position[ ((Norm@#) &) /@ #2 , 0.|0]) /@ {#1,#2,#3} ] [tPts, tangentVs, acceVs])
];

Clear[ curvatureNewGP ];

curvatureNewGP[fx_,fy_, tlist_, tPts_, tVelocities_, {cir_, evo_, rad_}] := Module[
{normalUnitVList, radiusList, centerVectors, centerPoints},
 normalUnitVList = Map[ Function[ Function[{-Last[#], First[#]}][Normalize[ # ]] ], tVelocities];
    radiusList = 1/ Table[ Evaluate@ FrenetSerretSystem[ {fx[t833],fy[t833]}, t833 ][[1,1]], {t833, tlist} ] ;
    centerVectors = normalUnitVList*radiusList;
    centerPoints = centerVectors + tPts;
    {If[cir, MapThread[Circle[#1, #2] & , {centerPoints, Abs[radiusList]}], {}],
     If[evo, {MapThread[Line[{#1, #2}] & , {tPts, centerPoints}], Point /@ centerPoints}, {{}, {}}],
     If[rad, {(Line[{{0, 0}, #1}] & ) /@ centerVectors, Point /@ centerVectors}, {{}, {}}]}];

Clear[ injectStyle ];

injectStyle[graPrims_, xstyle_] :=
Which[
Head[xstyle] =!= List, {xstyle, graPrims},
xstyle === {}, graPrims,
And[ Head[xstyle] === List, Length[xstyle] === 1 ], {Sequence @@ Flatten[{xstyle}] , graPrims},
True, Flatten /@ Transpose[
     {Take[Join @@ Table[xstyle, {Ceiling[Length[graPrims]/Length[xstyle]]}],
       Length[graPrims]], graPrims}]];

Options[PlaneCurvePlot] = Join[{
"PlotCurve" -> True,
"PlotDot" -> True,
"PlotParameterValue" -> False,
"NormalLineLength" -> 0,
"TangentLineLength" -> 0,
"CausticLineLength" -> 0,
"CausticOrigin" -> {0, 0},
"SecantConnection" -> {1, 1},
"ParallelCurveOffset" -> {0},
"ConchoidSetting" -> {{0, 0}, {0}},
"PedalPoint" -> None,
"InversionCircle" -> None,
"PlotOsculatingCircle" -> False,
"PlotInvolute" -> False,
"PlotEvolute" -> False,
"PlotRadial" -> False,
"CurveColorFunction" -> Function[Hue[0.65, #1, 0.8]],
"CurveStyle" -> {Thickness[0.005]},
"DotStyle" -> {{Red, PointSize[0.01]}},
"LineStyle" -> {{Green, Thickness[0.004]}},
"NormalLineStyle" -> {{Hue[0.6], Thickness[0.004]}},
"TangentLineStyle" -> {{Red, Thickness[0.004]}},
"SecantLineStyle" -> {{Red, Thickness[0.004]}},
"CausticLineStyle" -> {{Hue[0.5], Thickness[0.004]}}
}, Options[ParametricPlot]];

PlaneCurvePlot[{xf_Function, yf_Function}, {tmin_, tmax_, dt_:Automatic}, tlist_List:{Automatic}, xOpts:OptionsPattern[]] :=

With[{
xpedalPt = OptionValue["PedalPoint"],
invCircle = OptionValue["InversionCircle"],
xdotstyl = OptionValue["DotStyle"],
xlineStyle = OptionValue["LineStyle"],
tminN = N[tmin],
tmaxN = N[tmax],
dtN = If[dt === Automatic, N[tmax - tmin]/4., N[dt]]
},

Module[{
tVals,
tPoints,
tangentVecs,
acclVecs,
tDotGP,
tTextGP,
xnormalGP,
xtanlineGP,
xcausticLineGP,
xsecantLineGP,
xparallPointGP,
xconcoidLineGP,
xconcoidPointGP,
xpedalLineGP,
xpedalPointGP,
xinversPointGP,
xinvoluteLineGP,
xinvolutePointGP,
xoscCircleGP,
xevoluteLineGP,
xevolutePointGP,
xradialLineGP,
xradialPointGP,
xcurveGP},

tVals = If[tlist === {Automatic}, N@ Range[tmin, tmax, dtN], N[tlist]];
tPoints = Map[ Function[{xf[#1], yf[#1]}], tVals];

tangentVecs = If[
Or[
OptionValue["NormalLineLength"] != 0,
OptionValue["TangentLineLength"] != 0,
OptionValue["CausticLineLength"] != 0,
OptionValue["ParallelCurveOffset"] != {0},
OptionValue["PedalPoint"] =!= None,
OptionValue["PlotInvolute"],
OptionValue["PlotOsculatingCircle"],
OptionValue["PlotEvolute"],
OptionValue["PlotRadial"]
] ,
Map[Function[{Derivative[1][xf][#1], Derivative[1][yf][#1]} ], tVals],
 {} ];

acclVecs =
If[ Or[
OptionValue["PlotOsculatingCircle"],
OptionValue["PlotEvolute"],
OptionValue["PlotRadial"]] ,
Map[Function[{Derivative[2][xf][#1], Derivative[2][yf][#1]}], tVals],
{}
 ];

tTextGP = If[OptionValue@ "PlotParameterValue" === True, Table[Text[NumberForm[tVals[[t]], 4], tPoints[[t]], {-2, -1}], {t, 1, Length[tVals]}], {}];

tDotGP = If[OptionValue@ "PlotDot" === True, Point /@ tPoints, {}];
xnormalGP = If[OptionValue["NormalLineLength"] === 0, {}, normalLineGP[tPoints, tangentVecs, OptionValue["NormalLineLength"]]];

xtanlineGP = If[OptionValue["TangentLineLength"] === 0, {},
 tangentLineGP[tPoints, tangentVecs, OptionValue["TangentLineLength"]]];

xcausticLineGP = If[OptionValue["CausticLineLength"] === 0, {}, causticLineGP[tPoints, tangentVecs, N[OptionValue["CausticOrigin"]], N[OptionValue["CausticLineLength"]]]];

xsecantLineGP = If[OptionValue["SecantConnection"] === {1, 1}, {}, secantLineGP[tPoints, First[OptionValue["SecantConnection"]], Last[OptionValue["SecantConnection"]]]];

xparallPointGP = If[OptionValue["ParallelCurveOffset"] === {0}, {}, parallelDotGP[tPoints, tangentVecs, OptionValue["ParallelCurveOffset"]]];

 {xconcoidLineGP, xconcoidPointGP} = If[Last[OptionValue["ConchoidSetting"]] === {0}, {{}, {}}, conchoidLineAndDotGP[ First[OptionValue["ConchoidSetting"]], tPoints, Last[OptionValue["ConchoidSetting"]]]];

{xpedalLineGP, xpedalPointGP} = If[xpedalPt === None, {{}, {}}, pedalLineAndDotGP[xpedalPt, tPoints, tangentVecs]];

xinversPointGP = If[invCircle === None, {}, inversionDotGP[First@invCircle, Last@ invCircle, tPoints]];

{xinvoluteLineGP, xinvolutePointGP} = If[OptionValue["PlotInvolute"] === False, {{}, {}}, involuteLineAndDotGP[xf, yf, tVals, tPoints, tangentVecs]];

{xoscCircleGP, {xevoluteLineGP, xevolutePointGP}, {xradialLineGP, xradialPointGP}} =
If[acclVecs === {},
{{}, {{}, {}}, {{}, {}}},
curvatureGPs[tPoints, tangentVecs, acclVecs, {OptionValue["PlotOsculatingCircle"], OptionValue["PlotEvolute"], OptionValue["PlotRadial"]}]
];

xcurveGP = If[OptionValue[ "PlotCurve" ] === True,
First@ ParametricPlot[ {xf[t], yf[t]}, {t, tminN, tmaxN},
DisplayFunction -> Identity,
ColorFunction :> Function[ First[ Flatten[ {OptionValue["CurveColorFunction"]} ] ][#3] ],
PlotStyle :> OptionValue["CurveStyle"],
Evaluate[FilterRules[{xOpts}, Options[ ParametricPlot]]]
], {}];

Show[Graphics[ {
injectStyle[xradialLineGP, xlineStyle],
injectStyle[xoscCircleGP, xlineStyle],
injectStyle[xinvoluteLineGP, xlineStyle],
injectStyle[xevoluteLineGP, xlineStyle],
injectStyle[xtanlineGP , OptionValue["TangentLineStyle"]],
injectStyle[xnormalGP, OptionValue["NormalLineStyle"]] ,
injectStyle[xsecantLineGP, OptionValue["SecantLineStyle"]],
injectStyle[xcausticLineGP, OptionValue["CausticLineStyle"]],
injectStyle[xconcoidLineGP, xlineStyle],
injectStyle[xpedalLineGP, xlineStyle],
xcurveGP,
injectStyle[xinvolutePointGP, xdotstyl],
injectStyle[xevolutePointGP, xdotstyl],
injectStyle[xradialPointGP, xdotstyl],
injectStyle[xparallPointGP, xdotstyl],
injectStyle[xconcoidPointGP, xdotstyl],
injectStyle[xpedalPointGP, xdotstyl],

If[invCircle === None, {}, {Yellow, Circle @@ invCircle}],
injectStyle[xinversPointGP, xdotstyl],
injectStyle[tDotGP, xdotstyl],
tTextGP,
If[OptionValue["CausticLineLength"] === 0, {}, {Red, PointSize[0.02], Point[OptionValue["CausticOrigin"]]}],
If[Last[OptionValue["ConchoidSetting"]] === {0}, {}, {Blue, PointSize[0.02], Point[First[OptionValue["ConchoidSetting"]]]}],
If[xpedalPt === None, {}, {Blue, PointSize[0.02], Point[xpedalPt]}]
} ,
FilterRules[{xOpts}, Options[Graphics]], Axes -> True]]
] ] /;
And[NumericQ[tmin], NumericQ[tmax], Or[NumericQ[dt], SameQ[dt, Automatic]]];

PlaneCurvePlot[{fx_,fy_}, {t838_Symbol, tmin_, tmax_, dt_:Automatic}, tlist_List:{Automatic}, xOpts:OptionsPattern[]] :=
PlaneCurvePlot[ {Function[t838, fx ], Function[t838, fy]}, {tmin, tmax, dt}, tlist, xOpts] /;
And[NumericQ[tmin], NumericQ[tmax], Or[NumericQ[dt], SameQ[dt, Automatic]]];

End[];

EndPackage[];



BeginPackage["PlaneTiling`"]

Begin["`Private`"]

Clear[LatticeCoordinates];

LatticeCoordinates::"usage" =
"LatticeCoordinates[vectorA,vectorB,{m,n}] return a list of coordinates on a grid generated by vectors vectorA and vectorB with m and n points in each direction.
LatticeCoordinates[vectorA,vectorB,{m,n},s,alpha,{x,y}] scales, rotates, and translates the grid by s, alpha, and {x,y}.
Example:
Graphics[{Map[Point,LatticeCoordinates[{1,0},{1/2,Sqrt[3]/2},{7,5}],{2}]}, Frame -> True]
";

LatticeCoordinates[a_,b_,{m_,n_}]:=
Table[a*i+b*j,{j,0,n-1},{i,0,m-1}];

LatticeCoordinates[a_,b_,{m_,n_},s_]:=
Table[a*s*i+b*s*j,{j,0,n-1},{i,0,m-1}];

LatticeCoordinates[a_,b_,{m_,n_},s_,alpha_]:=
Function[
LatticeCoordinates[#1,#2,{m,n}]
]@@
((({{Cos@alpha,-Sin@alpha},{Sin@alpha,Cos@alpha}}) . #)&/@({a,b}*s))

LatticeCoordinates[a_,b_,{m_,n_},s_,alpha_,{x_,y_}]:=
Map[(#+{x,y})&,LatticeCoordinates[a,b,{m,n},s,alpha],{2}];

Clear[SquareLatticeCoordinates];

SquareLatticeCoordinates::"usage" =
"SquareLatticeCoordinates[n] return a list of coordinates that represents a square lattice.
n is the number of points in one dimention.
SquareLatticeCoordinates[n,s,alpha,{x,y}] scales, rotates, and translates the lattice by s, alpha, and {x,y}.
Example:
Graphics[Map[Point,SquareLatticeCoordinates[5],{2}], Frame -> True]
";

SquareLatticeCoordinates[n_ : 5] :=
 LatticeCoordinates[{1, 0}, {0, 1}, {n, n}]

SquareLatticeCoordinates[n_, s_] :=
 LatticeCoordinates[{1, 0}, {0, 1}, {n, n}, s]

SquareLatticeCoordinates[n_, s_, alpha_] :=
 LatticeCoordinates[{1, 0}, {0, 1}, {n, n}, s, alpha]

SquareLatticeCoordinates[n_, s_, alpha_, {x_, y_}] :=
 LatticeCoordinates[{1, 0}, {0, 1}, {n, n}, s, alpha, {x, y}]

Clear[TriangleLatticeCoordinates];

TriangleLatticeCoordinates[n_Integer] := TriangleLatticeCoordinates[n, 1, 0]

TriangleLatticeCoordinates[n_] := TriangleLatticeCoordinates[n, 1., 0]

TriangleLatticeCoordinates[n_, s_] := TriangleLatticeCoordinates[n, s, 0]

TriangleLatticeCoordinates[n_, s_, alpha_] :=
(Table[
Map[{{Cos@i, -Sin@i}, {Sin@i, Cos@i}} . # &, #, {2}], {i, alpha, alpha + 2*Pi - Pi/3, 2*Pi/3}] &)@
(LatticeCoordinates[{1, 0}, {Cos[2*Pi/3], Sin[2*Pi/3]}, {n, n}, s, 0, {s/2, Sqrt[3]/6*s}]);

TriangleLatticeCoordinates[n_, s_, alpha_, {x_, y_}] :=
 Map[# + {x, y} &, TriangleLatticeCoordinates[n, s, alpha], {3}]

Clear[StarLatticeCoordinates];

StarLatticeCoordinates[n_Integer] :=
 Map[Table[{Cos@t, Sin@t} + #, {t, 0, Pi - 2*Pi/6/2, 2*Pi/6}] &,
   TriangleLatticeCoordinates[n, 2], {3}]

StarLatticeCoordinates[n_?MachineNumberQ] :=
 Map[Table[{Cos@t, Sin@t} + #, {t, 0, Pi - 2*Pi/6/2, 2*Pi/6}] &,
   TriangleLatticeCoordinates[n, 2.],
 {3}]

StarLatticeCoordinates[n_, s_] :=
 Map[(#*s) &, StarLatticeCoordinates[n], {4}]

StarLatticeCoordinates[n_, s_, alpha_] :=
 Map[(({{Cos@alpha, -Sin@alpha}, {Sin@alpha, Cos@alpha}}*
        s) . #) &, StarLatticeCoordinates[n], {4}]

StarLatticeCoordinates[n_, s_, alpha_, {x_, y_}] :=
 Map[(({{Cos@alpha, -Sin@alpha}, {Sin@alpha, Cos@alpha}}*
         s) . # + {x, y}) &, StarLatticeCoordinates[n], {4}]

Clear[HexagonLatticeCoordinates];

HexagonLatticeCoordinates[n_ : 5] :=
Map[
Table[{Cos@t, Sin@t} + #, 
{t, 2*Pi/12, 2*Pi/12 + Pi - 2*Pi/6/2, 2*Pi/6}] &, 
TriangleLatticeCoordinates[n, Sqrt[3], 0], {3}]

HexagonLatticeCoordinates[n_, s_] :=
Map[(#*s) &, HexagonLatticeCoordinates[n], {4}]

HexagonLatticeCoordinates[n_, s_, alpha_] :=
 Map[
(({{Cos@alpha, -Sin@alpha}, {Sin@alpha, Cos@alpha}}* s) . #) &, 
HexagonLatticeCoordinates[n], {4}]

HexagonLatticeCoordinates[n_, s_, alpha_, {x_, y_}] :=
 Map[
(({{Cos@alpha, -Sin@alpha}, {Sin@alpha, Cos@alpha}}* s) . # + {x, y}) &, 
HexagonLatticeCoordinates[n], {4}]

Clear[cycledMatrix];

cycledMatrix[{m_Integer, n_Integer}, cycleLists : {_List ..}] :=
(Flatten[({Table[#1, {Quotient[#2, #3]}], Take[#1, Mod[#2, #3]]} &) @@ ({#, m, Length@#} &)@#] &) /@
Flatten[
(({(Flatten[#, 1] &)@Table[#1, {Quotient[#2, #3]}],
Take[#1, Mod[#2, #3]]} &) @@ ({#, n, Length@#} &)@ cycleLists), 1];

Clear[ColoredLatticePoints];

ColoredLatticePoints[a_,b_,{m_Integer,n_Integer},colors:{_List..}]:=
ColoredLatticePoints[a,b,{m,n},colors,1,0,{0,0}];

ColoredLatticePoints[a_,b_,{m_Integer,n_Integer},colors:{_List..},s_]:=
ColoredLatticePoints[a,b,{m,n},colors,s,0,{0,0}];

ColoredLatticePoints[a_,b_,{m_Integer,n_Integer},colors:{_List..},s_,alpha_]:=
ColoredLatticePoints[a,b,{m,n},colors,s,alpha,{0,0}];

ColoredLatticePoints[a_,b_,{m_Integer,n_Integer},colors:{_List..},s_,alpha_,{x_,y_}]:=
(Transpose[{#1,#2},{3,1,2}]&)@@{cycledMatrix[{m,n},colors],(Map[Point,LatticeCoordinates[a,b,{m,n},s,alpha,{x,y}],{2}])};

Clear[ColoredLatticeNetwork];

ColoredLatticeNetwork[n:(4|3),m_Integer,colors:{_List..}]:=
ColoredLatticeNetwork[n,m,colors,1,0,{0,0}];

ColoredLatticeNetwork[n:(4|3),m_Integer,colors:{_List..},s_]:=
ColoredLatticeNetwork[n,m,colors,s,0,{0,0}];

ColoredLatticeNetwork[n:(4|3),m_Integer,colors:{_List..},s_,alpha_]:=
ColoredLatticeNetwork[n,m,colors,s,alpha,{0,0}];

ColoredLatticeNetwork[n:(4|3),m_Integer,colors:{_List..},s_,alpha_,{x_,y_}]:=
Module[{lines,gp},

		lines=(Map[Line,#,{2}]&)@(Partition[#,2,1]&/@N@Table[{1,0}*s*i+(({Cos@#,Sin@#}&)@(2*Pi/n))*s*j,{j,-m,m},{i,-m,m}]);
		gp=(Transpose[#,{3,1,2}]&)@{cycledMatrix[{m*2,m*2+1},(RotateRight[#,m]&)@((RotateRight[#,m]&)/@colors)],lines};

	Translate2DGraphics[{x,y}]@Table[Rotate2DGraphics[{0,0},(2*Pi/n)*i+alpha]@gp,{i,0,n-1}]
		];

Clear[DirectedLatticeNetwork];

DirectedLatticeNetwork[allArguments___]:=
ColoredLatticeNetwork[allArguments]/.{{True,Line[pts_]}->Line@pts,{False,Line[pts_]}:>Line@Reverse@pts};

Clear[PolygonMotif];

(*Notes:
Input outputs:
Output has the form {Line[...],Line[...],...}, where there are GCD[p,q] elements.
The order of lines is their polygon rotated counter-clockwise.
The "sense" of coordinates inside every Line depends on whether q > p/2.
e.g.
in {8,3}, counter-clockwise; in {8,5}, clockwise.
Implementation note:
p2 and q2 are coprime, derived from p,q.
They are equal if p,q are coprime.*)

PolygonMotif[n_Integer,rest1___]:=
PolygonMotif[{n,1},rest1];

PolygonMotif[{p_Integer,q_Integer}]:=
PolygonMotif[{p,q},1,0,{0,0}];

PolygonMotif[{p_Integer,q_Integer},s_]:=
PolygonMotif[{p,q},s,0,{0,0}];

PolygonMotif[{p_Integer,q_Integer},s_,alpha_]:=
PolygonMotif[{p,q},s,alpha,{0,0}];

PolygonMotif[{p_Integer, q_Integer}, s_, alpha_, {x_, y_}] :=
Module[{p2, q2}, 
{p2, q2} = {p, q}/GCD[p, q];
Line /@ Table[{Cos[t + beta], Sin[t + beta]}*s + {x, y}, 
{beta, 0 + alpha, (2*Pi)/p2 + alpha - ((2*Pi)/p)/2, (2*Pi)/p}, 
{t, 0, (2*Pi)*q2, (2*Pi)/(p2/q2)}]];

(*
PolygonMotif
Old code.

In this old version, each Line connects to two points.
Originally this design decision is made because I thought this package will be Line based and accept such limitations.
Later I decided the package will not accept such limitation.
1998/06/01.

PolygonMotif[{p_Integer,q_Integer},s_,alpha_,{x_,y_}]:=
Map[Line,Transpose[{#,RotateLeft[#,q]}]]&@(Table[{Cos[t],Sin[t]}*s+{x,y},{t,0+alpha,2*Pi+alpha-2*Pi/p/2,2*Pi/p}])

*)

Clear[PolygonMotif2];

PolygonMotif2[n_Integer,rest1___]:=
PolygonMotif2[{n,1},rest1];

PolygonMotif2[{p_Integer,q_Integer}]:=
PolygonMotif2[{p,q},1,0,{0,0}];

PolygonMotif2[{p_Integer,q_Integer},s_]:=
PolygonMotif2[{p,q},s,0,{0,0}];

PolygonMotif2[{p_Integer,q_Integer},s_,alpha_]:=
PolygonMotif2[{p,q},s,alpha,{0,0}];

PolygonMotif2[{p_Integer,q_Integer},s_,alpha_,{x_,y_}]:=
Module[{p2,q2},{p2,q2}={p,q}/GCD[p,q];q2=Min[q2,p2-q2];
StarMotif[p,If[q===1,{1},{1, Cos[(Pi*q2)/p2]*Sec[Pi*(1/p - q2/p2)]}],s,alpha,{x,y}]];



Clear[StarMotif];

StarMotif[n_Integer, {}]:={};

StarMotif[n_Integer,li_?VectorQ]:=
StarMotif[n,li,1,0,{0,0}];

StarMotif[n_Integer,li_?VectorQ,s_]:=
StarMotif[n,li,s,0,{0,0}];

StarMotif[n_Integer,li_?VectorQ,s_,alpha_]:=
StarMotif[n,li,s,alpha,{0,0}];

StarMotif[n_Integer,li_?MatrixQ]:=
StarMotif[n,li,1,0,{0,0}];

StarMotif[n_Integer,li_?MatrixQ,s_]:=
StarMotif[n,li,s,0,{0,0}];

StarMotif[n_Integer,li_?MatrixQ,s_,alpha_]:=
StarMotif[n,li,s,alpha,{0,0}];

StarMotif[n_Integer,radii_?VectorQ,s_,alpha_,{x_,y_}]:=
StarMotif[n,Transpose[{radii,Table[2*Pi/n/Length@radii*i,{i,0,Length@radii-1}]}],s,alpha,{x,y}];

StarMotif[n_Integer,coords_?MatrixQ,s_,alpha_,{x_,y_}]:=
Polygon@(Flatten[#,1]&)@Transpose@Map[Table[{Cos[t+Last@#],Sin[t+Last@#]}*First[#]*s+{x,y},{t,alpha,2*Pi+alpha-2*Pi/n/2,2*Pi/n}]&,coords];

Clear[NestedPolygonMotif];

NestedPolygonMotif[n_Integer,m_Integer]:=
NestedPolygonMotif[n,m,1,0,{0,0}];

NestedPolygonMotif[n_Integer,m_Integer,s_]:=
NestedPolygonMotif[n,m,s,0,{0,0}];

NestedPolygonMotif[n_Integer,m_Integer,s_,alpha_]:=
NestedPolygonMotif[n,m,s,alpha,{0,0}];

NestedPolygonMotif[n_Integer, m_Integer, s_, alpha_, {x_, y_}] :=
(Polygon /@
NestList[(Append[#1, First[#1]] &)[(1/2*(First[#1] + Last[#1]) &) /@ Partition[#1, 2, 1]] &,
Table[{Cos[t], Sin[t]}*s + {x, y}, {t, 0 + alpha,  2*Pi + alpha, (2*Pi)/n}], m]) /.
 li_Line :> Drop[li, -1];

Clear[RotationSymbolMotif2];

(*
This is old RotationSymbolMotif, that uses n bars inside a circle to indicate n-fold rotation.
The advantage over traditional rotation symbol is that this maintains consistency over n-fold rotations, while traditional symbol breaks consistency for 2-fold rotation.
The drawback of the new symbol is that it cannot be filled, as to represent rotation centers that's on mirror line.
I decided to go with the traditional symbol.
It looks better in Mathematica.
Mathematica, or rather PostScript, does not render small circles on the screen well.
*)

RotationSymbolMotif2::"usage" =
"RotationSymbolMotif2[n] return a list of graphics primitives that represents n-fold rotation centered on {0,0}.
RotationSymbolMotif2[n,r,alpha,{x,y}] scales, rotates and translates it by r, alpha and {x,y}.
Example:
Graphics[{RotationSymbolMotif2[3]},Frame->True]
";

RotationSymbolMotif2[0]:={};

RotationSymbolMotif2[n_Integer]:=RotationSymbolMotif2[n,1,0,{0,0}];
RotationSymbolMotif2[n_Integer,s_]:=RotationSymbolMotif2[n,s,0,{0,0}];
RotationSymbolMotif2[n_Integer,s_,alpha_]:=RotationSymbolMotif2[n,s,alpha,{0,0}];

RotationSymbolMotif2[n_Integer, s_, alpha_, {x_, y_}] := 
  N@{(Line[{{x, y}, #}]) & /@ (
N@ Table[{Cos[t], Sin[t]}*s + {x, y}, {t, 0 + alpha, 2*Pi + alpha - 2*Pi/n/2, 2*Pi/n}]),
Line@(N@Table[{Cos[t], Sin[t]}*s + {x, y}, {t, 0 + alpha, 2*Pi + alpha, 2*Pi/12}])};

Clear[RotationSymbolMotif];

RotationSymbolMotif::"usage" =
"RotationSymbolMotif[n] return a list of graphics primitives that represents n-fold rotation centered on {0,0}.
RotationSymbolMotif[n,r,alpha,{x,y}] scales, rotates and translates it by r, alpha and {x,y}.
The graphics returned has area 1 for any n.
Example:
Graphics[{RotationSymbolMotif[2]}]
";

(*For a regular n-gon to have area 1, it must have radius Sqrt[2/(Sin[2*Pi/n]*n)]. I got this formula by using law of cosine and Heron's formula on isosceles*)

RotationSymbolMotif[0]:={};
RotationSymbolMotif[1]:={};

RotationSymbolMotif[2]:=With[{s=Sqrt[GoldenRatio/2],h=1/Sqrt[2*GoldenRatio]},Line@{{0,s},{-h,0},{0,-s},{h,0},{0,s}}];

RotationSymbolMotif[3]:=Line@{{0, 2/3^(3/4)},
  {1/3^(1/4), -(1/3^(3/4))},
  {-(1/3^(1/4)), -(1/3^(3/4))},
  {0, 2/3^(3/4)}};

RotationSymbolMotif[4]:=Line@{{1/2, 1/2}, {1/2, -(1/2)}, {-(1/2), -(1/2)}, {-(1/2), 1/2}, {1/2, 1/2}};

RotationSymbolMotif[6]:=Line@{{0, Sqrt[2]/3^(3/4)},
  {1/(Sqrt[2]*3^(1/4)),
   1/(Sqrt[2]*3^(3/4))},
  {1/(Sqrt[2]*3^(1/4)),
   -(1/(Sqrt[2]*3^(3/4)))},
  {0, -(Sqrt[2]/3^(3/4))},
  {-(1/(Sqrt[2]*3^(1/4))),
   -(1/(Sqrt[2]*3^(3/4)))},
  {-(1/(Sqrt[2]*3^(1/4))),
   1/(Sqrt[2]*3^(3/4))},
  {0, Sqrt[2]/3^(3/4)}};

RotationSymbolMotif[n_Integer]:=Line@(Table[{Sin[t],Cos[t]}*Sqrt[2/(Sin[2*Pi/n]*n)],{t,0,2*Pi,2*Pi/n}]);

RotationSymbolMotif[n_Integer,s_]:=RotationSymbolMotif[n,s,0,{0,0}];

RotationSymbolMotif[n_Integer,s_,alpha_]:=RotationSymbolMotif[n,s,alpha,{0,0}];

RotationSymbolMotif[n_Integer,s_,alpha_,{x_,y_}]:=Map[({{Cos[alpha],-Sin[alpha]},{Sin[alpha],Cos[alpha]}} . (#*s)+{x,y})&,RotationSymbolMotif[n],{2}];

Clear[ArrowMotif];

(* 2024-02-22
change ArrowMotif to use builtin Arrow. or simply remove it. check if elsewhere is using it.
 *)

ArrowMotif::"usage" =
"ArrowMotif[{{alpha1,s1}}] return a list of Line that represents an arrow from {0,0} to {1,0} with an arrow head line inclined alpha1 angle and length s1.
ArrowMotif[{{alpha1,s1},{alpha2,s2},...}] uses multiple lines to represent the arrow head.
ArrowMotif[arrowHeadMatix,s,alpha,{x,y}] scales, rotates and translates it by s, alpha and {x,y}.
See also: VectorMotif, DoubleLineMotif.
Example:
Graphics[{ArrowMotif[{{-30*Degree,.1},{30*Degree,.1}},1, 30*Degree,{1,1}]}, Frame->True]
";

ArrowMotif[headList_?MatrixQ]:={Line[{{0,0},{1,0}}],MapThread[Line[{{1,0},{-Cos[#1],Sin[#1]}*#2+{1,0}}]&,Transpose[Partition[Flatten[headList],2]]]};

ArrowMotif[headList_?MatrixQ,s_]:=ArrowMotif[headList,s,0,{0,0}];

ArrowMotif[headList_?MatrixQ,s_,alpha_]:=ArrowMotif[headList,s,alpha,{0,0}];

ArrowMotif[headList_?MatrixQ,s_,alpha_,{x_,y_}]:={Line[{{0,0},{1,0}}],MapThread[Line[{{1,0},{-Cos[#1],Sin[#1]}*#2+{1,0}}]&,Transpose[Partition[Flatten[headList],2]]]}/.li_Line:>Map[(({{Cos@alpha,-Sin@alpha},{Sin@alpha,Cos@alpha}}) . (#*s)+{x,y})&,li,{2}]

Clear[VectorMotif];

VectorMotif::"usage" =
"VectorMotif[{v1,v2}] return a list of Line that represents the vector {v1,v2}.
VectorMotif[{v1,v2},s,alpha,{x,y}] scales, rotates and translates it by s, alpha and {x,y}.
See also: ArrowMotif, DoubleLineMotif.
Example:
Graphics[{VectorMotif[{2,1}]},Frame->True]
";

VectorMotif[{a_,b_}]:=VectorMotif[{a,b},1,0,{0,0}];

VectorMotif[{a_,b_},s_]:=VectorMotif[{a,b},s,0,{0,0}];

VectorMotif[{a_,b_},s_,alpha_]:=VectorMotif[{a,b},s,alpha,{0,0}];

VectorMotif[{a_,b_},s_,alpha_,{x_,y_}]:=
ArrowMotif[{{-30*Degree,.1},{30*Degree,.1}},Sqrt[# . #]&@({a,b}),ArcTan[a,b]+alpha,{x,y}];

Clear[DoubleLineMotif];

DoubleLineMotif::"usage" =
"DoubleLineMotif[{a1,a2},{b1,b2},d] return a pair of parallel Lines from {a1,a2} to {b1,b2} with width d in between.
DoubleLineMotif[d] is equivalent to DoubleLineMotif[{0,0},{1,0},d].
DoubleLineMotif[d,s,alpha,{x,y}] scales, rotates and translates it by s, alpha and {x,y}.
See also: ArrowMotif, VectorMotif.
Example:
Graphics[{DoubleLineMotif[{0,0},{2,1},.1]},Frame->True]
";

DoubleLineMotif[a_?VectorQ,b_?VectorQ,d_]:=DoubleLineMotif[d,Sqrt[(b-a) . (b-a)],ArcTan@@(b-a),a];

DoubleLineMotif[d_]:={Line[{{0,-d/2},{1,-d/2}}],Line[{{0,d/2},{1,d/2}}]};

DoubleLineMotif[d_,s_]:=DoubleLineMotif[d,s,0,{0,0}];

DoubleLineMotif[d_,s_,alpha_]:=DoubleLineMotif[d,s,alpha,{0,0}];

DoubleLineMotif[d_,s_,alpha_,{x_,y_}]:=Map[({{Cos@alpha,-Sin@alpha},{Sin@alpha,Cos@alpha}} . #+{x,y})&,{Line[{{0,-d/2},{s,-d/2}}],Line[{{0,d/2},{s,d/2}}]},{3}];

Clear[RosetteMotif];

RosetteMotif::"usage" =
"RosetteMotif[graphics1,n,mirrorLine] return an n-fold centrally symmetric design based on graphics1.
mirrorLine is either True or False.
graphics1 must be a List of Point (not implemented yet), Line, Polygon (not implemented yet), or a Graphics object.
RosetteMotif[graphics1,n,mirror1,s,alpha,{x,y}] scales, rotates and translates it by s, alpha and {x,y}.
Example:
RosetteMotif[Line[{{1,0},{1,1}}],4,True,1,1,{3,2.53}]
";

RosetteMotif[gra_,n_Integer]:=RosetteMotif[gra,n,False,1,0,{0,0}];

RosetteMotif[gra_,n_Integer,mirror:(True|False)]:=RosetteMotif[gra,n,mirror,1,0,{0,0}];

RosetteMotif[gra_,n_Integer,mirror:(True|False),scale_]:=RosetteMotif[gra,n,mirror,scale,0,{0,0}];

RosetteMotif[gra_,n_Integer,mirror:(True|False),scale_,alpha_]:=RosetteMotif[gra,n,mirror,scale,alpha,{0,0}];

RosetteMotif[gra_,n_Integer,mirror:(True|False),scale_,alpha_,{x_,y_}]:=
	Block[{gra2},
		gra2=If[mirror,gra/.li_Line:>{li,Map[{First@#,-Last@#}&,li,{2}]},gra]/.Line[pt_List]:>Line[pt*scale];
		Table[gra2/.li_Line:>Map[{{Cos[t],Sin[t]},{-Sin[t],Cos[t]}} . #+{x,y}&,li,{2}],{t,alpha,2*Pi+alpha-2*Pi/n/2,2*Pi/n}]
		];

Clear[StainedGlassMotif];

StainedGlassMotif::"usage" =
"StainedGlassMotif[n, sideLength, randomParameter, {i,j}] return a j by i matrix of random Polygon that represents a stained glass.
The pattern is based on regular tiling of triangle, square, or hexagon of sideLength, depending on whether n is 3, 4, or 6.
randomParameter is a real number between 0 and 1 that controls the randomness of the image.
StainedGlassMotif[n,sl,rnd,{{xMin,xMax},{yMin,yMax}}] return graphics that covers a specified rectangle.
Example:
Graphics[{StainedGlassMotif[6,1,.5,{5,4}]/.poly_Polygon:>{RandomColor[], poly}}, AspectRatio->Automatic]
";

StainedGlassMotif[3, s_, r_, {{xMin_, xMax_}, {yMin_, yMax_}}] :=
   Module[{i, j},
    {i, j} = Ceiling[N[{(xMax - xMin)/(s/2), (yMax - yMin)/
           ((Sqrt[3]*s)/2)}]] + 4;
     Translate2DGraphics[(Plus @@ #1/2 & ) /@ {{xMin, xMax},
          {yMin, yMax}} - With[{a = {1, 0}*s*i,
          b = {1/2, Sqrt[3]/2}*s*j}, (a + b)/2 +
          {-(s/2), -((Sqrt[3]*s)/(2*3))}]][StainedGlassMotif[3, s, r,
       {i, j}]]];
StainedGlassMotif[3, s_, r_, {ii_Integer, jj_Integer}] :=
   Module[{pts}, pts = N[Table[{1, 0}*s*i + {1/2, Sqrt[3]/2}*s*j +
         {-(s/2), -((Sqrt[3]*s)/(2*3))}, {j, 0, jj}, {i, 0, ii}]];
     If[N[r] != 0, pts = Map[(1/2)*({Cos[#1], Sin[#1]} & )[
            RandomReal[{0, N[2*Pi]}]]*s*RandomReal[{0, N[r]}] + #1 & ,
        pts, {2}]]; (MapThread[{Polygon[Flatten[{#1, {First[#2]}}, 1]],
          Polygon[Flatten[{{Last[#1]}, #2}, 1]]} & , #1] & ) /@
      Map[Partition[#1, 2, 1] & , Partition[pts, 2, 1], {2}]];

StainedGlassMotif[4, s_, r_, {{xMin_, xMax_}, {yMin_, yMax_}}] :=
   Module[{i, j}, {i, j} = Ceiling[N[{xMax - xMin, yMax - yMin}/s]] +
       4; Translate2DGraphics[(Plus @@ #1/2 & ) /@ {{xMin, xMax},
          {yMin, yMax}} - ((1/2)*(-s + #1*s) & ) /@ {i, j}][
      StainedGlassMotif[4, s, r, {i, j}]]];
StainedGlassMotif[4, s_, r_, {ii_Integer, jj_Integer}] :=
   Module[{pts}, pts = N[Table[{s, 0}*i + {0, s}*j + {-(s/2), -(s/2)},
        {j, 0, jj}, {i, 0, ii}]]; If[N[r] != 0,
      pts = Map[(1/2)*({Cos[#1], Sin[#1]} & )[RandomReal[{0, N[2*Pi]}]]*
           s*RandomReal[{0, N[r]}] + #1 & , pts, {2}]];
     (MapThread[Polygon[Flatten[{#1, Reverse[#2]}, 1]] & , #1] & ) /@
      Map[Partition[#1, 2, 1] & , Partition[pts, 2, 1], {2}]];

StainedGlassMotif[6, s_, r_, {{xMin_, xMax_}, {yMin_, yMax_}}] :=
   Module[{i, j},
    {i, j} = Ceiling[N[{(xMax - xMin)/(Sqrt[3]*s), (yMax - yMin)/
           ((3*s)/2)}]] + 4; Translate2DGraphics[
       (Plus @@ #1/2 & ) /@ {{xMin, xMax}, {yMin, yMax}} -
        {(1/2)*Sqrt[3]*s*i, (3*s*j)/(2*2)}][StainedGlassMotif[6, s, r,
       {i, j}]]];
StainedGlassMotif[6, s_, r_, {ii_Integer, jj_Integer}] :=
   Module[{pts, pts2, gp},
    pts = N[Table[{Sqrt[3]*s, 0}*i + {(1/2)*(-Sqrt[3])*s, -(s/2)},
        {i, 0, ii}]]; pts = (Flatten[#1, 1] & )[
       Transpose[{pts, pts /. {x1_, y1_} :> {x1 + (Sqrt[3]*s)/2,
            y1 - s/2}}]]; pts2 = ({1, -1}*#1 & ) /@ pts;
     pts = (Flatten[#1, 1] & )[Table[Map[{0, 3*s}*i + #1 & ,
         {pts, pts2}, {2}], {i, 0, Ceiling[jj/2]}]];
     pts2 = If[N[r] == 0, pts,
       Map[(1/2)*({Cos[#1], Sin[#1]} & )[RandomReal[{0, N[2*Pi]}]]*s*
           RandomReal[{0, N[r]}] + #1 & , pts, {2}]];
     gp = (MapThread[Polygon[Flatten[{#1, Reverse[#2]}, 1]] & ,
         #1] & ) /@ Map[Partition[#1, 3, 2] & , (Flatten[#1, 1] & )[
         Transpose[({First[#1], Map[Rest, Last[#1], {2}]} & )[
           Transpose[(Partition[#1, 2] & )[Partition[pts2, 2, 1]]]]]],
        {2}]; If[OddQ[jj], Delete[gp, -1], gp]];

Clear[LineTransform];

LineTransform::"usage" =
"LineTransform[graphics1,replacementList] is an L-system like graphical function.
graphics1 is a list of Line or an Graphics object.
replacementList is a List of Line or graphic directives and primitives.
Each Line in graphics1 is replaced by the relation Line[{{0,0},{1,0}}]->replacementList.
LineTransform[graphics1,replacementList,Line[{p1,p2}]] use the relation Line[{p1,p2}]->replacementList.
Example:
LineTransform[Line[{{0,0},{1,0},{1,-1}}],{Line[{{0,0},{1/2,1/2},{1,0}}]}]
";

(*LineTransform code explained:
With two arguments: LineTransform[initialImage_, replacement_]. First we make sure that all line segments stand alone in the form of Line[{p1,p2}]. This is done by (N@initialImage/.line1:Line[{_,_,__}]:>Map[Line,Partition[First@line1,2,1]]).
The general plan is this: replace each line in initialImage by the specified motif. We rotate, scale, and translate the motif to each line's position to do this. alpha is the angle to rotate, and p1 is the translation needed.
With three arguments: LineTransform[initialImage_,replacement_,Line[{p3_,p4_}]]. Here we need to do an extra set of translations and rotations to put the motif in the right position, then move it onto each line in original image. This extra motion is combined into just one motion in the code.*)

LineTransform[initialImage_, replacement_] :=
   N[initialImage] /. line1:Line[{_, _, __}] :>
      Line /@ Partition[First[line1], 2, 1] /.
    Line[{p1_, p2_}] :> Block[{v, alpha, scale},
      v = p2 - p1; alpha = ArcTan @@ v; scale = N[Sqrt[v . v]];
       N[replacement] /. {line2_Line :>
          Map[{{Cos[alpha], -Sin[alpha]}, {Sin[alpha], Cos[alpha]}} . #1*scale +
             p1 & , line2, {2}], Point[pt_] :>
          Point[{{Cos[alpha], -Sin[alpha]}, {Sin[alpha], Cos[alpha]}} . pt*scale +
            p1]}];
LineTransform[initialImage_, replacement_, Line[{p3_, p4_}]] :=
   N[initialImage] /. line1:Line[{_, _, __}] :>
      Line /@ Partition[First[line1], 2, 1] /.
    Line[{p1_, p2_}] :> Block[{v, v2, alpha, scale},
      v = N[p2 - p1]; v2 = N[p4 - p3];
       alpha = N[ArcTan @@ v - ArcTan @@ v2];
       scale = N[Sqrt[v . v]/Sqrt[v2 . v2]]; N[replacement] /.
        {line2_Line :> Map[{{Cos[alpha], -Sin[alpha]}, {Sin[alpha], Cos[alpha]}} . (
                #1 - p3)*scale + p1 & , line2, {2}],
         Point[pt_] :> Point[{{Cos[alpha], -Sin[alpha]}, {Sin[alpha],
                Cos[alpha]}} . (pt - p3)*scale + p1]}];

Clear[Transform2DGraphics];



Transform2DGraphics::"usage" =
"Transform2DGraphics[f][graphics1] applies f to the coordinates of graphics primitves in graphics1.
f is a pure function having argument and output of the form {x,y}.
Some graphics primitves are not supported, such as three argument form of Circle, Disk, or Rectangle.
Example:
Transform2DGraphics[(#+{3,4})&][{PointSize[.05],Point[{0,0}]}]
";

Transform2DGraphics[f_Function][gra_] := 
  gra /. {poly_Polygon :> Map[f, poly, {2}], 
    li_Line :> Map[f, li, {2}], pt_Point :> Map[f, pt, {1}], 
    Text[str_, {a_, b_}, rest___] :> Text[str, f@{a, b}, rest], 
    Circle[{a_, b_}, radius_] :> Circle[f@{a, b}, radius], 
    Disk[{a_, b_}, radius_] :> Disk[f@{a, b}, radius], 
    Rectangle[{xm_, ym_}, {xM_, yM_}] :> 
     Polygon@(f /@ {{xm, ym}, {xM, ym}, {xM, yM}, {xm, yM}})};

Clear[Translate2DGraphics];

Translate2DGraphics::"usage" =
"Translate2DGraphics[{a1,a2}][graphics1] translates graphics1 by the vector {a1,a2}.
Example:
Graphics[Translate2DGraphics[{1,1}][{Line[{{0,0},{1,0},{0,1}}]}],Axes->True]
";

Translate2DGraphics[{0|0.,0|0.}][gra_]:=gra;
Translate2DGraphics[{a1_,a2_}][gra_]:=Transform2DGraphics[(#+{a1,a2})&][gra];

Clear[Rotate2DGraphics];

Rotate2DGraphics::"usage" =
"Rotate2DGraphics[{c1,c2},alpha][graphics1] rotates graphics1 by alpha around {c1,c2}.
Example:
Rotate2DGraphics[{0, 0}, alpha][Point[{1, 0}]] === Point[{Cos[alpha], Sin[alpha]}]
";

Rotate2DGraphics[{_,_},0|0.][gra_]:=gra;

Rotate2DGraphics[{0|0.,0|0.},alpha_][gra_]:=Transform2DGraphics[(({{Cos@alpha,-Sin@alpha},{Sin@alpha,Cos@alpha}}) . (#))&][gra];

Rotate2DGraphics[{a1_,a2_},alpha_][gra_]:=Transform2DGraphics[(({{Cos@alpha,-Sin@alpha},{Sin@alpha,Cos@alpha}}) . (#-{a1,a2})+{a1,a2})&][gra];

Clear[Reflect2DGraphics];
Reflect2DGraphics::"usage" =
"Reflect2DGraphics[][graphics1] reflects graphics1 by y-axes. Reflect2DGraphics[{a1,a2}][graphics1] reflects by the mirror line on vector {a1,a2}. Reflect2DGraphics[{a1,a2},{n1,n2}][graphics1] reflects by a line passing the point {n1,n2} and parallel to the vector {a1,a2}.
Example:
Reflect2DGraphics[{1, 0}][{Line[{{a1, a2}, {b1, b2}, {c1, c3}}]}] === {Line[{{a1, -a2}, {b1, -b2}, {c1, -c3}}]}
";

Reflect2DGraphics[][gra_]:=Reflect2DGraphics[{0,1}][gra];

Reflect2DGraphics[{0,0}][gra_]:=Transform2DGraphics[{-First@#,-Last@#}&][gra];

Reflect2DGraphics[{0,1}][gra_]:=Transform2DGraphics[{-First@#,Last@#}&][gra];

Reflect2DGraphics[{1,0}][gra_]:=Transform2DGraphics[{First@#,-Last@#}&][gra];

Reflect2DGraphics[{1,1}][gra_]:=Transform2DGraphics[{Last@#,First@#}&][gra];

Reflect2DGraphics[{a1_,a2_}][gra_]:=(Rotate2DGraphics[{0,0},ArcTan[a1,a2]])@Reflect2DGraphics[{1,0}]@(Rotate2DGraphics[{0,0},-ArcTan[a1,a2]][gra]);

Reflect2DGraphics[{a1_,a2_},{n1_,n2_}][gra_]:=(Translate2DGraphics[{n1,n2}])@(Reflect2DGraphics[{a1,a2}])@(Translate2DGraphics[-{n1,n2}][gra]);

Clear[GlideReflect2DGraphics];
GlideReflect2DGraphics::"usage" =
"GlideReflect2DGraphics[{a1,a2}][graphics1] glide reflects graphics1 by the vector {a1,a2}. GlideReflect2DGraphics[{a1,a2},{n1,n2}][graphics1] glide reflects by a line passing the point {n1,n2} and parallel to the vector {a1,a2}.
Example:
GlideReflect2DGraphics[{1, 1}, {1, 0}][{Line[{{0, 0}, {1, 0}, {0, 1}}]}] === {Line[{{2, 0}, {2, 1}, {3, 0}}]}
";

GlideReflect2DGraphics[{a1_,a2_}][gra_]:=Translate2DGraphics[{a1,a2}]@(Reflect2DGraphics[{a1,a2}][gra]);

GlideReflect2DGraphics[{a1_,a2_},{n1_,n2_}][gra_]:=Translate2DGraphics[{a1,a2}]@(Reflect2DGraphics[{a1,a2},{n1,n2}][gra]);

Clear[PolygonCounterOrientedQ];
(*Code based on Stan Wagon's book _Mathematica In Acton_ section 10.3.*)
PolygonCounterOrientedQ::"usage" =
"PolygonCounterOrientedQ[pointList] return True if the polygon formed by pointList is counterclockwise oriented, else it return False.
Example:
PolygonCounterOrientedQ[{{0,0},{1,0},{1,1},{0,1}}]
";

PolygonCounterOrientedQ[{{x1_,y1_},{x2_,y2_},{x3_,y3_}}]:=(1===Sign@Chop[Det@({{x1,y1,1},{x2,y2,1},{x3,y3,1}})/2.]);

PolygonCounterOrientedQ[pts:{{_,_}..}]:=Block[{n=Length@pts},PolygonCounterOrientedQ@Part[pts,(Position[#,First@Sort@#][[1,1]]&)@pts+{-1,0,1}/.{0->n,n+1->1}]];

Clear[LineIntersectionPoint,LineIntersectionPointCompiled];
LineIntersectionPoint::"usage" =
"LineIntersectionPoint[{A,B},{C,D}] return the intersection of lines AB and CD.
Example:
LineIntersectionPoint[{{0,0},{1,0}},{{0,0},{0,1}}]
";

If[$VersionNumber<3.,
	LineIntersectionPoint[{{a1_,a2_},{b1_,b2_}},{{c1_,c2_},{d1_,d2_}}]:={a2*b1*c1-a1*b2*c1-a2*b1*d1+a1*b2*d1-a1*c2*d1+b1*c2*d1+a1*c1*d2-b1*c1*d2,a2*b1*c2-a1*b2*c2-a2*c2*d1+b2*c2*d1-a2*b1*d2+a1*b2*d2+a2*c1*d2-b2*c1*d2}/(a2*c1-b2*c1-a1*c2+b1*c2-a2*d1+b2*d1+a1*d2-b1*d2),
	LineIntersectionPointCompiled=Compile[{a1,a2,b1,b2,c1,c2,d1,d2},{a2*b1*c1-a1*b2*c1-a2*b1*d1+a1*b2*d1-a1*c2*d1+b1*c2*d1+a1*c1*d2-b1*c1*d2,a2*b1*c2-a1*b2*c2-a2*c2*d1+b2*c2*d1-a2*b1*d2+a1*b2*d2+a2*c1*d2-b2*c1*d2}/(a2*c1-b2*c1-a1*c2+b1*c2-a2*d1+b2*d1+a1*d2-b1*d2)];
	LineIntersectionPoint[{{a1_,a2_},{b1_,b2_}},{{c1_,c2_},{d1_,d2_}}]:=LineIntersectionPointCompiled[a1,a2,b1,b2,c1,c2,d1,d2]];

Clear[LineIntersectionPoint2,LineIntersectionPoint2Compiled];
LineIntersectionPoint2::"usage" =
"LineIntersectionPoint2[{vec1,pt1},{vec2,pt2}] return the intersection of the following two lines: (1) line passing pt1 and parallel to vec1. (2) line passing pt2 and parallel to vec2. vec1 and vec2 must not be parallel and must not be zero vectors.
Example:
LineIntersectionPoint[{{1,0},{1,0}},{{0,1},{1,0}}]
";

If[$VersionNumber<3.,
	LineIntersectionPoint2[{{a1_,a2_},{b1_,b2_}},{{c1_,c2_},{d1_,d2_}}]:={(a2*b1*c1-a1*b2*c1-a1*c2*d1+a1*c1*d2),(a2*b1*c2-a1*b2*c2-a2*c2*d1+a2*c1*d2)}/(a2*c1-a1*c2),
	LineIntersectionPoint2Compiled=Compile[{a1,a2,b1,b2,c1,c2,d1,d2},{(a2*b1*c1-a1*b2*c1-a1*c2*d1+a1*c1*d2),(a2*b1*c2-a1*b2*c2-a2*c2*d1+a2*c1*d2)}/(a2*c1-a1*c2)];
	LineIntersectionPoint2[{{a1_,a2_},{b1_,b2_}},{{c1_,c2_},{d1_,d2_}}]:=LineIntersectionPoint2Compiled[a1,a2,b1,b2,c1,c2,d1,d2]];

Clear[PointOnLeftSideQ,PointOnLeftSideQCompiled];
PointOnLeftSideQ::"usage" =
"PointOnLeftSideQ[P,{A,B}] return True if point P is on the left side of the line AB. P, A, and B are coordinates of the form {real,real}.
Example:
PointOnLeftSideQ[{0,1},{{0,0},{1,0}}]
";

If[$VersionNumber<3.,
	PointOnLeftSideQ[{p1_,p2_},{{a1_,a2_},{b1_,b2_}}]:=N[(a2-b2)*p1+(-a1+b1)*p2+(-(a2*b1)+ a1*b2)]>0,
	PointOnLeftSideQCompiled=Compile[{p1,p2,a1,a2,b1,b2},N[(a2-b2)*p1+(-a1+b1)*p2+(-(a2*b1)+ a1*b2)]>0];
	PointOnLeftSideQ[{p1_,p2_},{{a1_,a2_},{b1_,b2_}}]:=PointOnLeftSideQCompiled[p1,p2,a1,a2,b1,b2]];

Clear[Wireframe2DGraphics];
Wireframe2DGraphics::"usage" =
"Wireframe2DGraphics[gra] changes Polygon and Rectangle in gra to the form Line[...], and Disk to Circle. Note: the three argument form of Rectangle is not converted. See also: ConvertCircleToLine.
Example:
Wireframe2DGraphics[{Rectangle[{0,0},{2,0}],Polygon[{{0,0},{0,-1},{1,-1}}]}]
";

Wireframe2DGraphics[gp_]:=gp/.{Disk->Circle,Polygon[{a_,rest1__}]:>Line[{a,rest1,a}],Rectangle[{xm_,ym_},{xM_,yM_}]:>Line@{{xm,ym},{xM,ym},{xM,yM},{xm,yM},{xm,ym}}};

Clear[ConvertCircleToLine];

ConvertCircleToLine::"usage" =
"ConvertCircleToLine[gra,maxLength] converts all Circle graphics primitives in gra to Line[...] form with all segments' length less than maxLength. See also: Wireframe2DGraphics.
Example:
ConvertCircleToLine[Circle[{0,0},1],.2]
";

ConvertCircleToLine[gp_,maxLength_]:=gp/.{Circle[{c1_,c2_},0,___]:>Point[{c1,c2}],Circle[{c1_,c2_},(r_)?NumberQ]:>With[{iStep=N[(2*Pi)/Ceiling[(2*r*Pi)/maxLength]]},Line@Table[r*{Cos[t],Sin[t]}+{c1,c2},{t,0,2*Pi,iStep}]],Circle[{c1_,c2_},{a_,b_}]:>With[{iStep=N[(2*Pi)/Ceiling[(2*Pi*Max[a,b])/maxLength]]},Line@Table[{a*Cos[t],b*Sin[t]}+{c1,c2},{t,0,2*Pi,iStep}]],Circle[{c1_,c2_},r_,{a_,b_}]:>With[{iStep=N[(b-a)/Ceiling[(2*r*Pi)/maxLength]]},Line@Table[r*{Cos[t],Sin[t]}+{c1,c2},{t,a,b,iStep}]]};

Clear[Cut2DGraphics];

Cut2DGraphics::"usage" =
"Cut2DGraphics[gra,{A,B}] cut graphics gra by the line AB.
Only Line and Polygon primitives are cut.
Others are untouched.
Example:
Cut2DGraphics[{Polygon[{{0,0},{1,0},{1,1},{0,1}}],Line[{{2,0},{2,1}}]},{{0,1/2},{1,1/2}}]
";

Cut2DGraphics[gp_,cuttingLine:{{_,_},{_,_}}]:=
gp/.{li_Line:>CutLine[li,cuttingLine],poly_Polygon:>CutPolygon[poly,cuttingLine]};

Clear[Slice2DGraphics];

Slice2DGraphics::"usage" =
"Slice2DGraphics[g,{A,B}] cut graphics g by the line AB.
The result has length two.
The first part contain primitives that lie on the left side of the cutting line and the last part are primitives on the right side.
Only Line and Polygon primitives are cut.
Point primitive are separated into sides.
Graphics directives or other primitives are untouched and duplicated in both part of the result.
Example:
Slice2DGraphics[{Polygon[{{0,0},{1,0},{1,1},{0,1}}],Line[{{2,0},{2,1}}]},{{0,1/2},{1,1/2}}]
";

(*0 is used as a temporary head.*)
Slice2DGraphics[gp_, cuttingLine : {{_, _}, {_, _}}] :=
({# /. (0[{a_List, _List}]) :> a, # /. (0[{_List, a_List}]) :> a} &)@
   ReplaceAll[ gp, {li_Line :> 0@SliceLine[li, cuttingLine],
     poly_Polygon :> 0@SlicePolygon[poly, cuttingLine],
     Point[pt_] :>
      0@If[PointOnLeftSideQ[pt,
         cuttingLine], {{Point@pt}, {}}, {{}, {Point@pt}}]}];

Clear[CookieStamp2DGraphics];

CookieStamp2DGraphics::"usage" =
"CookieStamp2DGraphics[gra,{pt1,pt2,...}]] cut and return graphics gra that lies inside the given polygon {pt1,pt2,...}. The polygon must be convex and not cross itself. The returned object is a List, not Graphics.
Example:
Graphics [CookieStamp2DGraphics[(Line/@Table[RandomReal[],{15},{15},{2}]),{{0,0},{1,1/2},{1/2,1}}] , Frame->True]
";

CookieStamp2DGraphics[gra_List,pts:{{_,_}..}]:=Fold[First@Slice2DGraphics[#1,#2]&,gra,((Partition[#,2,1]&)@(Append[#,First@#]&)@If[PolygonCounterOrientedQ@pts,pts,Reverse@pts])];

CookieStamp2DGraphics[gra_Graphics,pts:{{_,_}..}]:=CookieStamp2DGraphics[First@gra,pts];

Clear[CutPolygon];

CutPolygon::"usage" =
"CutPolygon[Polygon[{P1,P2,..}],{A,B}] cut polygon by the line AB.
The polygon can be nonconvex but must not cross itself.
The result polygons have the same orientation as the original.
See also: SlicePolygon.
Example:
CutPolygon[Polygon[{{0,0},{1,0},{1,1}}],{{1/2,0},{1/2,1}}]
";

If[$VersionNumber<3.,CutPolygon[Polygon[pts_],cuttingLine:{{_,_},{_,_}}]:=Module[{vertexes,cutPoints},
			vertexes=Map[(PointOnLeftSideQ[#,cuttingLine]@#)&,pts];

		If[SameQ@@(Head/@vertexes),Return@{Polygon@pts}];
			cutPoints=(If[SameQ@@(Head/@#),Null,LineIntersectionPoint[{#[[1,1]],#[[2,1]]},cuttingLine]]&)/@((Partition[#,2,1]&)@(Flatten[{#,First@#}]&)@vertexes);

			vertexes=(DeleteCases[#,Null,{1}]&)@(Flatten[#,1]&)@Transpose[{pts,cutPoints}];

		cutPoints=Sort@(DeleteCases[#,Null,{1}]&)@cutPoints;

		Polygon/@ReleaseHold@Fold[ReplaceAll[ReleaseHold@#1,{({a___,First@#2,b___,Last@#2,c___}:>Hold@Sequence[{First@#2,b,Last@#2},{Last@#2,c,a,First@#2}]),({a___,Last@#2,b___,First@#2,c___}:>Hold@Sequence[{First@#2,c,a,Last@#2},{Last@#2,b,First@#2}])}]&,{vertexes},Partition[cutPoints,2]]
	],
	CutPolygon[Polygon[pts_],cuttingLine:{{_,_},{_,_}}]:=Module[{vertexes,cutPoints},

			(*This is vertexes tagged with True False head to indicate side.*)
			vertexes=Map[(PointOnLeftSideQ[#,cuttingLine]@#)&,pts];

		(*for efficiency reasons, if all points lies on one side of cutting line, then the function is aborted immediately.*)
		If[SameQ@@(Head/@vertexes),Return@{Polygon@pts}];

		(*This is cut points intermingled with Null.*)
			cutPoints=(If[SameQ@@(Head/@#),Null,LineIntersectionPoint[{#[[1,1]],#[[2,1]]},cuttingLine]]&)/@((Partition[#,2,1]&)@(Flatten[{#,First@#}]&)@vertexes);

		(*This is vertexes and cut points in order of traversing the polygon sides.*)
			vertexes=(DeleteCases[#,Null,{1}]&)@(Flatten[#,1]&)@Transpose[{pts,cutPoints}];

		(*This is a list of cutPoints in order of knife action (may be reversed). (without Null)*)
		cutPoints=Sort@(DeleteCases[#,Null,{1}]&)@cutPoints;

		(*This is a list of cut polygons*)
		Polygon/@Fold[ReplaceAll[#1,{({a___,First@#2,b___,Last@#2,c___}:>Sequence[{First@#2,b,Last@#2},{Last@#2,c,a,First@#2}]),({a___,Last@#2,b___,First@#2,c___}:>Sequence[{First@#2,c,a,Last@#2},{Last@#2,b,First@#2}])}]&,{vertexes},Partition[cutPoints,2]]

		(*for versions prior to 3.0, the following line should be used instead*)
		(*Polygon/@ReleaseHold@Fold[ReplaceAll[ReleaseHold@#1,{({a___,First@#2,b___,Last@#2,c___}:>Hold@Sequence[{First@#2,b,Last@#2},{Last@#2,c,a,First@#2}]),({a___,Last@#2,b___,First@#2,c___}:>Hold@Sequence[{First@#2,c,a,Last@#2},{Last@#2,b,First@#2}])}]&,{vertexes},Partition[cutPoints,2]]*)
	]];

Clear[SlicePolygon];

SlicePolygon::"usage" =
"SlicePolygon[Polygon[{P1,P2,..}],{A,B}] cut polygon by the line AB.
The polygon can be nonconvex but must not cross itself.
The result list has length two.
The first part are polygons that lie on the left side of the cutting line.
The result polygons have the same orientation as the original.
Example:
SlicePolygon[Polygon[{{0,0},{1,0},{1,1}}],{{1/2,0},{1/2,1}}]
";

SlicePolygon[Polygon[pts_],cuttingLine:{{_,_},{_,_}}]:=
({Cases[#,True[b_]:>b,{1}],Cases[#,False[b_]:>b,{1}]}&)@
(((PointOnLeftSideQ[#[[1,2]],cuttingLine])@#&) /@ CutPolygon[Polygon@pts,cuttingLine]);

Clear[SlitLine];

SlitLine::"usage" =
"SlitLine[Line[{P1,P2,..}],{A,B}] cut the line by a segment with endpoints at A and B.
See also: CutLine, SliceLine.
Example:
SlitLine[Line[{{0,0},{1,1},{2,0}}],{{0,1/2},{1,1/2}}]
";

(*This module works on lines with two points only *)
SlitLine[Line[pts:{_List,_List}],cuttingLine:{{_,_},{_,_}}]:=ReplaceAll[(If[SameQ@@(PointOnLeftSideQ[#,cuttingLine]&/@pts),List@Line@pts,If[SameQ@@(PointOnLeftSideQ[#,pts]&/@cuttingLine),List@Line@pts,{Line@{First@pts,cutPoint},Line@{cutPoint,Last@pts}}]]),cutPoint->LineIntersectionPoint[pts,cuttingLine]];

(*This module works on lines that contains more than one point. If all points lies on one side of cutting line, then execution is terminated immediately with Return for efficiency reason.*)

If[$VersionNumber<3.,

	SlitLine[Line[pts_List/;(Length@pts>2)],cuttingLine:{{_,_},{_,_}}]:=Module[{temp},
		If[SameQ@@(temp=PointOnLeftSideQ[#,cuttingLine]&/@pts),Return@List@Line@pts];
		temp=MapThread[And[#1,#2]&,{(UnsameQ@@#&)/@Partition[temp,2,1],(PointOnLeftSideQ[First@cuttingLine,#]=!=PointOnLeftSideQ[Last@cuttingLine,#]&)/@Partition[pts,2,1]}];
		temp=Append[#,Null]&@MapThread[If[#2,0@LineIntersectionPoint[#1,cuttingLine],Null]&,{Partition[pts,2,1],temp}];
		temp=(DeleteCases[#,Null,{1}]&)@Flatten[Transpose@{pts,temp},1];
		(Line/@FixedPoint[(ReleaseHold@#/.{a__,0[pt_],b__}->Hold@Sequence[{a,pt},{pt,b}])&,List@temp,200(*safty limit*)])],

	SlitLine[Line[pts_List/;(Length@pts>2)],cuttingLine:{{_,_},{_,_}}]:=Module[{temp},
		If[SameQ@@(temp=PointOnLeftSideQ[#,cuttingLine]&/@pts),Return@List@Line@pts];
		temp=MapThread[And[#1,#2]&,{(UnsameQ@@#&)/@Partition[temp,2,1],(PointOnLeftSideQ[First@cuttingLine,#]=!=PointOnLeftSideQ[Last@cuttingLine,#]&)/@Partition[pts,2,1]}];
		temp=Append[#,Null]&@MapThread[If[#2,0@LineIntersectionPoint[#1,cuttingLine],Null]&,{Partition[pts,2,1],temp}];
		temp=(DeleteCases[#,Null,{1}]&)@Flatten[Transpose@{pts,temp},1];
		(Line/@ReplaceRepeated[List@temp,{a__,0[pt_],b__}:>Sequence[{a,pt},{pt,b}]])]];

(*(*testing code for SlitLine*)Show[
    Graphics[{Thickness[.02],
        SlitLine[
            Line@Table[Random[],{8},{2}],{{0,0},{1,1}/2}]/.a_Line:>{
              RGBColor[Random[],Random[],Random[]],a},{Red,Thickness[.001],
          Line@{{0,0},{1,1}/2}}}],AspectRatio->Automatic,PlotRange->{{0,1},{0,1}},Frame->True]*)

Clear[CutLine];

CutLine::"usage" =
"CutLine[Line[{P1,P2,..}],{A,B}] cut the line by the line AB.
See also: SlitLine, SliceLine.
Example:
CutLine[Line[{{1,1},{1,-1},{2,1},{2,-1},{3,1},{3,-1}}],{{0,0},{1,0}}]
";

(*This module works on lines with two points only *)
CutLine[Line[pts:{_List,_List}],cuttingLine:{{_,_},{_,_}}]:=
ReplaceAll[
(If[SameQ@@(PointOnLeftSideQ[#,cuttingLine]&/@pts),List@Line@pts,{Line@{First@pts,cutPoint},Line@{cutPoint,Last@pts}}]),
cutPoint->LineIntersectionPoint[pts,cuttingLine]];

(*This module works on lines that contains more than one point. If all points lies on one side of cutting line, then execution is terminated immediately with Return for efficiency reason.*)

If[$VersionNumber<3.,

	CutLine[Line[pts_List/;(Length@pts>2)],cuttingLine:{{_,_},{_,_}}]:=Module[{temp},
		If[SameQ@@(temp=Map[(PointOnLeftSideQ[#,cuttingLine])&,pts]),Return@List@Line@pts];
		temp=MapThread[(If[SameQ@@#2,Null,0@LineIntersectionPoint[#1,cuttingLine]]&),{Partition[pts,2,1],Partition[temp,2,1]}];
		Line/@FixedPoint[(ReleaseHold@#/.{a__,0[pt_],b__}->Hold@Sequence[{a,pt},{pt,b}])&,List@(DeleteCases[#,Null,{1}]&)@Flatten[Transpose@{pts,Append[temp,Null]},1],200(*safty limit*)]
	],

	CutLine[Line[pts_List/;(Length@pts>2)],cuttingLine:{{_,_},{_,_}}]:=Module[{temp},
		If[SameQ@@(temp=Map[(PointOnLeftSideQ[#,cuttingLine])&,pts]),Return@List@Line@pts];
		temp=MapThread[(If[SameQ@@#2,Null,0@LineIntersectionPoint[#1,cuttingLine]]&),{Partition[pts,2,1],Partition[temp,2,1]}];
		Line/@ReplaceRepeated[List@(DeleteCases[#,Null,{1}]&)@Flatten[Transpose@{pts,Append[temp,Null]},1],{a__,0[pt_],b__}:>Sequence[{a,pt},{pt,b}]]
	]];

Clear[SliceLine];

SliceLine::"usage" =
"SliceLine[Line[{P1,P2,..}],{A,B}] cut the line by the line AB.
The result list has length two.
The first part are lines that lie on the left side of the cutting line.
See also: CutLine, SlitLine.
Example:
SliceLine[Line[{{1,1},{1,-1},{2,1},{2,-1},{3,1},{3,-1}}],{{0,0},{1,0}}]
";

SliceLine[Line[pts_List],cuttingLine:{{_,_},{_,_}}]:=DeleteCases[#,Null,{2}]&@(If[PointOnLeftSideQ[First@pts,cuttingLine],#,Reverse@#]&)@Transpose@(Partition[Join[(CutLine[Line@pts,cuttingLine]),{Null,Null,Null}],2]);

Clear[WallpaperGroupData];
WallpaperGroupData::"usage" =
"WallpaperGroupData is a list of data about wallpaper groups, some of which are default values for WallpaperPlot. Each row contains: {ID Number,Conway notation,IUC notation,basis vectors,fundamental region,generators}.
Example:
TraditionalForm@TableForm@WallpaperGroupData";

WallpaperGroupData=Transpose@List[
			Range@17,
			{"o"|"o1","2222","**","xx","*2222","22*","22x","x*","2*22","442","*442","4*2","333","*333","3*3","632","*632"},
			{"p1","p2","pm"|"p1m","pg"|"p1g","pmm"|"p2mm","pmg"|"p2mg","pgg"|"p2gg","cm"|"c1m","cmm"|"c2mm","p4","p4m"|"p4mm","p4g"|"p4gm","p3","p3ml","p3lm","p6","p6m"|"p6mm"},
			{{{1,0},{1/10,7/10}},{{1,0},{1/10,7/10}},{{1,0},{0,-1 + GoldenRatio}},{{1,0},{0,-1 + GoldenRatio}},{{1,0},{0,-1 + GoldenRatio}},{{1,0},{0,-1 + GoldenRatio}},{{1,0},{0,-1 + GoldenRatio}},{{1/Sqrt[1 + 1/GoldenRatio^2],-(1/ (Sqrt[1 + 1/GoldenRatio^2]* GoldenRatio))},{1/Sqrt[1 + 1/GoldenRatio^2],1/(Sqrt[1 + 1/GoldenRatio^2]* GoldenRatio)}},{{1/Sqrt[1 + 1/GoldenRatio^2],-(1/ (Sqrt[1 + 1/GoldenRatio^2]* GoldenRatio))},{1/Sqrt[1 + 1/GoldenRatio^2],1/(Sqrt[1 + 1/GoldenRatio^2]* GoldenRatio)}},{{1,0},{0,1}},{{1,0},{0,1}},{{1,0},{0,1}},{{1,0},{1/2,Sqrt[3]/2}},{{1,0},{1/2,Sqrt[3]/2}},{{1,0},{1/2,Sqrt[3]/2}},{{1,0},{1/2,Sqrt[3]/2}},{{1,0},{1/2,Sqrt[3]/2}}},
			{Polygon[{{0,0},#1,#1+#2,#2}]&,Polygon[{{0,0},#1,#2}]&,Polygon[{{0,0},#1,#1+#2/2,#2/2}]&,Polygon[{{0,0},#1/2,#1/2+#2,#2}]&,Polygon[{{0,0},#1/2,#1/2+#2/2,#2/2}]&,Polygon[{{0,0},#1,#1+#2/4,#2/4}]&,Polygon[{{0,0},#1/2,#1/2+#2/2,#2/2}]&,Polygon[{{0,0},#1,#1+#2}]&,Polygon[{{0,0},#1,#1/2+#2/2}]&,Polygon[{{0,0},#1/2,#1/2+#2/2,#2/2}]&,Polygon[{{0,0},#1/2,#1/2+#2/2}]&,Polygon[{{0,0},#1/2,#2/2}]&,Polygon[{#1,(2*(#1+#2))/3,#2,(#1+#2)/3}]&,Polygon[{(#1+#2)/3,(2*(#1+#2))/3,#2}]&,Polygon[{(#1+#2)/3,#1,#2}]&,Polygon[{{0,0},#1/2,(#1+#2)/3,#2/2}]&,Polygon[{{0,0},#1/2,(#1+#2)/3}]&},
			{{},{Rotation[(#1+#2)/2,Pi]&},{Reflection[#1,#2/2]&},{GlideReflection[#1/2,#2/2]&},{Reflection[#1,#2/2]&,Reflection[#2,#1/2]&},{Reflection[#1,#2/4]&,Rotation[(#1+#2)/2,Pi]&},{GlideReflection[#1/2,#2/4]&,Rotation[(#1+#2)/2,Pi]&},{Reflection[#1+#2]&},{Reflection[-#1+#2,#1]&,Reflection[#1+#2]&},{Rotation[(#1+#2)/2,Pi/2]&},{Reflection[#1+#2]&,Rotation[(#1+#2)/2,Pi/2]&},{Reflection[-(#1/2)+#2/2,#1/2]&,Rotation[(#1+#2)/2,Pi/2]&},{Rotation[(#1+#2)/3,(2*Pi)/3]&},{Reflection[#1+#2]&,Rotation[(#1+#2)/3,(2*Pi)/3]&},{Rotation[(#1+#2)/3,(2*Pi)/3]&,Reflection[#2-#1,#1]&},{Rotation[{0,0},Pi/3]&},{Reflection[#1+#2]&,Rotation[{0,0},Pi/3]&}}
];

Clear[Translation];
Translation::"usage" =
"Translation[{x,y}] represents a translation by the vector {x,y}.";

Clear[GlideReflection];

GlideReflection::"usage" =
"GlideReflection[{x1,y1}] represents a glide reflection by the vector {x1,y1}.
GlideReflection[{x1,y1},{x2,y2}] is a glide reflection by the line parallel to {x1,y1} and passing {x2,y2}.
GlideReflection is used by WallpaperGroupData to represent group generators.";

Clear[Reflection];
Reflection::"usage" =
"Reflection[{x1,y1}] represents the reflection operation by the vector {x1,y1}. Reflection[{x1,y1},{x2,y2}] represents reflection by the line parallel to {x1,y1} and passing {x2,y2}. Reflection is used by WallpaperGroupData to represent group generators.";

Clear[Rotation];
Rotation::"usage" =
"Rotation[{x1,y1},alpha] represents a rotation of angle alpha around {x1,y1}. Rotation is used by WallpaperGroupData to represent group generators.";

Clear[getWallpaperGroupUnitMotif];

getWallpaperGroupUnitMotif[n_Integer,a_?VectorQ,b_?VectorQ,glideF_Function,mirrorF_Function,rotateF_Function]:=With[{z={0,0}},
		Part[Unevaluated[
			List[
		{{},{},{{},{},{},{}}},
		{{},{},{rotateF[2,#]&/@{z,a/2,a,b/2,b,a/2+b/2,a+b/2,b+a/2,a+b},{},{},{}}},
		{{},mirrorF/@{{z,a},{b/2,a+b/2},{b,a+b}},{{},{},{},{}}},
		{glideF/@{{z,a},{b/2,a+b/2},{b,a+b}},{},{{},{},{},{}}},
		{{},mirrorF/@{{z,a},{z,b},{a/2,a/2+b},{b/2,b/2+a},{a,a+b},{b,a+b}},{rotateF[2,#]&/@{z,a/2,a,b/2,b,a/2+b/2,a+b/2,b+a/2,a+b},{},{},{}}},
		{glideF/@{{z,b},{a/2,a/2+b},{a,a+b}},mirrorF/@{{b/4,b/4+a},{b*3/4,b*3/4+a}},{rotateF[2,#]&/@{z,a/2,a,b/2,b,a/2+b/2,a+b/2,b+a/2,a+b},{},{},{}}},
		{glideF/@{{a/4,a/4+b},{b/4,b/4+a},{a*3/4,a*3/4+b},{b*3/4,b*3/4+a}},{},{rotateF[2,#]&/@{z,a/2,a,b/2,b,a/2+b/2,a+b/2,b+a/2,a+b},{},{},{}}},
		{glideF/@{{a/2,a+b/2},{b/2,b+a/2}},mirrorF/@{{z,a+b}},{{},{},{},{}}},
		{glideF/@{{a/2,b/2},{a/2,a+b/2},{b/2,b+a/2},{a+b/2,b+a/2}},mirrorF/@{{z,a+b},{a,b}},{rotateF[2,#]&/@{z,a/2,a,b/2,b,a/2+b/2,a+b/2,b+a/2,a+b},{},{},{}}},
		{{},{},{rotateF[2,#]&/@{a/2,b/2,a+b/2,b+a/2},{},rotateF[4,#]&/@{z,a,b,a/2+b/2,a+b},{}}},
		{glideF/@{{a/2,b/2},{a/2,a+b/2},{b/2,b+a/2},{a+b/2,b+a/2}},mirrorF/@{{z,a+b},{z,a},{z,b},{a/2,a/2+b},{b/2,b/2+a},{a,b},{a,a+b},{b,a+b}},{rotateF[2,#]&/@{a/2,b/2,a+b/2,b+a/2},{},rotateF[4,#]&/@{z,a,b,a/2+b/2,a+b},{}}},
		{glideF/@{{z,a+b},{a/4,a/4+b},{b/4,b/4+a},{a*3/4,a*3/4+b},{b*3/4,b*3/4+a},{a,b}},mirrorF/@{{a/2,b/2},{a/2,a+b/2},{b/2,b+a/2},{a+b/2,b+a/2}},{rotateF[2,#]&/@{a/2,b/2,a+b/2,b+a/2},{},rotateF[4,#]&/@{z,a,b,a/2+b/2,a+b},{}}},
		{{},{},{{},rotateF[3,#]&/@{z,a,b,(a+b)*1/3,(a+b)*2/3,a+b},{},{}}},
		{glideF/@{{a/4,b/2},{b/4,a/2},{a/2,a+b/2},{b/2,b+a/2},{a*3/4,a/4+b},{b*3/4,b/4+a},{a+b/2,a*3/4+b},{b+a/2,b*3/4+a}},mirrorF/@{{z,a+b},{a/2,b},{b/2,a},{a,a/2+b},{b,b/2+a}},{{},rotateF[3,#]&/@{z,a,b,(a+b)*1/3,(a+b)*2/3,a+b},{},{}}},
		{glideF/@{{a/2,b/2},{a/2,a/2+b},{b/2,b/2+a},{a+b/2,b+a/2}},mirrorF/@{{z,a},{z,b},{a,b},{a,a+b},{b,a+b}},{{},rotateF[3,#]&/@{z,a,b,(a+b)*1/3,(a+b)*2/3,a+b},{},{}}},
		{{},{},{rotateF[2,#]&/@{a/2,b/2,a/2+b/2,a+b/2,b+a/2},rotateF[3,#]&/@{(a+b)*1/3,(a+b)*2/3},{},rotateF[6,#]&/@{z,a,b,a+b}}},
		{glideF/@{{a/4,b/2},{b/4,a/2},{a/2,b/2},{a/2,a/2+b},{b/2,b/2+a},{a/2,a+b/2},{b/2,b+a/2},{a*3/4,a/4+b},{b*3/4,b/4+a},{a+b/2,b+a/2},{a+b/2,a*3/4+b},{b+a/2,b*3/4+a}},mirrorF/@{{z,a+b},{z,a},{z,b},{a/2,b},{b/2,a},{a,b},{a,a/2+b},{b,b/2+a},{a,a+b},{b,a+b}},{rotateF[2,#]&/@{a/2,b/2,a/2+b/2,a+b/2,b+a/2},rotateF[3,#]&/@{(a+b)*1/3,(a+b)*2/3},{},rotateF[6,#]&/@{z,a,b,a+b}}}
			]
			],n]
	];

Clear[fundamentalMotifToUnitMotif];

fundamentalMotifToUnitMotif[n_Integer,motif_,v1_?VectorQ,v2_?VectorQ]:=(Fold[(If[Head@#2===Rotate2DGraphics,Table[Rotate2DGraphics[First@#2,Last@#2*i][#1],{i,0,2*Pi/Last@#2-1}],{#1,#2@#1}]&),motif,#]&)@(Through[(WallpaperGroupData[[n,6]]/.{GlideReflection->GlideReflect2DGraphics,Reflection->Reflect2DGraphics,Rotation->Rotate2DGraphics})[v1,v2]]);

Clear[MotifGraphics];
MotifGraphics::"usage" =
"MotifGraphics is an option for WallpaperPlot. MotifGraphics->gra specifies a fundamental motif to be used. gra must be a list of Point, Line, Polygon, and graphics directives, or a Graphics object. The default is MotifGraphics->(N@Line[{{2/12,1/12},{5/12,1/12},{5/12,2/12}}]).";

Clear[BasisVectors];
BasisVectors::"usage" =
"BasisVectors is an option for WallpaperPlot and RandomWallpaperPlot. BasisVectors->{{a1,a2},{b1,b2}} specifies the two basis vectors that generate the unit cell. Default is BasisVectors->Automatic. For exact default values, see WallpaperGroupData.";

(*PlotPoints*)

Clear[ShowSymmetryAndUnitCellGraphics];
ShowSymmetryAndUnitCellGraphics::"usage" =
"ShowSymmetryAndUnitCellGraphics is an option for WallpaperPlot. ShowSymmetryAndUnitCellGraphics->True specifies WallpaperPlot to draw symmetry elements and unit cell boundaries. The default is False.";

Clear[CenterGraphics];
CenterGraphics::"usage" =
"CenterGraphics is an option for WallpaperPlot. CenterGraphics->True specifies WallpaperPlot to center the graphics on {0,0}. The default is False.";

Clear[UnitCellGraphicsFunction];

UnitCellGraphicsFunction::"usage" =
"UnitCellGraphicsFunction is an option for WallpaperPlot.
UnitCellGraphicsFunction->f specifies that f is used to generate graphics primitives to represent unit cells.
f must be a pure function that accepts one argument, having the form {{x1,y1},{x2,y2},...}.
The list is the vertex set of some polygon.
f must return graphics primitives.
The default is UnitCellGraphicsFunction->({Hue[0.18],Polygon[#]}&).";

Clear[FundamentalRegionGraphicsFunction];

FundamentalRegionGraphicsFunction::"usage" =
"FundamentalRegionGraphicsFunction is an option for WallpaperPlot.
FundamentalRegionGraphicsFunction->f specifies that f is used to draw fundamental regions.
f must be a pure function that accepts one argument, having the form {{x1,y1},{x2,y2},...}.
The list is the vertex set of some polygon.
f must return graphics primitives.
The default is FundamentalRegionGraphicsFunction->({GrayLevel[0.4],Polygon[#]}&).";

Clear[GlideReflectionGraphicsFunction];

GlideReflectionGraphicsFunction::"usage" =
"GlideReflectionGraphicsFunction is an option for WallpaperPlot.
GlideReflectionGraphicsFunction->(f) specifies that f is used to draw glide reflection lines.
f must be a pure function that accepts two arguments in the form of {a1,a2} and {b1,b2}, representing the end points of a glide reflection line.
f must return a list of graphics primitives.
The default is GlideReflectionGraphicsFunction->({AbsoluteDashing[{5,5}],Line[{#1,#2}]}&).";

Clear[MirrorLineGraphicsFunction];

MirrorLineGraphicsFunction::"usage" =
"MirrorLineGraphicsFunction is an option for WallpaperPlot.
MirrorLineGraphicsFunction->f specifies that f is used to draw mirror lines.
f must be a pure function that accepts two arguments in the form of {a1,a2} and {b1,b2}, representing the end points of an mirror line.
f must return a list of graphics primitives.
The default is MirrorLineGraphicsFunction->(DoubleLineMotif[#1,#2,.05]&).
Example:
WallpaperPlot[17,MotifGraphics->{},ShowSymmetryAndUnitCellGraphics->True,GlideReflectionGraphicsFunction->({Blue,Line[{#1,#2}]}&),MirrorLineGraphicsFunction->({Red,Line[{#1,#2}]}&),PlotPoints->2]
";

Clear[RotationSymbolGraphicsFunction];

RotationSymbolGraphicsFunction::"usage" =
"RotationSymbolGraphicsFunction is an option for WallpaperPlot.
RotationSymbolGraphicsFunction->f specifies that f is used to draw rotation symbols.
f must be a pure function that accepts two arguments, an integer n and a coordinate {a,b}, representing the order and center of rotation respectively.
f must return a list of graphics primitives.
The default is RotationSymbolGraphicsFunction->({{GrayLevel[1],RotationSymbolMotif[#1,.1,0,#2]/.Line->Polygon},RotationSymbolMotif[#1,.1,0,#2]}&).
Example:
WallpaperPlot[16,MotifGraphics->{},ShowSymmetryAndUnitCellGraphics->True,RotationSymbolGraphicsFunction->({Text[#1,#2]}&),UnitCellGraphicsFunction->({}&)]
";

Clear[WallpaperPlot];

WallpaperPlot::"usage" =
"WallpaperPlot[n] plots a wallpaper having symmetry ID n. 1 <= n <= 17. WallpaperPlot[\"Conway notation\"] and WallpaperPlot[\"IUC notation\"] can also be used.
Example:
WallpaperPlot[17]
";

(* Options[WallpaperPlot]=Union[{MotifGraphics->(N@Line[{{2/12,1/12},{5/12,1/12},{5/12,2/12}}]),BasisVectors->Automatic,PlotPoints->{4,4},CenterGraphics->False,ShowSymmetryAndUnitCellGraphics->False,UnitCellGraphicsFunction->(Line@Append[#,First@#]&),FundamentalRegionGraphicsFunction->({}&),GlideReflectionGraphicsFunction->({AbsoluteDashing[{5,5}],Line[{#1,#2}]}&),MirrorLineGraphicsFunction->(DoubleLineMotif[#1,#2,.05]&)(\*({Red,Line[{#1,#2}]}&)*\),RotationSymbolGraphicsFunction->({{GrayLevel[1],RotationSymbolMotif[#1,.1,0,#2]/.Line->Polygon},RotationSymbolMotif[#1,.1,0,#2]}&)},(Options[Graphics]/.{(AspectRatio->_):>(AspectRatio->Automatic)})]; *)

Options[WallpaperPlot] =
  Union[{MotifGraphics -> (N@
       Line[{{2/12, 1/12}, {5/12, 1/12}, {5/12, 2/12}}]),
    BasisVectors -> Automatic, PlotPoints -> {4, 4},
    CenterGraphics -> False, ShowSymmetryAndUnitCellGraphics -> False,
     UnitCellGraphicsFunction -> (Line@Append[#, First@#] &),
    FundamentalRegionGraphicsFunction -> ({} &),
    GlideReflectionGraphicsFunction -> ({AbsoluteDashing[{5, 5}],
        Line[{#1, #2}]} &),
    MirrorLineGraphicsFunction -> (DoubleLineMotif[#1, #2, .05] &),
    RotationSymbolGraphicsFunction -> ({{GrayLevel[1],
         RotationSymbolMotif[#1, .1, 0, #2] /. Line -> Polygon},
        RotationSymbolMotif[#1, .1, 0, #2]} &)}, (Options[
      Graphics] /. {(AspectRatio -> _) :> (AspectRatio ->
         Automatic)})];

(*({Red,Line[{#1,#2}]}&)*)

WallpaperPlot[n_Integer, opts___Rule] :=
   Module[{lMotifGraphics, lBasisVectors, lPlotPoints,
     lShowSymmetryAndUnitCellGraphics, lCenterGraphics, UnitCellGF,
     FundamentalRegionGF, GlideReflectionGF, MirrorLineGF,
     RotationSymbolGF, unitCellGP, fundamentalRegionGP, wallpaperGP,
     symmetryElementsGP}, {lMotifGraphics, lBasisVectors, lPlotPoints,
       lShowSymmetryAndUnitCellGraphics, lCenterGraphics} =
      {MotifGraphics, BasisVectors, PlotPoints,
         ShowSymmetryAndUnitCellGraphics, CenterGraphics} /. {opts} /.
       Options[WallpaperPlot]; {UnitCellGF, FundamentalRegionGF,
       GlideReflectionGF, MirrorLineGF, RotationSymbolGF} =
      {UnitCellGraphicsFunction, FundamentalRegionGraphicsFunction,
         GlideReflectionGraphicsFunction, MirrorLineGraphicsFunction,
         RotationSymbolGraphicsFunction} /. {opts} /.
       Options[WallpaperPlot]; lPlotPoints =
      {First[Flatten[{lPlotPoints}]], Last[Flatten[{lPlotPoints}]]};
     If[lBasisVectors === Automatic, lBasisVectors =
       N[WallpaperGroupData[[n,4]]]]; If[Head[lMotifGraphics] ===
       Graphics, lMotifGraphics = First[lMotifGraphics]];
     {unitCellGP, fundamentalRegionGP, symmetryElementsGP} =
      If[lShowSymmetryAndUnitCellGraphics,
       {UnitCellGF[{{0, 0}, First[lBasisVectors],
          Plus @@ lBasisVectors, Last[lBasisVectors]}],
        FundamentalRegionGF[First[WallpaperGroupData[[n,5]] @@
           lBasisVectors]], getWallpaperGroupUnitMotif[n,
         First[lBasisVectors], Last[lBasisVectors],
         GlideReflectionGF[First[#1], Last[#1]] & ,
         MirrorLineGF[First[#1], Last[#1]] & , RotationSymbolGF]},
       {{}, {}, {}}]; wallpaperGP = fundamentalMotifToUnitMotif[n,
       N[lMotifGraphics], First[lBasisVectors], Last[lBasisVectors]];
     (Show[Graphics[#1, Sequence @@ FilterRules[
           Flatten[{opts, Options[WallpaperPlot]}], Options[Graphics]],
         AspectRatio -> Automatic]] & )[
      Map[Translate2DGraphics[#1][{unitCellGP, fundamentalRegionGP,
          wallpaperGP, symmetryElementsGP}] & ,
       With[{a = First[lBasisVectors], b = Last[lBasisVectors],
         c = If[lCenterGraphics, (1/2)*Plus @@ (lBasisVectors*
              lPlotPoints), 0]}, N[Table[a*i + b*j - c,
          {j, 0, Last[lPlotPoints] - 1}, {i, 0, First[lPlotPoints] -
            1}]]], {2}]]];

MapThread[(WallpaperPlot[#2|#3,opts___]:=WallpaperPlot[#1,opts])&,(Part[#,Range@3]&)@(Transpose[WallpaperGroupData])];

Clear[RandomWallpaperPlot];

RandomWallpaperPlot::"usage" =
"RandomWallpaperPlot[n,(opts)] plots a random wallpaper with nth symmetry group.
RandomWallpaperPlot[n,m,s,r,(opts)] base it on m-gons (m = 3, 4, 6) with side length s (0.05 <= s <= .5).
r controls the randomness of the wallpaper. (0 < r <=1)
Special options are: ColorFunction->(RandomColor) can be used to specify color.
PlotPoints->{4,4} can be specified to control the number of cells in the plot.
Example:
RandomWallpaperPlot[17,3,.1,1,PlotPoints->{3,3}]
";

Options[RandomWallpaperPlot] = Union[
{BasisVectors -> Automatic,
ColorFunction -> RandomColor,
PlotPoints -> {4, 4}},
Options[Graphics] /. {(AspectRatio -> _) :> AspectRatio -> Automatic}];

RandomWallpaperPlot[n_Integer, opts___] :=
   Module[{random1 = 1 - RandomReal[]/2}, RandomWallpaperPlot[n, 6,
      0.1, random1, opts]] /; 1 <= n <= 17 &&
     And @@ (#1 === Rule || #1 === RuleDelayed & ) /@
       Union[Head /@ {opts}];

RandomWallpaperPlot[n_Integer, m:3 | 4 | 6, sideLength_, r_,
    opts___] := Module[{basisVectors, colorFunction, plotPoints,
      allUnitCellsRange, fundamentalRegionRange, motifGP},
     basisVectors = N[BasisVectors /. {opts} /.
         Options[RandomWallpaperPlot]]; colorFunction =
       N[ColorFunction /. {opts} /. Options[RandomWallpaperPlot]];
      plotPoints = N[PlotPoints /. {opts} /.
         Options[RandomWallpaperPlot]]; plotPoints =
       {First[Flatten[{plotPoints}]], Last[Flatten[{plotPoints}]]};
      If[basisVectors === Automatic, basisVectors =
        N[WallpaperGroupData[[n,4]]]]; allUnitCellsRange =
       ({Min[#1], Max[#1]} & ) /@ Transpose[(Flatten[#1, 1] & )[
          ({{0, 0}, First[#1], First[#1] + Last[#1], Last[#1]} & ) /@
           Transpose[WallpaperGroupData][[4]]]];
      fundamentalRegionRange = ({Min[#1], Max[#1]} & ) /@
        Transpose @@ WallpaperGroupData[[n,5]] @@ basisVectors;
      motifGP = StainedGlassMotif[m, sideLength, r,
        fundamentalRegionRange]; motifGP =
       Fold[#1 /. p_Polygon :> First[SlicePolygon[p, #2]] & , motifGP,
        Partition[WallpaperGroupData[[n,5]] @@ basisVectors /.
          Polygon[{a_, rest__}] :> {a, rest, a}, 2, 1]];
      motifGP = motifGP /. p_Polygon :> {colorFunction[], p};
{
      Show[Graphics[{{GrayLevel[0.9], Polygon[
            ({{0, 0}, First[#1], First[#1] + Last[#1], Last[#1]} & )[
             basisVectors]]}, motifGP, First[WallpaperPlot[n,
            MotifGraphics -> {}, ShowSymmetryAndUnitCellGraphics ->
             True, BasisVectors -> basisVectors, PlotPoints -> {1, 1},
            DisplayFunction -> Identity]]}], AspectRatio -> Automatic,
        PlotRange -> {{-1, 1}, {-1, 1}}*0.1 + allUnitCellsRange,
        Sequence @@ FilterRules[Flatten[{opts, Options[
             RandomWallpaperPlot]}], Options[Graphics]], Frame -> True],
       WallpaperPlot[n, MotifGraphics -> motifGP, BasisVectors ->
         basisVectors, PlotPoints -> plotPoints,
        Sequence @@ FilterRules[Flatten[{opts}],
          Options[Graphics]]]
}
] /; 1 <= n <= 17 &&
     0.02 <= N[sideLength] <= 1 && 0 <= N[r] <= 1 &&
     And @@ (#1 === Rule || #1 === RuleDelayed & ) /@
       Union[Head /@ {opts}];

(*sideLength=Which[m===3,Sqrt[4/Sqrt[3]*area1],m===4,Sqrt@area1,m===6,Sqrt[2/(Sqrt[3]*3)*area1]];*)

End[];

EndPackage[];
