 
/* Stable:  draws stable manifold of K3 dynamical systems */
/* Writes to stdout a Postscript file */

#include "k3.h"
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

static double centerx = CENTERX, centery = CENTERY, radius = RADIUS;
static double length=LENGTH, sep = SEP;
static int drawtop = 1, drawbot = 1, verbose=VERBOSE, spin=0;

/* Very large arrays */
point ps[PSMAX];

int main(argc,argv)
	int argc;
	char *argv[];

{
	int i,n;
	double a,b;

	for(i=1; i<argc; i++)
	{	if(argv[i][0] != '-') usage();
		switch(argv[i][1])
	{
		case 'b':
		drawtop = 0;
		break;

		case 'c':
		if(argc <= (i+2)) usage();
		sscanf(argv[++i],"%lf",&centerx);
		sscanf(argv[++i],"%lf",&centery);
		break;

		case 'e':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&sep);
		setsep(sep);
		break;

		case 'l':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&length);
		break;

		case 'p':
		if(argc <= (i+2)) usage();
		sscanf(argv[++i],"%lf",&a);
		sscanf(argv[++i],"%lf",&b);
		setk3(a,b);
		break;

		case 'q':
		verbose = 0;
		break;

		case 'r':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&radius);
		break;

		case 's':
		spin = 1;
		break;

		case 't':
		drawbot = 0;
		break;

		default:
		usage();
	}
	}

	ps_open(argc,argv);
	ps_window(centerx,centery,radius);
	draw_manifold();
	ps_close();
	exit(0);
}

void usage()
{
	fprintf(stderr,"Usage:  stable [options]\n");
	fprintf(stderr,"  -b          draw only bottom\n");
	fprintf(stderr,"  -c [x y]    window center\n");
	fprintf(stderr,"  -e [eps]    accuracy\n");
	fprintf(stderr,"  -l [length] length of stable manifold drawn\n");
	fprintf(stderr,"  -p [a b]    K3 parameters\n");
	fprintf(stderr,"  -q          quiet\n");
	fprintf(stderr,"  -r [r]      window radius\n");
	fprintf(stderr,"  -s          spin picture\n");
	fprintf(stderr,"  -t          draw only top\n");
	fprintf(stderr,"Postscript file written to stdout\n");
	exit(1);
}

point initial_point()
{
	point p;

	p = lift(sep/10,-sep/10);
	return(p);
}

void draw_manifold()
{
	int i, n, fpath();
	double tot, t;

	n=1;

	ps[0] = initial_point();
	ps[1] = f(ps[0]);
	draw_path(n);
	tot=pathlength(n);

	while(length > tot)
	{	n=fpath(ps,n);
		tot = tot + pathlength(n);
		di(ps,n);
		draw_path(n);
	}
	if(verbose) 
		fprintf(stderr,"Total path length:  %.2lf\n",tot);
}

void draw_path(n)
	int n;
{
	int i;
	point p, q;

	for(i=0; i<n; i++)
	{	p = ps[i]; q = ps[i+1];
/*		if(rotate) {p=rotate(p); q=rotate(q);} */
		if(drawtop &&  top(p)) ps_line(p,q);
		if(drawbot && !top(p)) ps_line(p,q);
	}
}

double pathlength(n)
	int n;
{
	int i;
	double d;

	d=0;
	for(i=0; i<n; i++) d = d + dist(ps[i],ps[i+1]);
	return(d);
}

void di(ps,n)
{
	double d;

	if(verbose)
	{	d=pathlength(n);
		fprintf(stderr,"Points: %8d; Length %.2lf\n",
			n,d);
	}
}
