
BeginPackage["Midpoint Polygon Transformation`"]


(*Summer Research Program in Geometry Center,Minneapolis,1991*)
            (*  Funded by National Science Foundation*)
                      (*written by Arek Goetz*)

                  
MidPointPolygonMovie::usage=
"MidPointPolygonMovie[{{x1,y1},...,{xn,yn}},options] shows
frames of 'midpoint polygon' descendants of the polygon 
{{x1,y1},...,{xn,yn}}. In each corner of the frames there are
numbers which have the following meaning: 
upperleft - the frame number,    
upperright -the number of vertices of the polygon,    
lowerleft - iteration number,     
lowerright -absolute scaling factor.

MidPointPolygonMovie has the following options:

FirstFrameIteration -> 0,  
TotalFrameNumber  -> 2,
FractionPoint  -> 1/2,  
Blowing  -> Real,
IterationStep -> 1, 
Style ->Boundary,
ShowEllipse -> False,
VertexSize -> 5.";


FirstFrameIteration::usage="FirstFrameIteration,
is the index of the descandant shown in the first frame. 0
is the default value and means that the original
polygon is drawn in the first frame.";

TotalFrameNumber::usage="TotalFrameNumber is the number of
all frames.";
                    
FractionPoint::usage="FractionPoint defines the transformation
of a polygon by the formula:

T[{{x1,y1},...,{xn,yn}}] = {{fr x1+(1-fr)x2, fr
y1+(1-fr)y2},...,{fr xn+(1-fr)xn, fr yn+(1-fr)yn}}.
The default value is 1/2 which means that the vertices of
the transformed figure are the midpoints of the original
polygon.";
      
        
Blowing::usage="Blowing can be one of the words: Real,
Resized or ConstantPerimeter. They determine the scale of
the consecutive frames. 'Real' means that the scaling 
factor is always 1, 'Resized' causes the resizing of each
picture in such a way that the polygons tauch the edge of
frames at least at one point and are entirely showns in
the pictures. 'ConstantPerimeter' preserves the constant
perimeter in all frames.";

IterationStep::usage="IterationStep is a pure function
(1,#, #^2, 3# etc.) which determine how many more iterations should
be calculated between cosequetive frames. This function takes as an
argument the frame number. The function must give only natural numbers.
The default value is 1 which means that the consecutive frames will
show the consecutive iterations.";

Style::usage="Style determines the style of pictures.

Style -> Vertices - only vertices will be shown,
Style -> Boundary - the boundary of the polygons will be shown,
Style -> Filled   - the filled polygon will be shown,
Style -> VerticeOnBoundary - the vertices and the sides of polygons will be shown,
Style -> Spokes   - the vertices, the sides and the spokes will be shown.
A spoke is a line segment whose ends are the origin
and a vertex.";
                    
ShowEllipse::usage="ShowEllipse (True or False) calculates and
draws the 'midpoint ellipse'. The ellipse will appear in all frames.
The scaling factor is constant in this case and is also shown in all frames. This 
is to illustrate the main theorem.";

VertexSize::usage="VertexSize (a number 0...100) determines
the size of the vertices in all frames.";

SelfIntersecting::usage=
"The function SelfIntersecting[pol] returns True iff the
polygon pol is a selfintersecting one. A selfintrsecting polygon
is a polygon which has a pair of nonconsecutive sides that
have a common point."



RandomStellarPolygon::usage=
"RandomStellarPolygon[n,option] returns a random polygon whose 
descendants will eventually  become the affine images of
regular stellar n-gons. The centroids of these
polygons are at the origin. Basis[RandomStellarPolygon[n]]=
{0,0,...,ai,a(i+1),...,a(n-i-1),0,...,0}, where 1<=i<=n and the
coefficients ai,...,a(n-i-1) are chosen rangomly. The random
polygon is finally normalized so that
it touches the square {{1,1},{-1,1},{-1,-1},{-1,1}} at least
at one point and is entirely cotained in this square.
There is one option :  Precision ->16 which determines the precision
of calculations."
 
Basis::usage=
"The function Basis[pol] returns all complex coefficients
in the form of a list of the polygon pol written in the
following eigenvector basis:
 
v[j]={Cos[1j 2Pi/n]+I Sin[1j 2Pi/n],...,Cos[nj 2Pi/n]+I
Sin[nj 2Pi/n]},
                
where 1<=j<=n and n is the number of vertices of the polygon
pol. The vectors v[j], (1<=j<=n) are the all n regular, different
n-gons lying on the unit circle.";    

        
Perimeter::usage=
"Perimeter[pol] returns the perimeter of the polygon pol.";


MidPointPolygon::usage=
"MidPointPolygon[{{x1,y1},...,{xn,yn}}, fr] is a linear
transformation that returns the following polygon:
       
{{fr x1+(1-fr)x2, fr y1+(1-fr)y2},...,{fr xn+(1-fr)xn, fr
yn+(1-fr)yn}}.

If fr=1/2 the MidPointPolygon[pol] is a transformation which
returns the polygon whose vertices are the midpoints of the
sides of the polygon pol." 




Ellipse::usage=
"The function Ellipse[{{x1,y1},...,{xn,yn}}] returns the
complex coeficients a and b, and the lengths: c,d of the
semiexis in the form {a,b,c,d} such that the equation:
  
  z[û]= a Cos[û] + b Sin[û]                 (0<=û<2Pi)
 
 is the parametric equation of the 'midpoint ellipse'.";


ScalingFactor::usage=
"The function  ScalingFactor[{{x1,y1},...,{xn,yn}},fr]
 returns the scaling factor by which the consecutive,
 transformed polygons should be enlarged in order that
 the vertices lie on the 'midpoint ellipse'.";


Filled::usage = "Option for Style";
Boundary::usage = "Option for Style";
Spokes::usage = "Option for Style";
VerticesOnBoundary::usage = "Option for Style";
Vertices::usage = "Option for Style";

Resized::usage = "Option for Blowing";
Real::usage = "Option for Blowing";
ConstantPerimeter::usage = "Option for Blowing";




Begin["`private`"]
 messageSentence="`` is not a list of vertices of a polygon.";
 sentenceReal="`` is not a real number.";
                             (*max[obj] returns the radius of the minimal disk
                              centered  at the origin containing the object.*)
 max[obj_]:=Max[Abs[Flatten[obj]]];
 RealQ[s_]:=NumberQ[s] && Im[s]==0;

                (*checkFunction[f,n] checks if the image of the function f is 
                 a subset of natural numbers for the natural arguments to the
                  number n.*)  
 checkFunction[f_,n_]:=
  Block[{i,fun}, 
         fun=True;
         i=1;
         While[fun && i<=n,
                   If[IntegerQ[f[i]] && f[i]>0,i++,fun=False]
              ];
         Return[fun]
       ]; 
                 (*The function check[x] returns True iff x is a list of at least 
                   three points.*)
check[pol_]:=
 Block[{a,i},
       a=True;
       If[MatrixQ[pol],Do[
                          If[Not[VectorQ[pol[[i]]]
                                  && Length[pol[[i]]]==2
                                  && RealQ[pol[[i,1]]]
                                  && RealQ[pol[[i,2]]]
                               ],
                             a=False
                            ],
                          {i,1,Length[pol]}
                         ],  
                        a=False
          ];
        If[Length[pol]<3,a=False];
        Return[a]
      ];
      
     

Options[MidPointPolygonMovie]:= 
{
 FirstFrameIteration -> 0,
 TotalFrameNumber ->    2,
 FractionPoint          -> 1/2,
 Blowing                    -> Real,
 IterationStep          -> 1,
 Style                         -> Boundary,
 ShowEllipse            -> False,
 VertexSize              -> 5
};
  
MidPointPolygonMovie::notmat=messageSentence;
MidPointPolygonMovie::FFIoption="`` is not a nonnegative number.";
MidPointPolygonMovie::TFNoption="`` is not a positive number.";
MidPointPolygonMovie::FPoption="`` is not a real number between 0 and 1.";
MidPointPolygonMovie::Boption="`` is not one of the following words:
Real, Resized nor ConstantPerimeter.";
MidPointPolygonMovie::ISoption="`` is not an expression which can
be treated as a pure function like: 1, #, #^2, 3# etc. ,which takes only natural values.";
MidPointPolygonMovie::Soption="`` should be one of the words:
Vertices,Boundary,Filled,VerticesOnBoundary or Spokes.";
MidPointPolygonMovie::SEoption="`` should be True or False.";
MidPointPolygonMovie::VSoption="`` is not a nonegative real
number less than 100.";
MidPointPolygonMovie::polresize="Cannot resize the descendants of 
the polygon ``. 
The follwing pictures will be the real images.";




MidPointPolygonMovie[polygon_,options___]:=
Block[
 {
  firstFrameIteration,
  totalFrameNumber,
  fractionPoint,
  blowing,
  iterationStep,
  style,
  showEllipse,
  pointsize,
  number,
  k,
  nest,
  obj,
  pol,
  scale,
  end,
  i,
  qq,
  a,
  b,
  o,
  el
 },
     
     
      
 {firstFrameIteration,
  totalFrameNumber, 
  fractionPoint,
  blowing,
  iterationStep,
  style,
  showEllipse,
  pointsize
 } 
  
    = 
           
 {FirstFrameIteration,
  TotalFrameNumber,
  FractionPoint,
  Blowing,
  IterationStep&,
  Style,
  ShowEllipse,
  VertexSize
 } /. {options} /. Options[MidPointPolygonMovie]; 
       
end=False;     
Which[
      Not[check[polygon]],
           Message[MidPointPolygonMovie::notmat,polygon];
           end=True,
           
      Not[IntegerQ[firstFrameIteration] && firstFrameIteration >=0],
           Message[MidPointPolygonMovie::FFIoption,firstFrameIteration];
           end=True,

      Not[IntegerQ[totalFrameNumber] && totalFrameNumber >=1],
           Message[MidPointPolygonMovie::TFNoption,totalFrameNamber];
           end=True,
     
      Not[RealQ[fractionPoint] && Abs[fractionPoint]<1],
           Message[MidPointPolygonMovie::FPoption,FractionPoint];
           end=True,
          
      Not[MemberQ[{Real,Resized,ConstantPerimeter},blowing]],
          Message[MidPointPolygonMovie::Boption,blowing];
          end=True,
          
      Not[checkFunction[iterationStep,totalFrameNumber]],
          Message[MidPointPolygonMovie::ISoption,IterationStep];
          end=True,
          
      Not[MemberQ[{Vertices,Boundary,Filled,VerticesOnBoundary,Spokes},
                                   style
                                 ]
             ],
          Message[MidPointPolygonMovie::Soption,style];
          end=True,
          
      Not[MemberQ[{True,False},showEllipse]],
          Message[MidPointPolygonMovie::SEoption,showEllipse];
          end=True,
          
          
      Not[0 <= pointsize <= 100],
          Message[MidPointPolygonMovie::VSoption,pointsize];
          end=True
     ];
     
     
If[Not[end],
 
   number=Length[polygon];       (*number of vertices.*)
                                                        
   massCenter[obj_]:=Sum[ obj[[i]],{i,1,number}]/number;
                            (*The function    massCenter[{{x1,y1},...,{xn,yn}}]   
                            gives the position of  the center of mass of the given
                            points.*)
                            
   massCenteredObject[obj_]:=obj-
                             Table[massCenter[obj],{Length[obj]}];
             (*The function  massCenteredObject[{{x1,y1},...,{xn,yn}}] 
             changes the coordinates in such a way that the center of
             mass is at the origin  (translation).
                  Example:
                   massCenteredObject[{{1,1},{1,-1}}]={{0,1},{0,-1}} *)
  
  
  pol=massCenteredObject[polygon];
  blow=1;
  nest=firstFrameIteration;
  spokeFunction[qq_]:=Line[{{0,0},qq}];
 
  If[max[pol]==0 && (blowing == Resized || blowing == ConstantPerimeter),
       blowing=Real;
       Message[MidPointPolygonMovie::polresize,polygon]
    ]; 
   


    
    
  If [showEllipse, el=Ellipse[pol];
                   scale=ScalingFactor[pol,fractionPoint];
                   blowing="";
                   If[scale==Infinity,blowing=Real;
                                      showEllipse=False;
                                      scale=1; 
                      Message[MidPointPolygonMovie::polresize,polygon]
                      ];
                   blow=scale^firstFrameIteration//N;
                   pol=pol blow;           
                   a=el[[1]];
                   b=el[[2]];
                   ellipsePoints=Table[{Re[a Cos[Pi/40 o]+
                                           b Sin[Pi/40 o]
                                          ],
                                        Im[a Cos[Pi/40 o]+
                                           b Sin[Pi/40 o]
                                          ]
                                        },{o,0,80}
                                      ]//N
     ];
 
    midPointPolygon[obj_]:=
      Transpose[{fractionPoint Transpose[obj][[1]]+
                 (1-fractionPoint) RotateLeft[Transpose[obj][[1]]],
                 fractionPoint Transpose[obj][[2]]+
                 (1-fractionPoint) RotateLeft[Transpose[obj][[2]]]
               }];
            

 
   pol=Nest[midPointPolygon,pol,firstFrameIteration];

   If [ blowing==Resized, blow=blow /max[pol];
                          pol=pol/max[pol]
      ];

   If [blowing==ConstantPerimeter, blow=2/Perimeter[pol];
                                   pol=blow pol
      ];      
      
   If[showEllipse,Print["Each picture is resized by the the scaling factor:",
                                            scale 
                                       ]
      ];
      
                                      
   Do[ 
        If [blowing==Resized, blow=blow /max[pol];
                              pol=pol/max[pol]
           ];
        
        If [blowing==ConstantPerimeter, blow=blow 2/Perimeter[pol];
                     pol=pol 2/Perimeter[pol]
           ];
         
           
        Show[Graphics[{ 
                       Text[nest,{-1,-1},{-1.7,-1.7}],
                       Text[k,{-1,1},{-1.7,1.7}],                
                       Text[number,{1,1},{1.7,1.7}],
                       Text[blow//N,{1,-1},{1.7,-1.7}],
                         
                       If[showEllipse,Line[ellipsePoints],Text[k,{-1,1},{-1.7,1.7}]],
                                       
                       Which[blowing==Real,Text["Real",{0,1},{0,1.7}],
                             blowing==Resized,Text["Resized",{0,1},{0,1.7}],
                             showEllipse,     Text["Const. scaling factor",
                                                   {0,1},{0,1.5}
                                                                 ],
                             blowing==ConstantPerimeter,Text["Constant perimeter",
                                                             {0,1},{0,1.7}
                                                            ]
                            ],                                                              
                    
                    
                        
                    
                       Which[  style==Vertices,Flatten[{
                                                                                  PointSize[pointsize/200],
                                                                                  Map[Point,pol]
                                                                                }],
                                     style==Boundary,Flatten[{
                                                                                    Line [pol],
                                                                                    Line[{First[pol],Last[pol]}]
                                                                                  }],
                                     style==Filled,Polygon[pol],
                            
                                     style==VerticesOnBoundary,
                                                   Flatten[{
                                                                  PointSize[pointsize/200],
                                                                  Map[Point,pol],
                                                                  Line [pol],
                                                                  Line[{First[pol],Last[pol]}]
                                                                }],
                                      style==Spokes,Flatten[{
                                                                                 PointSize[pointsize/200],
                                                                                 Map[Point,pol],
                                                                                 Line [pol],
                                                                                 Line[{First[pol],Last[pol]}],
                                                                                 Map[spokeFunction,pol]
                                                                                }]                           
                                  ]                                     
                          },                                            
                      AspectRatio->Automatic,
                      PlotRange ->{{-1,1},{-1,1}},
                      Framed -> True                  
                     ]
           ]; 
             
      If [showEllipse, pol=scale^iterationStep[k] pol; 
                       blow=scale^iterationStep[k] blow//N
         ];    
       nest=nest+iterationStep[k];
       If[k<totalFrameNumber,
            pol=Nest[midPointPolygon,pol,iterationStep[k]]
         ],
      {k,totalFrameNumber}
     ] (* end of 'Do'*)
    ] (* end of 'If'*)
  ] (*end of 'Block' in MidPointPolygonMovie*)
  
 
      
(*The function intersecting[{a1,a2},{b1,b2},{c1,c2},{d1,d2}] returns
  True iff the line segments:
  {{a1,a2},{b1,b2}} and {{c1,c2},{d1,d2}} are disjoint.*)  
      
intersecting[{a1_,a2_},{b1_,b2_},{c1_,c2_},{d1_,d2_}]:=
 Block[
       {alfa,beta,mat,al,be},
       mat=(a1-b1)(d2-c2)-(d1-c1)(a2-b2);
       al=100;
       be=100;
       alfa=(d1-b1)(d2-c2)-(d1-c1)(d2-b2);
       beta=(a1-b1)(d2-b2)-(d1-b1)(a2-b2);
       If [mat != 0, al=alfa/mat; be=beta/mat];
       If [mat == 0 && alfa==0 && beta==0 &&
           (Sqrt[(a1-b1)^2+(a2-b2)^2]+
            Sqrt[(c1-d1)^2+(c2-d2)^2]
           )^2 >= (c2+d2-a2-b2)^2+
                  (c1+d1-a1-b1)^2 
                      ,al=0; be=0];
       Return[0 <= al <= 1 && 0 <= be <= 1]
      ];
      
      
norm[obj_]:=Sqrt[obj.obj];

      

SelfIntersecting::notmat=messageSentence;
SelfIntersecting[p_]:=
Block[
      {num,point,i,j,k,pol},
      If[check[p],
         num=Length[p];
         pol=Join[p,p];
         point=False;
         i=1;
         j=2;
         While[Not[point] && j<=num/2,
               While[Not[point] && i<=num, 
                     point=intersecting[pol[[i]],pol[[i+1]],
                                        pol[[i+j]],pol[[i+j+1]]
                                       ];
                     i++
                    ];i=1;j++
              ];
         Return[point],
         Message[SelfIntersecting::notmat,p]
        ]
     ];  
      
 


RandomStellarPolygon::range="`` cannot be the number of the vertices
of a RandomStellarPolygon.";
RandomStellarPolygon::Poption="`` is not a natural number.";
Options[RandomStellarPolygon]:={Precision -> 16};
RandomStellarPolygon[n_,option___ ]:=
 Block[{k1,k2,k,i,j,g,m,precision},
 precision=Precision /. {option}/. Options[RandomStellarPolygon];
      If[IntegerQ[precision] && precision > 0,      
        If[ IntegerQ[n] && n>=4,
             k1=Random[Integer,{1,(n-1.1)/2}];
             k2=Random[Integer,{1,(n-1.1)/2}];
             k=Min[k1,k2];
             pol1=RotateRight[Join[Table[I Random[Real,{-1,1},precision]+
                                           Random[Real,{-1,1},precision],
                                         {n-2k-1}
                                        ],Table[0,{1+2k}]
                                  ],k
                             ];
             pol1=N[pol1.Table[Table[Cos[2 i j Pi/n]+ I Sin[2 i j Pi/n],
                                     {i,n}
                                    ],
                               {j,n}
                              ],precision
                   ];
             pol1=Table[{Re[pol1[[g]]],Im[pol1[[g]]]},{g,n}];
           
             Return[pol1/max[pol1]],
      
             Message[RandomStellarPolygon::range,n]
             ],Message[RandomStellarPolygon::Poption,option]
           ]
              ];
     
     
           
Basis::notmat=messageSentence;       
Basis[pol_]:=
 Block[{num,w},
       If[check[pol],
           num=Length[pol];
           Return[N[1/num Table[Table[E^(-2 Pi/num i j I),{i,num}],{j,num}].
                                   Table[pol[[w,1]]+ I pol[[w,2]],{w,num}],
                              Precision[pol[[1,1]]]
                             ]
                        ],
           Message[Basis::notmat,pol]
         ]
      ];   
      






Perimeter::notmat=messageSentence;
Perimeter[obj_]:=
 Block[{i},
       If[check[obj],
             Return[Sum[norm[obj[[i+1]]-obj[[i]]],{i,1,Length[obj]-1}]+
                        norm[obj[[Length[obj]]]-obj[[1]]]//N
                         ], Message[Perimeter::notmat,obj]
           ]
              ];
                          


MidPointPolygon::notmat=messageSentence;
MidPointPolygon::type=sentenceReal;
        
MidPointPolygon[obj_,fractionPoint_]:=
      (
       If[check[obj] && RealQ[fractionPoint],
           Return[Transpose[{fractionPoint Transpose[obj][[1]]+
                            (1-fractionPoint) RotateLeft[Transpose[obj][[1]]],
                            fractionPoint Transpose[obj][[2]]+
                            (1-fractionPoint) RotateLeft[Transpose[obj][[2]]]
                           }]
                 ]
         ];
       If[Not[check[obj]],Message[MidPointPolygon::notmat,obj]];
       If[Not[RealQ[fractionPoint]],
           Message[MidPointPolygon::type,fractionPoint]
         ]
      );
      
      
 
 
Ellipse::notmat=messageSentence;
Ellipse::type=sentenceReal;
Ellipse[pol_]:=
 Block[{num,base,q},
      If[check[pol], 
          num=Length[pol];
          base=Basis[pol];
          q=1;While[q<num/2 && 
                    Chop[base[[q]]]==0 &&
                    Chop[base[[num-q]]]==0,
                    q++
                   ];
          Return[ Chop[{(base[[q]]+base[[num-q]]), 
                        I(base[[q]]-base[[num-q]]),
                        Abs[base[[q]]]+Abs[base[[num-q]]],
                        Abs[Abs[base[[q]]]-Abs[base[[num-q]]]]
                      }]
                ],
          Message[Ellipse::notmat,pol]
         ]
        ];



      
ScalingFactor::notmat=messageSentence;
ScalingFactor::type=sentenceReal;
ScalingFactor::nexist="ScalingFactor[``] does not exist since `
` is of the form {{A1,A2},{B1,B2},...,{A1,A2},{B1,B2}}."
ScalingFactor[pol_,frac_]:=
 Block[{num,base,q},
        If[check[pol] && RealQ[frac],
            num=Length[pol];
            base=Basis[pol];
            q=1;While[q<num/2 && 
                      Chop[base[[q]]]==0 &&
                      Chop[base[[num-q]]]==0,
                      q++
                     ];
             If[q<num/2,If[ frac !=.5, Return[1/Abs[frac+(1-frac) 
                                             E^(I q 2 Pi/num)]],
                                       Return[1/Cos[q Pi/num ]]
                          ],
                         Return[Infinity];
                         Message[ScalingFactor::nexist,pol,pol]
               ]
          ];
         If[Not[check[pol]],Message[ScalingFactor::notmat,pol]];
         If[Not[RealQ[frac]],Message[ScalingFactor::type,frac]]
      ];


 



End[]

EndPackage[]




