(* Content-type: application/mathematica *) (*** Wolfram Notebook File ***) (* http://www.wolfram.com/nb *) (* CreatedBy='Mathematica 6.0' *) (*CacheID: 234*) (* Internal cache information: NotebookFileLineBreakTest NotebookFileLineBreakTest NotebookDataPosition[ 145, 7] NotebookDataLength[ 10824, 341] NotebookOptionsPosition[ 10239, 317] NotebookOutlinePosition[ 10573, 332] CellTagsIndexPosition[ 10530, 329] WindowFrame->Normal ContainsDynamic->False*) (* Beginning of Notebook Content *) Notebook[{ Cell["\<\ (* 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]] 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)); \ \>", "Input", CellChangeTimes->{{3.46494495229196*^9, 3.464944991330719*^9}, { 3.464945023946058*^9, 3.464945054147516*^9}, {3.464945121703919*^9, 3.464945129169017*^9}, {3.464945507425739*^9, 3.464945541443512*^9}, { 3.464945829655068*^9, 3.464945832317475*^9}, 3.464945955995211*^9, { 3.464946234744772*^9, 3.464946249867779*^9}, {3.464946298990825*^9, 3.464946307119427*^9}}], Cell[BoxData[""], "Input", CellChangeTimes->{{3.464944939957135*^9, 3.464944942321242*^9}}], Cell[BoxData[""], "Input", CellChangeTimes->{{3.464944725319935*^9, 3.464944881297636*^9}, { 3.464945086903809*^9, 3.464945110664107*^9}, {3.464945144815256*^9, 3.464945178917764*^9}, {3.464945217818471*^9, 3.464945241245181*^9}, { 3.464945301340907*^9, 3.464945404150374*^9}, {3.464945463374961*^9, 3.464945467200841*^9}, {3.464945564848787*^9, 3.46494556653454*^9}, { 3.464945606938235*^9, 3.464945711547429*^9}, {3.464945854945627*^9, 3.464945913164083*^9}, {3.464946008641682*^9, 3.464946009100906*^9}}], Cell[CellGroupData[{ Cell[BoxData[{ RowBox[{ RowBox[{"vert", "=", RowBox[{"Table", "[", RowBox[{ RowBox[{"Exp", "[", " ", RowBox[{"2", " ", "Pi", " ", "I", " ", RowBox[{"k", "/", "6"}]}], "]"}], ",", RowBox[{"{", RowBox[{"k", ",", "1", ",", "6"}], "}"}]}], "]"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"center", "=", "0"}], ";"}], "\[IndentingNewLine]", RowBox[{"kakutani", "[", RowBox[{"vert", ",", "center", ",", "30", ",", ".01"}], "]"}], "\[IndentingNewLine]", "numhits", "\[IndentingNewLine]", RowBox[{"plotpolyC", "[", "vert", "]"}], "\[IndentingNewLine]"}], "Input", CellChangeTimes->{{3.464944725319935*^9, 3.464944881297636*^9}, { 3.464945086903809*^9, 3.464945110664107*^9}, {3.464945144815256*^9, 3.464945178917764*^9}, {3.464945723583538*^9, 3.464945725029148*^9}, { 3.464945989615157*^9, 3.464946025916149*^9}, {3.464946078026521*^9, 3.46494607943268*^9}, {3.464946146818266*^9, 3.464946206297543*^9}, { 3.464946256921222*^9, 3.464946258462572*^9}, {3.464946321164427*^9, 3.464946339733349*^9}}], Cell[BoxData[ RowBox[{"{", RowBox[{"7", ",", "4", ",", "4", ",", "5", ",", "2", ",", "8"}], "}"}]], "Output", CellChangeTimes->{{3.464945234135864*^9, 3.464945248449414*^9}, { 3.464945975089427*^9, 3.464946098055294*^9}, {3.464946182685716*^9, 3.464946211082052*^9}, {3.464946263631994*^9, 3.464946279661697*^9}, { 3.464946329575737*^9, 3.464946345111436*^9}}] }, Open ]], Cell[BoxData[""], "Input", CellChangeTimes->{{3.464946372405023*^9, 3.46494637719781*^9}}], Cell[CellGroupData[{ Cell[BoxData[{ RowBox[{ RowBox[{"vert", "=", RowBox[{"{", RowBox[{"0", ",", "3", ",", RowBox[{"3", "+", "I"}], ",", " ", RowBox[{"4", "+", "I"}], ",", " ", RowBox[{"4", "+", RowBox[{"3", "I"}]}], ",", " ", RowBox[{"2", "+", " ", RowBox[{"5", "I"}]}], ",", " ", RowBox[{ RowBox[{"-", "1"}], " ", "+", " ", RowBox[{"3", "I"}]}], ",", " ", "I"}], "}"}]}], ";"}], "\[IndentingNewLine]", RowBox[{ RowBox[{"center", "=", RowBox[{"2", " ", "+", " ", RowBox[{"3", "I"}]}]}], ";"}], "\[IndentingNewLine]", RowBox[{"kakutani", "[", RowBox[{"vert", ",", "center", ",", "1000", ",", ".01"}], "]"}], "\[IndentingNewLine]", "numhits", "\[IndentingNewLine]", RowBox[{"plotpolyC", "[", "vert", "]"}]}], "Input", CellChangeTimes->{{3.464946372405023*^9, 3.464946474071542*^9}}], Cell[BoxData[ RowBox[{"{", RowBox[{ "82", ",", "45", ",", "45", ",", "104", ",", "351", ",", "271", ",", "87", ",", "15"}], "}"}]], "Output", CellChangeTimes->{{3.464946435477893*^9, 3.464946486483311*^9}}] }, Open ]] }, WindowSize->{936, 750}, WindowMargins->{{150, Automatic}, {Automatic, 52}}, FrontEndVersion->"6.0 for Linux x86 (32-bit) (June 19, 2007)", StyleDefinitions->"Default.nb" ] (* End of Notebook Content *) (* Internal cache information *) (*CellTagsOutline CellTagsIndex->{} *) (*CellTagsIndex CellTagsIndex->{} *) (*NotebookFileOutline Notebook[{ Cell[568, 21, 6348, 207, 3470, "Input"], Cell[6919, 230, 92, 1, 32, "Input"], Cell[7014, 233, 523, 7, 32, "Input"], Cell[CellGroupData[{ Cell[7562, 244, 1076, 23, 143, "Input"], Cell[8641, 269, 375, 7, 30, "Output"] }, Open ]], Cell[9031, 279, 91, 1, 32, "Input"], Cell[CellGroupData[{ Cell[9147, 284, 856, 23, 121, "Input"], Cell[10006, 309, 217, 5, 30, "Output"] }, Open ]] } ] *) (* End of internal cache information *)