(*This notebook has Outer Biliards routines for any regular polygon (including the Digital Filter Map and the Dual Center Map). It will automatically  generate the First Famil S[k] tiles. The variable 'npoints'  below specifies the number of vertices. Change npoints to your value and Evaluate. The only round-off is in the actual Tangent Map routines but aside from that all calculations are exact, so the polygons,star polygons, star points and scales are exact, and any orbit will be exact if the initil point and 'corner' sequence is known. This FirstFanily.nb is designed to introduce the singularity set W which we call the 'web' using N = 14 as a sample.. For documentation see the PDFs.*)

npoints =14; dplaces =35; HalfN =If[OddQ[npoints], Floor[npoints/2], (npoints/2)-1];  (*dplaces is the number of decimal places*) (*HelfN =N/2-1 for N even and Floor[N/2] for N odd*)

Mac=Table[{Sin[2 Pi k/npoints],Cos[2 Pi k/npoints]},{k,1,npoints}];

RadiusFromHeight[height_,n_]:=height/Cos[Pi/n];
RadiusFromSide[side_,n_]:=side/(2*Sin[Pi/n]);
SideFromRadius[radius_,n_]:=2*radius*Sin[Pi/n];
HeightFromRadius[radius_,n_]:=radius*Cos[Pi/n];

If[OddQ[npoints],  r0 = RadiusFromHeight[1,npoints]; s0 = SideFromRadius[r0,npoints]; h0=1; rD = RadiusFromSide[s0,2*npoints];
hD = HeightFromRadius[rD,2*npoints]; Mc =r0*Mac; M0 = RotateRight[Mc]; Mc=Append[Mc, Mc[[1]]]];McTemp=Mc;

If[EvenQ[npoints],Mc=RotationTransform[-Pi/npoints][Mac]; r0 = RadiusFromHeight[1,npoints];rD = r0;Mc = r0*Mc;M0=RotateRight[Mc];
Mc=Append[Mc,Mc[[1]]]; s0=SideFromRadius[rD,npoints];h0 = 1.0;rM0=RadiusFromSide[s0,npoints/2]]; McTemp=Mc;

TangentMapRoutines=Module[{},Mc= N[Mc];Axy=Range[npoints];Mxy=Axy;Bxy=Axy; (*here it makes sense to use high precision round-off*)
Edges=Table[{Mc[[i]],Mc[[i+1]]},{i,1,npoints}];

Dxy=Table[{Edges[[i,2,1]]-Edges[[i,1,1]],Edges[[i,2,2]]-Edges[[i,1,2]]},{i,npoints}];(*Deltax and Deltay for each edge*)

Do[If[Abs[Dxy[[i,2]]]≤Abs[Dxy[[i,1]]],Axy[[i]]="True",Axy[[i]]="False"],{i,1,npoints}];(*we will swap x and y when slope>1*)

Do[If[Axy[[i]]=="True",Mxy[[i]]=Dxy[[i,2]]/Dxy[[i,1]];
If[Abs[Mxy[[i]]]<10^(-10),Mxy[[i]]=0],Mxy[[i]]=Dxy[[i,1]]/Dxy[[i,2]]];
If[Abs[Mxy[[i]]]<10^(-10),Mxy[[i]]=0],{i,1,npoints}]; (*the slopes of the edges*)

Do[If[Axy[[i]]=="True",Bxy[[i]]=Mc[[i,2]]-Mxy[[i]]*Mc[[i,1]],Bxy[[i]]=Mc[[i,1]]-Mxy[[i]]*Mc[[i,2]]],{i,1,npoints}];
(*intercepts*)
Do[f[i_,x_]:=N[Mxy[[i]],dplaces]*x+N[Bxy[[i]],dplaces],{i,1,npoints}];(*one function for each edge*)

Fc[p_]:=Module[{d=Range[npoints]},Do[d[[i]]=(p[[1]]-Mc[[i,1]])^2+(p[[2]]-Mc[[i,2]])^2,{i,1,npoints}];
Extract[Flatten[Position[d,Min[d]]],1]];(*Use nearest vertex for first corner*)

Cros[p_,i_]:=Dxy[[i,2]]*(Mc[[i,1]]-p[[1]])-Dxy[[i,1]]*(Mc[[i,2]]-p[[2]]);(*cross procuct to test for valid tangent*)

Nc[p_]:=Module[{},temp=cc;
While[Cros[p,cc]>10^(-14),cc=Mod[cc+1,npoints,1]];
If[temp==cc&&Unequal[cc,Fc[p]],Return[0],Return[cc]]];


Np[p_,i_]:=If[i>0,Flatten[{{-p+2*Mc[[i]]}}],p];(*next point routine:pnew=-p+2c*)

NextPoint[p_]:=N[Np[p,Nc[p]],dplaces];(*apply NP to corner to get next point in orbit*)

V[p_,k_]:=Module[{},cc=Fc[N[p,dplaces]];
NestList[NextPoint,p,k]] (*End TangentMapRoutines*)]; Mc = McTemp; (*use approximations only for Tangent map*)

Tau[p_]:= Last[V[p,1]];
 
(*The Tangent Map. For orbits use NestList with Tau but V above is faster*)
Tau2[p_]:= Last[V[p,2]];  (*Tau^2 is called the 'Return Map'*)

StepSequenceRoutines=Module[{},Ncp[p_]:=Module[{},For[i=cc,i≤npoints+1,i++,i=Mod[i,npoints,1];(*redo Nc to track skip sequences*)If[Cros[p,i]<-10^(-14),SD[[index]]=Mod[i-cc,npoints];
index=index+1;cc=i;Return[cc](*otherwise increment i*)]]];

Gp[p_]:=Np[p,Ncp[p]];(*For Vp,use old Np and Fc with Ncp*)(*redo Gc to track skip sequences*)

Vp[p_,k_]:=Module[{K},index=1;SD=Range[k+1];cc=Fc[p];
K=NestList[Gp,p,k+1];Lastq=Last[K];SD=Drop[SD,1]];
];

Cr[L_,xmin_,xmax_,ymin_,ymax_]:=Module[{smax=Length[L],SD=Range[Length[L]]},For[s=1,s≤smax,s++,SD[[s]]=s;
If[L[[s,1]]>xmax||L[[s,2]]>ymax,SD[[s]]=0];
If[L[[s,1]]<xmin||L[[s,2]]<ymin,SD[[s]]=0]];
Return[SD]];  (*used in Crop below *)

Crop[L_,xmin_,xmax_,ymin_,ymax_]:=Module[{},SD=Cr[L,xmin,xmax,ymin,ymax];
Ix=DeleteCases[SD,0];Return[L[[Ix]]];]; (*crops points to whatever rectangular region you specify*)

CropBox[A_]:=Crop[A,left,right,bottom,top]; (* crops to current values of left,right, bottom, top*)
CropB[v_]:=Select[v,#[[1]]>left&&#[[1]]<right&&#[[2]]<top&&#[[2]]>bottom&];


RotateCorner[w_,nc_,center_]:=Module[{q=Range[nc]},q[[1]]=w;For[s=2,s≤nc,s++,q[[s]]=RotationTransform[(1-s)*2*Pi/nc,center][w]];Return[q]];
  
(*Note the use of RotationTransform here which will not compile for versions before 6.0, but it is easy to rewrite this without using RotationTransform. Since most Mathematica code is optimized it is usually better to use their routines when possible.*)

Period[L_]:=Module[{n=Length[L],f=First[L]},For[i=2,i≤n,i++,If[Abs[L[[i,1]]-f[[1]]]<10^(-13)&&Abs[L[[i,2]]-f[[2]]]<10^(-13),Return[i-1]]]];
(* Period[K] will search the elements of K looking for a match for K[[1]]. You might occasionally need to adjust the match criteria.*)
Periodq[L_,q_]:=Module[{n=Length[L]},For[i=2,i≤n,i++,If[Abs[L[[i,1]]-q[[1]]]<10^(-13)&&Abs[L[[i,2]]-q[[2]]]<10^(-13),Return[i-1]]]];
(*this is used in FindPeriod below for extended searches, using fixed initial point q *)

FindPeriod[startpoint_,steps_,length_]:=Module[{q=startpoint,post=0},q1=q;Unprotect[K];
For[j=1,j≤steps,j++,K=V[q1,length];
post=Periodq[K,q];Print[j,"   ",post];q1=Last[K];
K= .;If[post>1,Abort[]]]];

Datumcrop[startpoint_,steps_,length_,name_]:=Module[{q=startpoint},Unprotect[K];w=Range[steps];Wcrop={};
For[j=1,j≤steps,j++,K=V[q,length];
w[[j]]=Crop[K,left,right,bottom,top];
Print[j,"  "Length[w[[j]]]];
Wcrop=Join[Wcrop,w[[j]]];q=Last[K];K= .];Print[q];
Print[Length[Wcrop]];Save[name<>ToString[j-1]<>"Wcrop",Wcrop]];

Projections = Module[{},Tx = Range[HalfN];Wc=Range[HalfN];
Tx=Table[Table[Mod[j,npoints,1],{j,1,npoints*k,k}],{k,1,HalfN}];
     Wc=Table[Table[Mc[[Tx[[k]][[j]]]],{j,1,npoints}],{k,1,HalfN}];(*remappings of vertices of Mom*)

IND[q_,maxn_]:=Module[{Ind=Range[maxn]},V[q,1];Ind[[1]]=Mod[cc,npoints,1];
St=Vp[q,maxn];
For[s=2,s≤maxn,s++,Ind[[s]]=Mod[Ind[[s-1]]+St[[s-1]],npoints,1]];
Return[Ind]];
PIM[q_,n_,j_]:=Module[{Pim=Range[n]},Pim[[1]]=q;
For[s=2,s≤n,s++,Pim[[s]]=Pim[[s-1]]+2*(Sum[(-1)^k*Wc[[j]][[Ind[[k]]]],{k,2*s-3,2*s-2}])];
Return[Pim]];

(*To view the Wc mappings in Grid Form: GraphicsGrid[{Table[Graphics[poly[Wc[[k]]]],{k,1,HalfN}]},Frame->All]
To view the projections in Grid Form: First set k = 250 or some other value;Px=Table[Graphics[{poly[Mom],Blue,Line[PIM[q1,k,j]]}],{j,1,HalfN}];
 Px[[1]] should be plotted with points so:
Px[[1]]=Graphics[{poly[Mom],AbsolutePointSize[1.0],Blue,Point[PIM[q1,k,1]]}];
Then partition the Px any way you desire, for example:
GraphicsGrid[{{Px[[1]],Px[[2]]},{Px[[3]],Px[[4]]}},Frame->All] *)
];
WebRoutines = Module[{},web[stepsize_,length_,depth__,offset_]:=Module[{jmax=2*(length-offset)/stepsize,m,z=offset/stepsize},Jxy= Range[npoints];For[s=1,s≤npoints,s++,Jxy[[s]]={};
If[s>1,m=s-1,m=npoints];h[x_]:=N[f[m,x]];
If[Axy[[m]]=="True",x=Mc[[s,1]];
If[Dxy[[m,1]]>0,dx=stepsize,dx=-stepsize];
x=x+z*dx;(*to get sign right*)For[j=1,j<jmax,j=j+2,x=x+dx;
Jxy[[s]]=Join[Jxy[[s]],V[{x,h[x]},depth]]],y=Mc[[s,2]];If[Dxy[[m,2]]>0,dy=stepsize,dy=-stepsize];
y=y+z*dy;
For[j=1,j<jmax,j=j=j+2,y=y+dy;
Jxy[[s]]=Join[Jxy[[s]],V[{h[y],y},depth]];
(*If*)]   (*For*)]  (*For*)]  (*web*)]; 

WebPoints[stepsize_,length_,depth__]:=Module[{},web[stepsize,length,depth,0];WP={};
For[i=1,i≤npoints,i++,WP=Union[WP,Jxy[[i]]]];
Return[WP]];

Gr[depth_]:=Show[Graphics[{AbsolutePointSize[1.0],poly[M0],Blue,Point[WebPoints[.01,6,depth]]},PlotRange->{{left,right},{bottom,top}}]];  (*Set box[center,length] to your region of interest and change Gr to match your needs*)
(*GraphicsGrid[{{Gr[0],Gr[1],Gr[2],Gr[3]},{Gr[4],Gr[5],Gr[6],Gr[7]}},Frame->All]*) (*four by 4 grid to different depths*)

(*End WebRoutines*)];
  
ExtendedWebRoutines = Module[{},
V[p_,k_]:=Module[{},cc=Fc[p];FixedPointList[NextPoint,p,k]];
Rt[p_]:=ReflectionTransform[{1,0}][p];
T2=RotationTransform[-4*Pi/7][T1];

FEdge[edgenumber_,length_,increment_,offset_]:=Module[{},m=Mod[edgenumber+1,npoints,1];
Ds[m_]:=If[m>1,Dxy[[m-1]]/EuclideanDistance[Mc[[m]],Mc[[m+1]]],Dxy[[npoints]]/EuclideanDistance[Mc[[npoints-1]],Mc[[npoints]]]];
Pstart[m_]:=Mc[[m]]+Ds[m]*offset;
Pstop[m_]:=Mc[[m]]+Ds[m]*(length+offset);
Pk[m_,k_]:=Pstart[m]+Ds[m]*increment*k;
H=Table[Pk[m,k],{k,1,Floor[(length-offset)/increment]}];
Return[H]];

W0[length_,increment_,offset_]:=Module[{},W1={};
For[i=1,i≤npoints,i++,W1=Union[W1,FEdge[i,length,increment,offset]]];
Return[W1]];

(*H0=Table[{x,Mc[[6]][[2]]},{x,Mc[[6]][[1]],GenStar[[1]],-.005}];*)
(*T1=Flatten[Table[V[H0[[k]],500],{k,1,Length[H0]}],1];N[Length[T1]/(500*1000)]*)
(*Show[Graphics[{poly[Mom],AbsolutePointSize[1.0],Blue,Point[T1]},PlotRange->{{left,right},{bottom,top}}]]*)
(*Graphics[{poly[Mom],AbsolutePointSize[1.0], Point[T1], Magenta, Point[H0]}]*)
(*H1=W0[2,.01,1];*)
(*Graphics[{poly[Mom],AbsolutePointSize[1.0],Point[H1],AbsolutePointSize[5.0]}]*)
(*H7=FEdge[3,.5,.001,0];Length[H7]*)
(*T1=Flatten[Table[V[H0[[k]],1000],{k,1,Length[H0]}],1]; N[Length[T1]/(4500*1000)]*)

(*End ExtendedWebRoutines*)];

box[center_,r_]:=Module[{},left=center[[1]]-r;right=center[[1]]+r;
top=center[[2]]+r;bottom=center[[2]]-r];
thinbox[center_,r_]:=Module[{},left=center[[1]]-r;right=center[[1]]+r;
top=center[[2]]+r/3;bottom=center[[2]]-r/3];
Mbox[Cx_]:=Module[{},left=Cx[[1]][[1]][[1]];
right=Cx[[1]][[2]][[1]];top=Cx[[1]][[2]][[2]];
bottom=Cx[[1]][[1]][[2]];];
(* box[{0,0},3] will set left,right,bottom,top to match a square of side 6 centered at {0,0} so you can say:
Show[Grahics[......, PlotRange->{{left,right},{bottom,top}}]], Mbox is for pasting (Control-C) a Mathematica-generated box to get a generic box with left,right,botom top instead of two corners*)

HeightFromSide[side_,n_]:=N[side/(2*Tan[Pi/n]),dplaces];
EvenScale[n_]:=(1-(Cos[Pi/n])^2)/Cos[Pi/n]^2;
OddScale[n_]:=(1-Cos[Pi/n])/Cos[Pi/n];
GenerationScale[n_]:=If[OddQ[n],OddScale[n],EvenScale[n]];
slope[p1_,p2_]:=(p1[[2]]-p2[[2]])/(p1[[1]]-p2[[1]]);
Deltaxy[d_,slope_]:={d*Cos[ArcTan[slope]],d*Sin[ArcTan[slope]]};
Distance[x_,y_]:=Sqrt[(x[[1]]-y[[1]])^2+(x[[2]]-y[[2]])^2];

App[W_]:=Append[W,W[[1]]];poly[M_]:=Line[App[M]];

DigitalFilter = Module[{}, w= N[2*Pi/npoints];a=2*Cos[w];
f[x_]:=Mod[x+1,2]-1;
Df[{x_,y_}]:={y,f[-x+a*y]}; 
Shear[{x_,y_}]:={x,x *Cos[w]+y *Sin[w]};
IShear[{x_,y_}]:={x,-x *Cot[w]+(y /Sin[w])};
DfToTr[A_]:=If[Length[A]==2,RotationTransform[Pi/2][IShear[A]],RotationTransform[Pi/2][IShear/@A]];
TrToDf[A_]:=If[Length[A]==2,Shear[RotationTransform[-Pi/2][A]],Shear/@RotationTransform[-Pi/2][A]];
S[{x_,y_}]:=Which[-x+a*y≥1,1,-x+a*y<-1,-1,True,0];];

(*Example - to find the first 5 terms in the S-sequence of {-.5,.9}
Orbit=NestList[Df,{-.5,.9},5];Ssequence=S/@Orbit;
Example to generate a web in Df space and then map it to Trapezoid space - which is also Tangent Map space. 
This will transform ellipses into circles*)
(*H0 = Table[{x,1},{x,-1,1,.01}]; DfWeb=Flatten[Table[NestList[Df,H0[[k]],200],{k,1,Length[H0]}],1];
TrWeb=DfToTr[DfWeb];Graphics[Point[TrWeb]]*)

(*Dual Center Routines are below.  Dc extends complex z by +/- one unit and rotates by w. So N and -N will be constructed
from the interva [-1,1] as shown below. This symmetry is inherent in the outer-billiards map.*) 

DualCenter = Module[{}, w= N[2*Pi/npoints]; Dc[z_]:=Exp[-I*w]*(z-Sign[Im[z]]);

H0=Table[x,{x,-1,1,.005}]; (*iterate this interval under Dc to generate N and -N*)
Web=Flatten[Table[NestList[Dc,H0[[k]],npoints],{k,1,Length[H0]}],1];Webb={Re[#],Im[#]}&/@Web;  
Graphics[{AbsolutePointSize[1.0],Blue,Point[Webb]},Axes->True];
(* The parameters of N will be rN = RadiusFromSide[1,npoints];hN=HeightFromRadius[rN,npoints];cN={1/2,hN}; The trnsformation is ;
Tz[x_]:=(TranslationTransform[-cM][x])/hM;*)
];

FirstFamilyRoutines =Module[{},

(*N-odd*)If [OddQ[npoints],

h0=1;r0= RadiusFromHeight[h0,npoints];star =Table[{-Tan[k*Pi/npoints],-1}, {k,1,HalfN}];  scale[k_]:=Tan[Pi/npoints]/Tan[k*Pi/npoints];
GenScale = scale[HalfN];GenStar = star[[HalfN]];CFR[r_]:=(1-r)*GenStar; 
OuterDual[x_]:=RotationTransform[-Pi/npoints][x*r0];cS[k_]:=OuterDual[star[[k]]];
hD = 2 + GenScale;rD = RadiusFromHeight[hD,2*npoints];
hS[k_]:=hD*GenScale/scale[k]; rS[k_]:= RadiusFromHeight[hS[k],(2*npoints)];
S[k_]:=RotateCorner[cS[k]+{rS[k],0},npoints*2,cS[k]];S[0] =M0;MFF=Table[S[k],{k,0,HalfN}];

hDS[k_]:=hD*Tan[Pi/(2*npoints)]*Tan[k*Pi/(2*npoints)]; (*this is hD*GenScale[2N}/scale[k] of 2N*)
rDS[k_]:=If[OddQ[k],RadiusFromHeight[hDS[k],npoints],RadiusFromHeight[hDS[k],2*npoints]];
cDS[k_]:=If[OddQ[k],CFR[hDS[k]],CFR[hDS[k]]];
DS[k_]:=If[OddQ[k],RotateCorner[cDS[k]+{0,rDS[k]},npoints,cDS[k]],RotateCorner[cDS[k]+{rDS[k],0},2*npoints,cDS[k]]];
DFF=Table[DS[k],{k,1,npoints-4}]; 
FirstFamily=Join[MFF,DFF]];

(*N twice-odd*) If[OddQ[npoints/2],

h0 = 1;hD=1;
star=Table[{-Tan[k*Pi/npoints],-1}, {k,1,HalfN}];scale[k_]:=Tan[Pi/npoints]/Tan[k*Pi/npoints];
OuterDual[x_]:=RotationTransform[-Pi/npoints][x*rD];cS[k_]:=OuterDual[star[[k]]];
GenStar=star[[HalfN]];GenScale= scale[HalfN];
 rS[k_]:=If[OddQ[k],RadiusFromHeight[cS[k][[2]]+1,HalfN+1],rD*GenScale/scale[k]];
S[k_]:=If[OddQ[k],RotateCorner[cS[k]+{0,rS[k]},HalfN+1,cS[k]],RotateCorner[cS[k]+{rS[k],0},npoints,cS[k]]];
hS[k_]:=If[OddQ[k],HeightFromRadius[rS[k],npoints/2],HeightFromRadius[rS[k],npoints]];

S[0]= Mc;FFR=Table[S[k],{k,0,HalfN-1}];
cM=cS[HalfN-1];cDS[k_]:=ReflectionTransform[{1,0},cM][cS[k]];
DS[k_]:=ReflectionTransform[{1,0},cM][S[k]]; hDS[k_]:=hS[k];rDS[k_]:=rS[k];
FirstFamily=Join[FFR,ReflectionTransform[{1,0},cM]/@FFR]
];

(*N twice-even *) If[EvenQ[npoints/2],

h0 = 1;hD=1;star =Table[{-Tan[k*Pi/npoints],-1}, {k,1,HalfN}];scale[k_]:=Tan[Pi/npoints]/Tan[k*Pi/npoints];
GenStar=star[[HalfN]];GenScale= scale[HalfN];
OuterDual[x_]:=RotationTransform[-Pi/npoints][x*rD];cS[k_]:=OuterDual[star[[k]]];
If [npoints==4,cM0 = {-1,0},cM0 = cS[HalfN-1]];
cDS[k_]:=ReflectionTransform[{1,0},cM0][cS[k]];
rS[k_]:=rD*GenScale/scale[k];
hS[k_]:=HeightFromRadius[rS[k],npoints]; hDS[k_]:=hS[k];rDS[k_]:=rS[k];
Temp[k_]:=RotateCorner[cS[k]+{rS[k],0},npoints,cS[k]];
S[k_]:=RotationTransform[-Pi/npoints,cS[k]][Temp[k]];S[0]= Mc;D= S[HalfN]; 
DS[k_]:=ReflectionTransform[{1,0},cM0][S[k]];

FFR=Table[S[k],{k,0,HalfN-1}];FFD=Table[DS[k],{k,0,HalfN-2}];
FirstFamily=Join[FFR,FFD];
]];

V[p_, k_] := Module[{}, cc = Fc[N[p, dplaces]];
   NestList[NextPoint, p, k]]; (*Duplicate definition*)


(* Matrix Formats for star points and scales *)
(*MatrixForm[Table[-star[[k]][[1]], {k,1,HalfN}]]*)
(*MatrixForm[Table[scale[k], {k,1,HalfN}]]*)

(*Graphic to show  the forward and trailing extended edges of N 
Pgon[p_,q_]:=Table[{Cos[2 Pi k q/p],Sin[2 Pi k q/p]},{k,1,p}];
d1= EuclideanDistance[{0,0},GenStar];
P0 =RotationTransform[-Pi/2*npoints][Pgon[npoints,HalfN]*d1];
P1=RotationTransform[-Pi/npoints][Pgon[npoints,HalfN]*d1];
P2= {RotationTransform[-Pi/npoints][Pgon[npoints/2,HalfN/2]*d1],RotationTransform[Pi/npoints][Pgon[npoints/2,HalfN/2]*d1]};
Graphics[{Cyan,poly/@FirstFamily,Blue,If[OddQ[npoints],poly[P0]],If[EvenQ[npoints/2],poly[P1]],If[OddQ[npoints/2], poly/@P2] }]  *)


Graphics[{poly /@ FirstFamily,  Text["N = ", {-.2, .2} ], Text[npoints, {.2, .2}]}, Axes -> True, Ticks -> False]

(*Step 1: Standard definitions of S[1] and S[2] as simply S1 and S2:*)

S1 = S[1]; S2 = S[2]; cS1 = cS[1]; cS2 = cS[2]; hS1 = hS[1]; hS2 = hS[2];

(*Step 2: Midpoints and Star points*)

 MidS1={cS1[[1]],-1}; MidS2={cS2[[1]],-1};

StarS1=Table[MidS1+{hS1*Tan[k*Pi/(2*npoints)],0},{k,1,4}];
StarS2=Table[MidS2+{hS2*Tan[k*Pi/(2*npoints)],0},{k,1,npoints-1}];


StarS1R=ReflectionTransform[{1,0},cS1][StarS1];
StarS2R=ReflectionTransform[{1,0},cS2][StarS2];

V[p_, k_] := Module[{}, cc = Fc[N[p, dplaces]];
   NestList[NextPoint, p, k]]; (*Duplicate definition*)

(* To find the first 500 points in the orbit of a point p set K = V[p,500]and then Period[K] will check for possible period*)

(* Below are instructions to find the first 2 invariant regions for N = 14 . Since N is twice-odd the first regions is bounded  by [Floor[N/4] ]= 3 but the matching star point is star[4], so this initial region goes from star[1] to star[4] and we will color it blue: The second invariant region will go from star[5] to star[4] and we will color it Magenta. As always there are multiple ways to generate these invariant regions, but the normal methods are to iterate 'test' points with relatively dense orbits or iterate a collection of points on the horizontal cross section of the region to generate what we call a local web. Here we do two local webs for the first two invariant regions of N = 14.  *)

s1 = star[[1]]; s4 = star[[4]]; s5 = star[[5]];
     H0=Quiet[W0[-s4[[1]]-s0/2, .01,-s1[[1]]-s0/2]];     H1=Quiet[W0[-s5[[1]]-s0/2, .01,-s4[[1]]-s0/2]];  (*about 1400 points in each*)
    (*The Quiet command will keep Mathematica from reminding us that these are exact points and iteration may downgrade them*)
     
               (*For H0 and H1 these will be points spaced by .01 on each edge of N . Now iterate each point using the NestList command in Mathematica in its shortened 'V' form:    *)

BlueWeb = Flatten[Table[V[H0[[k]], 1], {k, 1, Length[H0]}], 50];   MagentaWeb = Flatten[Table[V[H1[[k]], 1], {k, 1, Length[H1]}], 50];

Graphics[{Black,poly[Mc], Cyan, poly/@FirstFamily,Blue,AbsolutePointSize[1.00],Point[BlueWeb] ,Magenta,Point[MagentaWeb]}]