
/* Postscript plotting routines */

#include "cx.h"
#include "circles.h"
#include "lim.h"
#include <math.h>
#include <string.h>
#include <stdio.h>

/* Defaults */
static double linewidth=LINEWIDTH;
static int box=BOX, win=WIN;
double diam, wxmax, wxmin, wymax, wymin;
static char clipstr[300]=CLIP;

/* Calls to change parameters */
void set_linewidth(w)
	double w;
{	 linewidth=w;
}

/* Reverse box option */
void ps_box()
{	box = (box==0);
}

/* Set window */
void ps_win(xmin,ymin,xmax,ymax)
	double xmin, ymin, xmax, ymax;
{
	char str[200];

	win=1;
	wxmin=xmin; wxmax=xmax;
	wymin=ymin; wymax=ymax;
	sprintf(str,
"newpath\n %.5lf %.5lf moveto %.5lf %.5lf lineto\n %.5lf %.5lf lineto %.5lf %.5lf lineto\nclosepath clip\n",
	xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin);
	strcat(clipstr,str);
}

/* Determine if line is inside window */
int inside_line(circ)
	circle circ;
{
	double x, y;
	x = meetx(circ,wymin);
	if(x >= wxmin && x <= wxmax) return(1);
	x = meetx(circ,wymax);
	if(x >= wxmin && x <= wxmax) return(1);
	y = meety(circ,wxmin);
	if(y >= wymin && y <= wymax) return(1);
	y = meety(circ,wxmax);
	if(y >= wymin && y <= wymax) return(1);
	return(0);
}

/* Determine if circle is inside window */
int inside(circ)
	circle circ;
{
	double r;

	r = circ.r;
	if (r < 0) return(inside_line(circ));
	if (circ.c.x+r < wxmin) return(0);
	if (circ.c.y+r < wymin) return(0);
	if (circ.c.x-r > wxmax) return(0);
	if (circ.c.y-r > wymax) return(0);
	return(1);
}

/* Clip circle */
void ps_clip_circle(c,r)
	complex c;
	double r;
{
	char str[80];
	sprintf(str,
	"newpath\n %.5lf %.5lf %.5lf 0 360 arc\n",
	c.x,c.y,r);
	strcat(clipstr,str);
	strcat(clipstr,"gsave stroke grestore clip\n");
}

/* Scale and box */
void ps_window(xmin,ymin,xmax,ymax)
	double xmin, ymin, xmax, ymax;
{
	int marginx, marginy;
	double s, sx, sy, tx, ty;

/* Window overrides */
	if(win)
	{	xmax = wxmax; xmin = wxmin;
		ymax = wymax; ymin = wymin;
	}
	diam = fabs(xmax);
	if(fabs(xmin) > diam) diam = fabs(xmin);
	if(fabs(ymax) > diam) diam = fabs(ymax);
	if(fabs(ymin) > diam) diam = fabs(ymin);

	sx = PAGEX/(xmax-xmin);
	sy = PAGEY/(ymax-ymin);
	s  = (sx < sy) ? sx : sy;
	marginx = (PAGEX-s*(xmax-xmin))/2;
	marginy = (PAGEY-s*(ymax-ymin))/2;

	printf("%%%%BoundingBox:  %d %d %d %d\n",
		CORNX+marginx,
		CORNY+marginy,
		CORNX+PAGEX-marginx,
		CORNY+PAGEY-marginy);
	printf("%d %d translate\n",
		CORNX+marginx,
		CORNY+marginy);
	printf("%.2lf setlinewidth\n",linewidth);
	printf("%%%%\n");
	printf("%%%% Change to complex coordinates\n");
	printf("%.5lf %.5lf scale\n",s,s);
	printf("currentlinewidth %.5lf div setlinewidth\n",s);
	printf("%.5lf %.5lf translate\n",-xmin,-ymin); 
	if(box)
	{	printf("%%%% Box around figure %%%%\n");
		printf(
"newpath\n%.5lf %.5lf moveto %.5lf %.5lf lineto\n %.5lf %.5lf lineto %.5lf %.5lf lineto\nclosepath stroke\n",
		xmin,ymin,xmin,ymax,xmax,ymax,xmax,ymin);
	}
	printf("%%%% Clipping box %%%%\n");
	printf("%s",clipstr);
	printf("newpath\n");
	wxmin=xmin;
	wxmax=xmax;
	wymin=ymin;
	wymax=ymax;
}

void ps_plot_circles(stack,nstack)
	circle stack[];
	int nstack;
{
	int i;
	double bord, r, size, xmin, xmax, ymin, ymax;

	if(nstack == 0) return;
	xmin = xmax = stack[0].c.x;
	ymin = ymax = stack[0].c.y;
	for(i=0; i<nstack; i++) 
	{	r=stack[i].r;
		if(r<0) r=0;
		if(stack[i].c.x-r < xmin) xmin = stack[i].c.x-r;
		if(stack[i].c.x+r > xmax) xmax = stack[i].c.x+r;
		if(stack[i].c.y-r < ymin) ymin = stack[i].c.y-r;
		if(stack[i].c.y+r > ymax) ymax = stack[i].c.y+r;
	}

	printf("%%%% Kleinian group limit set in plane\n");
	printf("%%%%\n");
	size = (xmax-xmin > ymax-ymin) ? xmax-xmin : ymax-ymin;
	bord = BORDER*size;
	ps_window(xmin-bord,ymin-bord,xmax+bord,ymax+bord);
	ps_def_circle();
	for(i=0; i<nstack; i++)
		if(inside(stack[i]))
		ps_circle(stack[i]);
}

void ps_plot_circles_s2(stack,nstack)
	circle stack[];
	int nstack;
{
	int i;
	double r;

	printf("%%%% Kleinian group limit set on sphere\n");
	printf("%%%%\n");
	ps_window(-1.1,-1.1,1.1,1.1);
	ps_def_ellipse();
	printf("0 0 1 0 360 arc stroke\n");
	for(i=0; i<nstack; i++)
		ps_s2_circle(stack[i]);
}

/* Plot circles on front and back */
void ps_see_thru(stack,nstack,gray)
	circle stack[];
	int nstack;
	double gray;
{
	static matrix refl = 
		{0.0, 0.0, 1.0, 0.0, 1.0, 0.0, 0.0, 0.0};
	circle circ;
	int i;

	printf("%%%% Kleinian group limit set on sphere\n");
	printf("%%%%\n");
	ps_window(-1.1,-1.1,1.1,1.1);
	ps_def_ellipse();
	printf("0 0 1 0 360 arc stroke\n");
	printf("%.2lf setgray\n",gray);
	for(i=0; i<nstack; i++)
	{	circ = image_circle(refl,stack[i]);
		ps_s2_circle(circ);
	}
	printf("0 setgray\n");
	for(i=0; i<nstack; i++)
		ps_s2_circle(stack[i]);
}

void ps_circle(circ)
	circle circ;
{
	if(circ.r <= 0) ps_line(circ);
		else 
	printf("%.5lf %.5lf %.5lf c\n",
		circ.c.x,circ.c.y,circ.r);
}

void ps_line(circ)
	circle circ;
{
	complex a, b, v;

	v.x    = 2*diam * cos(-PI*circ.r);
	v.y    = 2*diam * sin(-PI*circ.r);
	a = sub(circ.c,v);
	b = add(circ.c,v);
	printf("%.5lf %.5lf %.5lf %.5lf l\n",
		a.x,a.y,b.x,b.y);
}

void ps_open(argc,argv)
	int argc;
	char *argv[];
{
	int i;

	printf("%%!\n");
	printf("%%%% Lim Version 2.0\n");
	printf("%%%% Command line:\n");
        printf("%%%% ");
        for(i=0; i<argc; i++)
                printf("%s ",argv[i]);
	printf("\n");
}

void ps_close()
{
	printf("showpage\n");
}

/* Print circle on S^2 */
void ps_s2_circle(circ)
	circle circ;
{
	vector v;
	double r;

	v = circle_to_s2(circ,&r);
	ps_ellipse(v,r);
}
	
/* Convert circle on S^2 to ellipse arc in x-y plane */
/* Circle specified by normal vector (need not be unit) and radius */
void ps_ellipse(v,r)
	vector v;
	double r;
{
	double angle,c,end,l,rminor,start,theta,w,x,y,z;
	x = v.x; y=v.y; z=v.z;

	l = vabs(v);
	x = x/l;
	y = y/l;
	z = z/l;
	c = sqrt(1-r*r);
	rminor = z*r;
	if(x==0.0 && y==0.0) theta = 0;
	else
	theta = 180*atan2(y,x)/PI;
	angle = 0;
	if(x == 0.0 && y == 0.0 && z < 0.0) return;
	if(x != 0.0 || y != 0.0)
	{	w = c*z/(r*sqrt(x*x+y*y));	
		if( w < -1.0) return;
		if( w < 1.0) angle = 180*acos(w)/PI;
	}
	if(angle < 1.0)
	printf("%.4lf %.4lf %.4lf %.4lf %.0lf ef\n",
	c*x, c*y, rminor, r, theta);
		else
	printf("%.4lf %.4lf %.4lf %.4lf %.0lf %.0lf %.0lf es\n",
	c*x, c*y, rminor, r, theta, angle, 360-angle);
}

/* Define ps line and full circle */
void ps_def_circle()
{
	printf("/l {moveto lineto stroke} def\n");
	printf("/c {0 360 arc stroke} def\n");
}

/* Define ps angled ellipse */
void ps_def_ellipse()
{
printf("/ellipsedict 9 dict def\n");
printf("ellipsedict /mtrx matrix put\n");
printf("/ellipse \n");
printf("  { ellipsedict begin\n");
printf("    /endangle exch def\n");
printf("    /startangle exch def\n");
printf("    /originangle exch def\n");
printf("    /yrad exch def\n");
printf("    /xrad exch def\n");
printf("    /y exch def\n");
printf("    /x exch def\n");
printf("    /savematrix mtrx currentmatrix def\n");
printf("    x y translate\n");
printf("    originangle rotate\n");
printf("    xrad yrad scale\n");
printf("    0 0 1 startangle endangle arc\n");
printf("    savematrix setmatrix\n");
printf("   end\n");
printf("  } def\n");
printf("/es { ellipse stroke} def\n");
printf("/ef { 0 360 ellipse stroke} def\n");
}
