(* findangle finds the interior angle at v given three
vertices in order u,v,w, where we assume the domain is on the
left as we gor from u to v *)



  plotpoly[vertlist_]:=
          Show[Graphics[
             {
               GrayLevel[1],Thickness[.01], {Polygon[vertlist] } ,
               GrayLevel[0],Thickness[.01], {Line[vertlist] }
              },
             {PlotRange -> All, Axes -> None, AspectRatio -> Automatic}]]

  ctor[z_]:=(
     output={};
     Do[
     new={Re[z[[k]]],Im[z[[k]]]};
     AppendTo[output,new];
     ,{k,1,Length[z]}];
    output
  );


  plotpolyC[vert_]:=(
    vert1=ctor[vert];
    vert2=AppendTo[vert1,vert1[[1]] ];
    plotpoly[vert2]
   )


findangle[u_,v_,w_]:=(
    a=Abs[u-v];
    b=Abs[w-v];
    c=Abs[w-u];
    z=(w-v)/(v-u);
    If[(Im[z]==0)&&(Re[z] < 0), beta= Pi];
    If[(Im[z]==0)&&(Re[z]>0), beta=0];
    alpha=ArcCos[(c^2 - a^2 -b^2)/(-2 a b)];
    If[Im[z]>0, beta=alpha -Pi];
    If[Im[z]<0, beta= Pi - alpha];
    N[1+(beta/Pi)]
);

(* find angle list takes a listof complex numbers, treates
them as vertices of a planar polygona nd finds interior
angles *)

findanglelist[data_]:=(
   L =Length[data];
   output={findangle[data[[L]],data[[1]], data[[2]] ]};
   Do[
    AppendTo[output, findangle[data[[k-1]],data[[k]], data[[k+1]]]];
    ,{k,2,L-1}];
    AppendTo[output, findangle[data[[L-1]],data[[L]], data[[1]]]];
    output
);

F[t_,theta_,alpha_]:=Product[ Abs[1-Exp[ I (t-theta[[k]]) ] ]^(alpha[[k]]-1),
             {k,1,Length[theta]}];

mySCint3[a1_,b1_,theta_,alpha_]:=(
   lambda=.500000000000000000000000000000;
   num=30;
   num2=40;
   ylist=N[Table[a1+ lambda^(k-1) (b1-a1)/2 ,{k,1,num}],100];
    outsum1=0;
   Do[outsum1=outsum1+mySCint4[ylist[[k+1]],ylist[[k]],theta,alpha,num2],{k,1,num-1}];
   ylist=N[Table[b1- lambda^(k-1) (b1-a1)/2 ,{k,1,num}],100];
   Do[outsum1=outsum1+mySCint4[ylist[[k]],ylist[[k+1]],theta,alpha,num2],{k,1,num-1}];
   outsum1
);

mySCint4[a_,b_,theta_,alpha_,n_]:=(
   xlist=Table[a+(k-1)(b-a)/(2 n) ,{k,1,2n+1}];
   outsum=F[xlist[[1]],theta,alpha];
   outsum=outsum+F[xlist[[2n+1]],theta,alpha];
   Do[outsum=outsum+ 2 F[xlist[[2 k+1]],theta,alpha ],{k,1,n-1}];
   Do[outsum=outsum+ 4 F[xlist[[2 k]] ,theta,alpha],{k,1,n}];
   outsum=outsum (b-a)/(3(2n))
);

distance2edge[vertex1_,vertex2_,z_]:=(
     length=Abs[vertex2-vertex1];
     newz= length (z-vertex1)/(vertex2-vertex1);
     newx=Re[newz];
     newy=Im[newz];
     If[(0<=newx)&&(newx<=length)&&(newy>0),newy,Min[Abs[z-vertex1],Abs[z-vertex2]] ]
);
distance2edge2[vertex1_,vertex2_,z_]:=(
     length=Abs[vertex2-vertex1];
     newz= length (z-vertex1)/(vertex2-vertex1);
     newx=Re[newz];
     newy=Im[newz];
     If[(0<=newx)&&(newx<=length),Abs[newy],Min[Abs[z-vertex1],Abs[z-vertex2]] ]
);

(* distance to boundary find the distance from a complex number
to a polygon given by a listof vertices *)

distance2bdy[vert_,z_]:=(
     dist={};
     L=Length[vert];
     Do[
         AppendTo[dist, distance2edge[vert[[k]],vert[[k+1]],z] ];
     ,{k,1,L-1}];
     AppendTo[dist,  distance2edge[vert[[L]],vert[[1]],z] ];
     minimum=Min[dist];
     position=0;
     Do[ If[minimum==dist[[k]], position=k];
       ,{k,1,L}];
     {minimum,position}
);
distance2bdy2[vert_,z_]:=(
     dist={};
     L=Length[vert];
     Do[
         AppendTo[dist, distance2edge2[vert[[k]],vert[[k+1]],z] ];
     ,{k,1,L-1}];
     AppendTo[dist,  distance2edge2[vert[[L]],vert[[1]],z] ];
     minimum=Min[dist];
     position=0;
     Do[ If[minimum==dist[[k]], position=k];
       ,{k,1,L}];
     {minimum,position}
);
(* kakutani performs num  random walks starting at center  inside a polygon given by vertices
   vert until  each reaches within epsilon of the boundry. It returns a
   list of integers, which say how often the kth edge was hit.  *)

kakutani[vert_,center_,num_,epsilon_]:=(
        L =Length[vert];
        out=Table[0,{k,1,L}];
        debug={};
        pathlist={vert};
        Do[(
           debugz={center};
           z=center;
           dist=distance2bdy[vert,z];
           While[dist[[1]]>epsilon,(
              z=z+dist[[1]] Exp[I 2 Pi Random[] ];
              dist=distance2bdy[vert,z];
              AppendTo[debug,{z,dist}];
              AppendTo[debugz,z];
           )];
           out[[dist[[2]]]]=out[[dist[[2]]]]+1;
           AppendTo[pathlist,debugz];
         ),
        {k,1,num}];
        If[Min[out]>0,(
            numhits=out;
            out=N[2 Pi out/Apply[Plus,out]];
            thetapoints={0};
            Do[AppendTo[thetapoints,thetapoints[[k]]+out[[k]]],{k,1,Length[out]-1}];
            plotSCa[thetapoints,findanglelist[vert],vert[[2]]-vert[[1]]]
        ),out]     
);

(* plotSCa does not shade interior of polygon and alighns image in same
   orientation as original  *)
   
plotSCa[theta_,alpha_,firstedge_]:=(
    angle={ Arg[firstedge]  };
    Do[AppendTo[angle,angle[[k-1]]+ Pi (1- alpha[[k]]) ],{k,2,Length[alpha]}];
    sidelength={};
    newtheta=Append[theta,theta[[1]] ];
    Do[If[newtheta[[k]]<newtheta[[k-1]],newtheta[[k]]=
            newtheta[[k]]+ 2 Pi],{k,2,Length[newtheta]}];
    (*Plot[F[t,theta,alpha],{t,0,2 Pi}] *)
    Do[
       newside=mySCint3[newtheta[[k]],newtheta[[k+1]],theta,alpha];
       AppendTo[sidelength,newside]
    ,{k,1,Length[theta]}];
    vertices={{0,0}};
    verticesC={0};
    Do[(
       newpoint=vertices[[k,1]]+I vertices[[k,2]] +
                          sidelength[[k]]Exp[I angle[[k]] ];
       AppendTo[vertices,{Re[newpoint],Im[newpoint]}];
       AppendTo[verticesC,newpoint];
    ), {k,1,Length[theta]}];

   Show[Graphics[
             {
               GrayLevel[1],Thickness[.01], {Polygon[vertices] } ,
               GrayLevel[0],Thickness[.01], {Line[vertices] }
              },
             {PlotRange -> All, Axes -> None, AspectRatio -> Automatic}]]
);
mySCint[a_,b_,theta_,alpha_,n_]:=Sum[F[a+k(b-a)/n,theta,alpha],{k,1,n-1}]((b-a)/(n-1));

vert=Table[Exp[ 2 Pi I k/6],{k,1,6}];
center=0;
kakutani[vert,center,30,.01]
numhits
plotpolyC[vert]

vert={0,3,3+I, 4+I, 4+3I, 2+ 5I, -1 + 3I, I};
center=2 + 3I;
kakutani[vert,center,1000,.01]
numhits
plotpolyC[vert]