/* DLA Main Line */

#include "cx.h"
#include <math.h>
#include <stdio.h>
#define PI 3.14159265358979

#define SPIKERR 0.01
#define SPIKE 1.41421
#define MAXAGG 10000
#define FLAGG 1
#define NGEN 10
#define RADIUS 1.1
#define NCIRC 0
#define NPTS 200
#define PICEPS 0.01
#define FID "caps.data"

/* Capacity file */
static char fid[80]=FID;

/* Aggregate data */
complex theta[MAXAGG];
double a[MAXAGG];
static int nagg = 0;
static double piceps=PICEPS, radius=RADIUS,
	spike = SPIKE;
static int flagg=FLAGG, ngen=NGEN, npts=NPTS, ncirc=NCIRC;

/* Functions */
complex map(), f(), randt();
double capf(), cp(), cv(), diam(), dmap(), df(), finda();
double spread();

main(argc,argv)
	int argc;
	char *argv[];
{
	int i, n, seed;

	srandom(time(NULL)); 	/* Random seed */
	for(i=1; i<argc; i++)
	{	if(argv[i][0] != '-') usage();
		switch(argv[i][1])
	{
		case 'b':
		ps_box();
		break;

		case 'c':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%d",&ncirc);
		break;

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

		case 'f':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%s",fid);
		break;

		case 'g':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%d",&ngen);
		break;

		case 'n':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%d",&npts);
		break;

		case 'o':
		flagg=0;
		break;

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

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

		case 'S':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%d",&seed);
		srandom(seed);
		break;

		default:
		usage();
	}
	}
	grow();
	di(argc,argv);
}

eprintf(s) char *s; {fprintf(stderr,s);}

usage()
{
	eprintf("Usage:  dla [options] (defaults below)\n");
	eprintf("  -b [reverse box option] (boxed)\n");
	eprintf("  -c [number of outer circles] (0)\n");
	eprintf("  -e [picture epsilon] (0.01)\n");
	eprintf("  -f [capacities file name] (caps.data)\n");
	eprintf("  -g [number of generations] (10)\n");
	eprintf("  -o omit lump plot\n");
	eprintf("  -n [number of points on buds] (200)\n");
	eprintf("  -r [outer radius] (1.1)\n");
	eprintf("  -s [size of buds] (sqrt(2))\n");
	eprintf("  -S [random seed]\n");
	exit();
}

/* Display results */
di(argc,argv)
	int argc;
	char *argv[];
{
	double d;
	int plot;

	di_caps();
	plot= (ncirc != 0) || flagg;
	if(!plot) return;

	ps_open(argc,argv);
	if(ncirc > 0) d=1.1*diam(radius);
	else d=1.1*diam(1.0);
	ps_window(-d,-d,d,d);
	piceps = piceps*d;
	if(flagg) di_growth();
	if(ncirc > 0) di_outer(radius,ncirc);
	ps_close();
}

/* Display Riemann map out to given radius */
di_outer(radius,ncirc)
	double radius;
	int ncirc;
{
	int i;
	double r;

	for(i=1; i<=ncirc; i++)
	{	r = 1.0 +(radius-1.0)*i/ncirc;
		di_map(r);
	}
}

/* Draw Riemann map applied to circle */
di_map(radius)
	double radius;
{
	int i;
	complex z;
	double deriv, theta;

	ps_newline();
	theta = 0.0;
	while (theta < 2*PI)
	{	z=polar(radius,theta);
		ps_draw(map(z));
		deriv = dmap(z);
		theta += (deriv>1) ? piceps/deriv : piceps;
	}
	z = polar(radius,0.0);
	ps_draw(map(z));
}

/* Draw Riemann map applied to arc */
di_spike(radius,ang1,ang2)
	double radius,ang1,ang2;
{
	int i;
	complex z;
	double deriv, delta, theta;

	ps_newline();
	z=polar(radius,ang1);
	ps_draw(map(z));

	delta = (ang2-ang1)/npts;
	theta = ang1+delta;
	while (theta < ang2-delta)
	{	z=polar(radius,theta);
		ps_draw(map(z));
		theta += delta;
	}
	z = polar(radius,ang2);
	ps_draw(map(z));
}

di_caps()
{
	int i;
	double c;
	FILE *cap_file;

	cap_file = fopen(fid,"w");
	c = 1.0;
	fprintf(cap_file,"%4d.  %10.3lf\n",1,c);
	for(i=0; i<nagg; i++)
	{	c *= capf(a[i]);
		fprintf(cap_file,"%4d.  %10.3lf\n",i+2,c);
	}
	fclose(cap_file);
}

/* Display growing Riemann images */
di_growth()
{
	complex w, z;
	double ang1, ang2, d;
	int i, k;

	k=nagg;
	nagg=0;
	di_map(1.0);
	for(i=0; i<k; i++)
	{	ang1=arg(theta[i])-cp(a[i]);
		ang2=ang1+2*cp(a[i]);
		nagg=i+1;
		di_spike(1.0,ang1,ang2);
	}
}

/* Compute diameter of image */
double diam(r)
	double r;
{
	double s, t;
	complex z;
	int i;

	s=0.0;
	for(i=0; i<=100; i++)
	{	z=polar(r,2*i*PI/100);
		t=infnorm(map(z));
		if(t > s) s = t;
	}
	return(s);
}

/* Grow */
grow()
{
	int i;
	for(i=1; i<ngen; i++) aggregate();
	fprintf(stderr,"%d generations complete.\n",ngen);
}

/* Add to aggregate */
aggregate()
{
	if(nagg >= MAXAGG)
	{	eprintf("Exceeded max of %d\n",MAXAGG);
		exit();
	}
	theta[nagg]=randt();
	a[nagg]=finda((theta[nagg]));
	nagg++;
}
/* Find the right size for new lump */
double finda(t)
        complex t;
{
        double a, dist, fac;

        a=spike/dmap(t);
        fac = 2.0;
	dist=spread(t,a);
        do
        {       if(dist<spike) do
		{	a = a*fac;
			dist=spread(t,a);
		} while(dist<spike);
		else do
		{	a = a/fac;
			dist=spread(t,a);
		} while(dist>spike);
		fac = sqrt(fac);
        } while(fabs(dist-spike)/spike > SPIKERR);
        return(a);
}

/* Distance between endpoints of proposed lump */
double spread(t,a)
	complex t;
	double a;
{
	double delta, dist;
	complex z1, z2;

        delta=cv(a);
        z1=mult(t,polar(1.0,delta));
        z2=mult(t,polar(1.0,-delta));
        dist=cabs(sub(map(z1),map(z2)));
	return(dist);
}

/* Riemann mapping to outside of aggregate */
complex map(z)
	complex z;
{
	int i;
	for(i=nagg-1; i>=0; i--)
		z = f(z,theta[i],a[i]);
	return(z);
}

double dmap(z)
	complex z;
{
	int i;
	double d;

	d = 1.0;
	for(i=nagg-1; i>=0; i--)
		d *= df(z,theta[i],a[i]);
		z = f(z,theta[i],a[i]);
	return(d);
}

/* Random point on the circle */
complex randt()
{
	long random();
	complex t;

	t=polar(1.0,2*PI*random()/1073741824.0);
	return(t);
}
