(*this module writes a list of complex numbers into a list of {real,imaginary} pairs*)
cz[sool_]:=Module[{nofz,tofz},nofz=Dimensions[sool][[1]];
tofz=Table[{Re[sool[[ll]]],Im[sool[[ll]]]},{ll,1,nofz}];Return[tofz]];

(*caculates the zeros, poles and sigular points of a function ff*)

zps[LL1_,LL2_]:=Module[{ff,lam,PP,nder,dim},
dim=2;
ff=Sum[1/((z+2*(dim-Cos[2*Pi*n1/LL1]-Cos[2*Pi*n2/LL2]))),{n1,0,LL1-1},{n2,0,LL2-1}]/(LL1*LL2);
lam=Factor[ff];
PP=Denominator[lam];QQ=Numerator[lam];poles=z/.NSolve[QQ==0,z,WorkingPrecision->30];zeros=z/.NSolve[PP==0,z,WorkingPrecision->30];
nder=Expand[D[PP,z]*QQ-D[QQ,z]*PP];singp=z/.NSolve[nder==0,z,WorkingPrecision->30]]

zps[4,4];

temp=Show[ListPlot[cz[zeros],PlotStyle->{RGBColor[0,1,0],PointSize[0.02]}],
ListPlot[cz[poles],PlotStyle->{RGBColor[1,0,0],PointSize[0.02]}],
ListPlot[cz[singp],PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}]]

zps[8,8];

temp=Show[ListPlot[cz[zeros],PlotStyle->{RGBColor[0,1,0],PointSize[0.02]}],
ListPlot[cz[poles],PlotStyle->{RGBColor[1,0,0],PointSize[0.02]}],
ListPlot[cz[singp],PlotStyle->{RGBColor[0,0,1],PointSize[0.02]}]]

xmin=-8;xmax=0;
ymin=-0.8;ymax=0.8;
dellx=2;delly=0.2;

xla=Join[Table[{x,x,{0.0,0.}},{x,xmin,xmax,dellx}],
Table[{x," ",{0.025,0.},{GrayLevel[0.],AbsoluteThickness[1.2]}},{x,xmin-dellx,xmax+dellx,dellx/2}]];

xlatop=Join[Table[{x," ",{0.0,0.}},{x,xmin,xmax,dellx}],
Table[{x," ",{0.025,0.},{GrayLevel[0.],AbsoluteThickness[1.2]}},{x,xmin-dellx,xmax+dellx,dellx/2}]];


yla=Join[Table[{x,x,{0.0,0.}},{x,ymin,ymax,delly}],
Table[{x," ",{0.025,0.},{GrayLevel[0.],AbsoluteThickness[1.2]}},{x,ymin-delly,ymax+delly,delly/2}]];

ylr=Join[Table[{x," ",{0.0,0.}},{x,ymin,ymax,delly}],
Table[{x," ",{0.025,0.},{GrayLevel[0.],AbsoluteThickness[1.2]}},{x,ymin-delly,ymax+delly,delly/2}]];


final=Show[temp,Frame->True,AspectRatio->1,Axes->False,PlotRange->{{xmin-dellx/4,xmax+dellx/4},{ymin-delly/4,ymax+delly/4}},FrameTicks->{xla,yla,xlatop,ylr},
TextStyle->{FontSize->22,FontWeight->"Plain",FontFamily->"Courier"},PlotLabel->StyleForm["8x8"],FrameLabel->{StyleForm["Re M^2"],StyleForm["Im M^2"]},FrameStyle->Thickness[0.007]]

SetDirectory["~/courses/cp09/week7"];
Export["8by8.eps",final,"EPS",ImageSize->500,
ImageOffset->{30,30}];

(*puts the three sets of points in 3 files
real   imaginary
real   imaginary  etc...*)
inpf=poles;ttt=cz[inpf];
dini=Dimensions[cz[inpf]][[1]];

stmp=OpenWrite["polesby8.d"];


Do[WriteString[stmp,ttt[[ll,1]]];WriteString[stmp,"  "];
WriteString[stmp,ttt[[ll,2]]];
WriteString[stmp,"\n"],{ll,1,dini}];

Close[stmp];

inpf=zeros;ttt=cz[inpf];
dini=Dimensions[cz[inpf]][[1]];

stmp=OpenWrite["zerosby8.d"];


Do[WriteString[stmp,ttt[[ll,1]]];WriteString[stmp,"  "];
WriteString[stmp,ttt[[ll,2]]];
WriteString[stmp,"\n"],{ll,1,dini}];

Close[stmp];

inpf=singp;ttt=cz[inpf];
dini=Dimensions[cz[inpf]][[1]];

stmp=OpenWrite["singpby8.d"];


Do[WriteString[stmp,ttt[[ll,1]]];WriteString[stmp,"  "];
WriteString[stmp,ttt[[ll,2]]];
WriteString[stmp,"\n"],{ll,1,dini}];

Close[stmp];