\\ this is a file written in pari,
\\ The aim is to get pari to write xfig pictures of
\\ domains for congruence subgroups etc, which is useful for
\\ illustrating ideas in talks.
\\ pari is used since it's free and easy to install.
\\ output is in xfig since then you can easily add notes to
\\ the picture, export as ps and ps_tex, and include in a latex file

\\ WARNING: has not been extensively tested.
\\ not using local variables in pari can cause problems

\\ search for "USER FUNCTIONS" for a description of how to use this

\\ here is an example:

\\? read("xfigpictures.pari");
\\? HeadAndCols("example.fig",["ff1199","ff3377","ff5555","ff7733","ff9900"])
\\? makeFilledImageOfPolygon("example.fig",["oo",-1/2+I*sqrt(3)/2,1/2+I*sqrt(3)/2],vector(10,i,s*t^(i-5)),vector(10,i,i+20),2,0)
\\? makeFilledImageOfPolygon("example.fig",["oo",-1/2+I*sqrt(3)/2,1/2+I*sqrt(3)/2],vector(5,i,t^(i-3)),vector(5,i,i+31),2,0)

\\ open the file example.fig with the xfig program to view

\\ the above example uses t,s,i as defined next:
\\ some standard matrices it's useful to have around: 

t=[1,1;0,1]; i=[1,0;0,1]; s=[0,-1;1,0];
\\ remove this line from the file if you want to use these letters for something else!


\\ change the following number, "picScale" to get different sized pictures
\\ this is the number of pixels for one unit, or something
\\ like that; in xfig, 1inch is usually 1200 pixels.

picScale=3000.0;


\\ TO DO: add functions for drawing farey symbol pictures
\\ including labels


\\ ===========================================

\\ header for xfig file:
{
makeheader(file)=
\\ see xfig webpages for description of xfig
\\ most of these comments will vanish from xfig file
\\ after resaving in xfig.
write(file,"#FIG 3.2
Landscape # orientation
Center # justification
Inches # units
Letter # paper size
100.00 # export magnification
Single # whether multi page 
-2     # this means no transparant background for gifs
# Fundamental domains picture
1200 2 # resolution and coordinate system (xfig default/only option)

")}

\\ scaling
rd(vv)=round(vv*picScale)

\\ the code for an xfig arc needs the center and three points
\\ on the arc.  The center can be computed from the three points.
\\ assuming the arcs are all geodesics, the center will always
\\ have y coord 0, assuming it's an arc, not a straight line
\\ the second point can also be determined from the end points,
\\ assuming the arc is a geodesic

{
makearc(file,p1,p3)=
\\ compute center, radius, etc of circle:
if(real(p1-p3)!=0,
 p1=p1*1.0;
 p3=p3*1.0;
 if(imag(p1-p3)==0,
  cx=real(p1+p3)/2,
  perpline = I*(p1-p3);
 if(imag(perpline)!=0,
    la=imag(p1+p3)/imag(perpline)/2);
  cx=(p1+p3)/2 - la*perpline;
 );
 radius=sqrt(norm(p1-cx));
 tt=(p1+p3)/2-cx;
\\ this is mid point on arc
 p2=cx+radius*tt/sqrt(norm(tt));
\\ this is the centre point
 cx=rd(real(cx));
 x1=rd(real(p1));y1=rd(imag(p1));
 x2=rd(real(p2));y2=rd(imag(p2));
 x3=rd(real(p3));y3=rd(imag(p3));
 thedirection=if(real(p3-p1)>0,0,1);
 write(file,"
 5 1 0 1 0 7 50 -1 -1 0 0 ",thedirection," 0 0 ",
 cx," ",0," ",
 x1," ",-y1," ",
 x2," ",-y2," ",
 x3," ",-y3,"\n"
 ),
 write(file,"
 2 1 0 1 0 0 50 -1 20 0 0 0 -1 0 0 2 \n",
	 rd(real(p1))," ",-rd(imag(p1))," ",rd(real(p3))," ",-rd(imag(p3)))
);
}


\\ this is for making the points of an arc for
\\ making a polyline to be concatinated to other lines
\\ for filling.

{
makearcpoints(p1,p3)=
\\ compute center of circle:
if(real(p1-p3)!=0,
 p1=p1*1.0;
 p3=p3*1.0;
 if(imag(p1-p3)==0,
  cx=real(p1+p3)/2,
  perpline = I*(p1-p3);
 if(imag(perpline)!=0,
    la=imag(p1+p3)/imag(perpline)/2);
  cx=(p1+p3)/2 - la*perpline;
 );
 radius=sqrt(norm(p1-cx));
\\ tt is the bisector line between
\\ the lines from the midpoint of the
\\ circle to the two arc points p1, p3
\\ it's used to compute the point p2 half
\\ way between p1 and p3 on the arc.
\\ the point is used for xfig arc, but
\\ it's not needed for using xfig polyline.
 tt=(p1+p3)/2-cx;
 if(tt!=0,
 p2=cx+radius*tt/sqrt(norm(tt)),
 p2=cx+I*radius);
\\ start and end angles:
 thet1 = -I*log((p1 - cx)/radius);
 thet2 = -I*log((p3 - cx)/radius);
\\ angle of arc:
 anglea = -I*log((p1-cx)/(p3-cx));
\\ there is a problem of whether to take pi or -pi
\\ hopefully this fixes things if they seem wrong; should be 
\\ a better way to do this.
 if(imag(exp(I*(thet1 - anglea/2)))<0,anglea=-anglea);
 vector(20,k,cx + radius*exp(I*(thet1 - (k-1)/20*anglea)))
 ,
 [p1,p3]
)}


\\ the following will just add one line to a file
\\ max is the maximum height of a point (inches),
\\ used for vertical lines.
{
makeimageofstandardtriangle(file,matrx,mx)=
write(file,"\n# triangle under matrix ",matrx,":");
\\ points as ratios:
e1p1=matrx*[1,0]~;
e1p2=matrx*[-1/2+I*(sqrt(3)/2),1]~;
e2p1=e1p2;
e2p2=matrx*[1/2+I*(sqrt(3)/2),1]~;
e3p1=e2p2;
e3p2=e1p1;
e1p2i=e1p2[1]/e1p2[2];
if(e1p1[2]!=0,
  e1p1i=e1p1[1]/e1p1[2];  
  makearc(file,e1p1i,e1p2i),
  makearc(file,I*mx + real(e1p2i),e1p2i)
  );
  e2p1i=e2p1[1]/e2p1[2];  e2p2i=e2p2[1]/e2p2[2];
  makearc(file,e2p1i,e2p2i);
e3p1i=e3p1[1]/e3p1[2];  
if(e1p1[2]!=0,
  e3p2i=e3p2[1]/e3p2[2];
  makearc(file,e3p1i,e3p2i),
  makearc(file,e3p1i,I*mx + real(e3p1i))
  )
}

\\ for filling a triangle, use polyline

{
makeImageofTrianglePoints(matrx,mx)=
\\ points as ratios:
e1p1=matrx*[1,0]~;
e1p2=matrx*[-1/2+I*(sqrt(3)/2),1]~;
e2p1=e1p2;
e2p2=matrx*[1/2+I*(sqrt(3)/2),1]~;
e3p1=e2p2;
e3p2=e1p1;
e1p2i=e1p2[1]/e1p2[2];
if(e1p1[2]!=0,
  e1p1i=e1p1[1]/e1p1[2];  
  listofpoints=makearcpoints(e1p1i,e1p2i),
  listofpoints=makearcpoints(I*mx + real(e1p2i),e1p2i)
  );
  e2p1i=e2p1[1]/e2p1[2];  e2p2i=e2p2[1]/e2p2[2];
  listofpoints=concat(listofpoints,makearcpoints(e2p1i,e2p2i));
e3p1i=e3p1[1]/e3p1[2];  
if(e1p1[2]!=0,
  e3p2i=e3p2[1]/e3p2[2];
  listofpoints=concat(listofpoints,makearcpoints(e3p1i,e3p2i)),
  listofpoints=concat(listofpoints,makearcpoints(e3p1i,I*mx + real(e3p1i)))
  );
listofpoints
}


\\ for a filled in polygon in upperhalf plane:
{
makeImageofPolyPoints(polygon,matrx,mx)=
\\ points as ratios:
listofpoints=[];
polygona=concat(polygon,[polygon[1]]);
for(kk=1,length(polygon),
\\ assume all points of polygon are different,
\\ so infinity does not occur twice in a row.
     if(polygona[kk]=="oo",
        end1a=matrx*[1,0]~,
        end1a=matrx*[polygona[kk],1]~);
     if(polygona[kk+1]=="oo",
        end2a=matrx*[1,0]~,
        end2a=matrx*[polygona[kk+1],1]~);
     if(end1a[2]==0,
        end1=real(end2a[1]/end2a[2]) + I*mx,
        end1=end1a[1]/end1a[2]);
     if(end2a[2]==0,
        end2=real(end1a[1]/end1a[2]) + I*mx,
        end2=end2a[1]/end2a[2]);
     listofpoints=concat(listofpoints,makearcpoints(end1,end2))
);
listofpoints
}

\\ this takes a list of complex numbers and puts them
\\ into x y coordinates in an xfig file.
{
pointsToXfig(file,colour,listofPoints)=
 write1(file,"2 3 0 1 0 ",colour," 50 -1 20 0 0 0 -1 0 0 ",
       length(listofPoints)+1,"\n");
 for(j=1,length(listofPoints),
     write1(file," ",rd(real(listofPoints[j]))," ",-rd(imag(listofPoints[j]))));
 write1(file," ",rd(real(listofPoints[1]))," ",-rd(imag(listofPoints[1])));
 write(file,"\n");
}

\\ this is for making images of the standard triangle to a file
\\ the standard triangle is the one with vertices at infinity, and
\\ +/- 1/2 + I*sqrt(3)/2.  
\\ I might delete these parts with the standard triangle later,
\\ and just have everything done using polyline instead of arcs.
{
makeimageSofstandardtriangle(file,matrices,mx)=
 makeheader(file);
 for(i=1,length(matrices),
     makeimageofstandardtriangle(file,matrices[i],mx)
 )
}

\\ this is for filled images of the standard triangle
{
makeFilledimageSofstdTri(file,matrices,colours,mx)=
 makeheader(file);
 for(jj=1,length(matrices),
     listofpointsA=makeImageofTrianglePoints(matrices[jj],mx);
     pointsToXfig(file,colours[jj],listofpointsA)
 )
}


\\ USER FUNCTIONS BELOW:
\\ this is the main function of this file:

\\ For the following function:
\\ file: string, probably ending in ".fig", assumed to be empty.
\\ polygon: list of points in upperhalf plane
\\ matrices: list of matrices to apply to the polygon
\\ colours: list of colours, which should be a list of the
\\  same length as the list of matrices, and should consist of
\\  numbers between 0 and 31, assuming you use the standard xfig colours
\\  (see xfig program for the palet of what they are) 
\\ mx: instead of plotting infinity, this is used for max height of
\\  imaginary part of point plotted, e.g., 2 is a reasonable value.
\\ Head: 0 or 1 depending on whether a header is required (1 means it is)
\\  this should be 1 the first time the file is written to, 0 after
\\  the file should be empty the first time; we don't delete any file,
\\  only add to them.


{
makeFilledImageOfPolygon(file,polygon,matrices,colours,mx,Head)=
  if(Head,makeheader(file));
  for(jj=1,length(matrices),
     listofpointsA=makeImageofPolyPoints(polygon,matrices[jj],mx);
     pointsToXfig(file,colours[jj],listofpointsA)
 )
}


\\ This function makes the header and also a list of user defined colours.
\\ if you want more than the standard 32 colours, you need this, then
\\ use makeFilledImageOfPolygon with Head set to 0 subsequently.
\\ colours should be entered as a hexidecimal string, corresponding to
\\ red blue green, e.g., "ff8800" for orange.

{
HeadAndCols(file,colourlist)=
makeheader(file);
for(ii=1,length(colourlist),write(file,"0 ",31+ii," #",colourlist[ii]))
}

\\ example usage:
\\ ? read("xfigpictures.pari");
\\ ? makeimageSofstandardtriangle("example1.fig",[t*s*t*t*s,s,s*t,s*t*t*t*t,s*t*t*t,t*s*t*t,t,i,t*t],2)
\\ open file with xfig to see picture.

\\ another example:
\\ ? HeadAndCols("example2.fig",["ff1100","ff3300","ff5500","ff7700","ff9900"])
\\ ? makeFilledImageOfPolygon("example2.fig",[0,I,2],[i,t,s,s*t],[32,33,34,35],2,0)

\\ another example:
\\ ? HeadAndCols("example3.fig",["ff1100","ff3300","ff5500","ff7700","ff9900"])
\\ ? makeFilledImageOfPolygon("example3.fig",["oo",I,1/2+I*sqrt(3)/2],[i,s,t,s*t,t*s],[1,2,3,4,5],2,0)



\\ other useful functions:

\\ function for multiplying a matrix by a list of fractions:

{
mult1(matrx,frc)=
 v1=if(frc==oo||frc=="oo",matrx*[1,0]~,matrx*[frc,1]~);
 if(v1[2]==0,"oo",v1[1]/v1[2])
}


mult(matrx,poly)=vector(length(poly),i,mult1(matrx,poly[i]))