(*intersection of a line through p1,p1 and vertical line through x*)
tlx[p1_,p2_,x_]:=If[p1[[1]]==p2[[1]],0,intersect[p1,p2-p1,{x,0},{0,1}][[2]]];
(*lattice points on vertical line through x
between a line through p1,p2 and a line through p3,p4*)
latline[p1_,p2_,p3_,p4_,x_]:=(y1=tlx[p1,p2,x]; y2=tlx[p3,p4,x];Tuples[{{x},Range[Ceiling[Min[y1,y2]],Max[y1,y2]]}]);
(*returns list of latice points contained on a line determined by pts*)
line2lat[pts_]:=(lpts={};Do[y=tlx[pts[[1]],pts[[2]],x]; If[IntegerQ[y],AppendTo[lpts,{x,y}]],{x,Ceiling[Sort[pts][[1]][[1]]],Sort[pts][[2]][[1]]}];If[pts[[1]][[1]]!=pts[[2]][[1]],lpts,Tuples[{{pts[[1]][[1]]},Range[Ceiling[Min[pts[[1]][[2]],pts[[2]][[2]]]],Max[pts[[1]][[2]],pts[[2]][[2]]]]}]]);
(*returns list of latice points contained in a triangle determined by pts*)
tri2lat[pts_]:=(tri=Sort[pts]; lpts={};If[tri[[1]][[1]]!=tri[[2]][[1]],Do[lpts=Union[lpts,latline[tri[[1]],tri[[2]],tri[[1]],tri[[3]],x]],{x,Ceiling[tri[[1]][[1]]],tri[[2]][[1]]}]];If[tri[[2]][[1]]!=tri[[3]][[1]],Do[lpts=Union[lpts,latline[tri[[2]],tri[[3]],tri[[1]],tri[[3]],x]],{x,Ceiling[tri[[2]][[1]]],tri[[3]][[1]]}]]; lpts);
(*returns list of latice points contained in a quadrilateral determined by pts*)
quad2lat[pts_]:=Union[tri2lat[pts[[{1,2,3}]]],tri2lat[pts[[{1,3,4}]]]];
(*returns list of latice points contained in a line/triangle/quadrilateral determined by pts*)
getlat[pts_]:=Switch[Length[Union[pts]],2,line2lat[Union[pts]],3,tri2lat[Union[pts]],4,quad2lat[pts]]

(*old/slow/worse version of getlat, leaving it here incase I need it,returns list of latice points contained in convex hull of a list of points*)
getlatOLD[pts_]:=(lpts={};poly=ConvexHullMesh[pts];
Do[If[RegionMember[poly,{x,y}],AppendTo[lpts,{x,y}]],{x,Ceiling[Min[pts[[All,1]]]],Floor[Max[pts[[All,1]]]]},{y,Ceiling[Min[pts[[All,2]]]],Floor[Max[pts[[All,2]]]]}];lpts);

(*scaled coefficients of a,bth derivitive of x^iy^j*)
d[i_,j_,a_,b_]:=Binomial[i,a]*Binomial[j,b];
(*scaled coefficients of all derivitives (up to order n) of x^iy^j*)
d[pt_,n_]:=(der={};Do[AppendTo[der,d[pt[[1]],pt[[2]],a,b]],
{b,0,n},{a,0,n-b}];der);
(*returns a matrix of d[pt_,n_] vectors for each pt in pts*)
dm[pts_,n_]:=(dmat={};
Do[AppendTo[dmat,d[pt,n]],{pt,pts}];
Transpose[dmat]);
(*Returns the maximum multiplicity of a polynomial made from getlat[pts] at the {1,1}, together with a polynomial giving the maximum multiplicity. sorted: {maxm mult,poly}*)
maxmult[pts_]:=(mm=0;lat=getlat[pts];
m=Floor[Max[lat[[All,1]]]]+Floor[Max[lat[[All,2]]]];If[Length[lat]<1||m==0,dat={0,If[Length[lat]<1,{{0,0}},lat],{0},{0}},
Do[If[MatrixRank[dm[lat,n-1]]<Length[lat],mm=n,Break[]],{n,1,m}];{mm,Factor[prettypolly[If[mm>0,NullSpace[dm[lat,mm-1]][[1]],Table[0,{Length[lat]}]],lat]]}]);

(*self explanatory*)
showfan[f_]:=Graphics[Map[Arrow[{{0,0},Normalize[#]}]&,f]];
(*area of polygon, but 0 if its a line or point*)
area2[pts_]:=If[Length[Union[pts]]<3,0,Area[Polygon[pts]]];
(*these turn the lattice points and coefficients from maxmult into readable polynomials*)
pretty[c_,pt_]:=If[c==0,0,c*X^pt[[1]]Y^pt[[2]]];
prettypolly[cs_,pts_]:=Total[MapThread[pretty,{cs,pts},1]];
(*sorts the vectors of the fan clockwise*)
clock[f_]:=SortBy[f,-N[Arg[Complex[#[[1]],#[[2]]]]]&];
(*these just shift points around*)
sh[pts_,x_,y_]:=Map[(#+{x,y})&,pts];
mv[pts_]:=sh[pts,-Min[pts[[All,1]]],-Min[pts[[All,2]]]];
(*rotates by pi/4 and normalizes so greatest nonzero entry is 1*)
perp[v_]:=(nrm=Max[Abs[v]];({{0,1},{-1,0}}.v)/nrm);
(*checks if they're parallel*)
ispar[v1_,v2_]:=v1[[1]]*v2[[2]]==v1[[2]]*v2[[1]];
(*find intersection of lines through pi with direction vi*)
intersect[p1_,v1_,p2_,v2_]:=If[ispar[v1,v2],"parallel lines don't intersect",(p1+v1*(Inverse[Transpose[{v1,v2}]].(p2-p1))[[1]])];
(*reorders the fan (assumed already clockwise) so the first 2 coorespond to the sides  the normal polygons should be scaled by*)
pickaside[f_]:=(ang=N[MapThread[VectorAngle,{f,RotateLeft[f]}],20];a2=MapThread[Max,{ang+RotateLeft[ang],ang+RotateRight[ang]}];m=Min[a2]; j=FirstPosition[a2,m]);
(*takes fan, outputs vectors in the direction of the sides that can be scaled using the length of the first 2 sides*)
fan2sides[f_]:=Map[perp,RotateLeft[f,pickaside[f]-1]];
(*takes side vectors, outputs vertices of polygon scaled by l1,l2 and shifted into first quadrant*)
sides2poly[{s1_,s2_,s3_,s4_},l1_,l2_]:=(p1=l1*s1;p2=p1+l2*s2; p3=intersect[p2,s3,{0,0},s4];mv[{{0,0},p1,p2,p3}]);
(*takes fan, outputs a normal polygon scaled to have side lengths l1,l2*)
fan2poly[f_,l1_,l2_]:=sides2poly[fan2sides[f],l1,l2];
(*takes fan, outputs points of a graph of max multiplicity, together with all the data from maxmult, and areas. sorted: {i,j,max mult,poly, self intersection} ({i,j,mult} being the coordinates of the point)
generates polygons for all sets of scaling factors {i,j} with i+j≤n*)
fan2graph[f_,n_]:=(mults={};Do[po=fan2poly[f,k,l]; If[Length[Union[po]]>1,bt=maxmult[po];AppendTo[mults,{k,l,bt[[1]],bt[[2]],2*area2[po]-bt[[1]]^2}]],{k,0,n},{l,0,n-k}];mults);
(*these take a crossection of the graph of fan2graph (along x+y=1) keeping all the data, organized: ({scaled coordinates},{original coordinates},poly,self intersection)*)
scale[{v1_,v2_,v3_,v4_,v5_}]:=({{v1/(v1+v2),v3/(v1+v2)},{v1,v2,v3},v4,v5});
fan2crosec[f_,n_]:=Map[scale,fan2graph[clock[f],n]]

(*takes data organized like fan2crosec and produces a list of all irreducible factors of all polynomials containned in the data*)
getirr[p_]:=DeleteCases[If[p[[2]][[3]]>0,Map[#[[1]]&,FactorList[p[[3]]]],{}],_Integer|X|Y];
irrlist[d_]:=Union[Flatten[Map[getirr,d]]];
(*find multiplicity of a polynomial at {1,1}*)
mult[p_]:=(mon=MonomialList[p/.{X->(X+1),Y->(Y+1)},{X,Y},"DegreeLexicographic"][[-1]];Exponent[mon,Y]+Exponent[mon,X]);
(*find smallest polynomial normal to f containing pts*)
bestfit[f_,pts_]:=(pol={};
Do[AppendTo[pol,intersect[MaximalBy[pts,f[[i]].#&][[1]],perp[f[[i]]],MaximalBy[pts,RotateLeft[f][[i]].#&][[1]],perp[RotateLeft[f][[i]]]]],{i,Length[f]}];pol);
(*get lattice points of all monomials in a polynomial (X^i*Y^j->{i,j})*)
poly2pts[p_]:=Map[#[[1]]&,CoefficientRules[p,{X,Y}]];
(*takes a fan and polynomial and uses bestfit[f,pol] to produce data similar to the output of fan 2 crosec but with added self intersection numbers for triangles cooresponding to the triangles created by removing a side of the fan*)
reversal[f_,pol_]:=(lp=poly2pts[pol]; rf=RotateLeft[clock[f],pickaside[clock[f]]-1]; bff=bestfit[rf,lp];s1=Max[Abs[bff[[4]]-bff[[1]]]]; s2=Max[Abs[bff[[2]]-bff[[1]]]];mu=mult[pol];  ar=area2[bff];ar1=If[ispar[rf[[2]],rf[[3]]]||ispar[rf[[3]],rf[[4]]]||ispar[rf[[4]],rf[[2]]],"not a polygon",area2[bestfit[rf[[2;;4]],lp]]];ar2=If[ispar[rf[[1]],rf[[3]]]||ispar[rf[[3]],rf[[4]]]||ispar[rf[[4]],rf[[1]]],"not a polygon",area2[bestfit[rf[[{1,3,4}]],lp]]];{{s1/(s1+s2),mu/(s1+s2)},{s1,s2,mu},pol,2*ar-mu^2,2*ar1-mu^2,2*ar2-mu^2});
(*puts everything together*)
makedata[f_,n_]:=(f2c=fan2crosec[f,n]; rev=Map[reversal[f,#]&,irrlist[f2c]];{f,f2c,rev})

(*just "simplifies" by representing rays with coprime coordinates*)
simpfan[f_]:=Map[(#/GCD[#[[1]],#[[2]]])&,f];
(*makes {{0,1},{1,0},{-1,-1}} with added rays*)
fans1[x_,y_]:=(f1={};Do[If[GCD[a,b]==1(**),
AppendTo[f1,clock[simpfan[{{0,1},{1,0},{-1,-1},{-a,b}(**)}]]]],{a,1,x},{b,1,y}];f1);
(*makes {{2,3},{1,-2},{-8,1}} with added rays*)
fans2[x_,y_]:=(f2={};Do[If[GCD[a,b]==1 && !(2*a==b)(**),
AppendTo[f2,clock[simpfan[{{2,3},{1,-2},{-8,1},{a,-b}(**)}]]]],{a,1,x},{b,1,y}];f2);
(*makes {{11,3},{-4,-3},{-2,3}} with added rays*)
fans3[x_,y_]:=(f3={};Do[If[GCD[a,b]==1 && !(3*a==2*b)(**),
AppendTo[f3,clock[simpfan[{{11,3},{-4,-3},{-2,3},{-a,b}(**)}]]]],{a,1,x},{b,1,y}];f3);
(*makes {{13,10},{-13,-9},{1,0}} with added rays*)
fans4[x_,y_]:=(f3={};Do[If[GCD[a,b]==1 (**),
AppendTo[f3,clock[simpfan[{{13,10},{-13,-9},{1,0},{a,b}(**)}]]]],{a,1,x},{b,1,y}];f3);
(*makes symetric fans*)
fanss[mi_,ma_]:=(sf={};Do[If[!(GCD[a,b]==1),Continue[]];AppendTo[sf,simpfan[{{a,b},{1,0},-{a,b},-{1,0}}]],{a,mi,ma},{b,mi,ma}];sf);
(*makes a random fan with a ray in each quadrant*)
ranfan[n_]:=simpfan[MapThread[Times,{RandomInteger[{1,n},{4,2}],{{1,1},{1,-1},{-1,-1},{-1,1}}}]];

(*
as an example of how one mught use the above functions to generate some data, the following code would generate 50 random fans and save all the data needed to make their plots into dr.txt
*)
fans=Table[ranfan[20],{i,50}]
Do[makedata[fan,12]>>>dr.txt,{fan,ranfans}]