#define SWITCH
/*	laplace - solve Laplace's equation for specified boundary conditions
				in two dimensions

	history...
		28 Jun 86	-g switch sets edge of grid to zero potential, 
					-m switch specifies a margin around the specified boundaries.
		21 Jun 86	Allowing boundary to be between grid points.
		31 May 86	showing boundries and calculated potentials are optional.

	approximate performance...

h=w=32, n=64
Current time is 20:58:25.76

Current time is 21:00:03.36


h=w=64, n=128
Current time is 21:01:11.29

Current time is 21:07:08.87

...producing 42K file

		...on Z-100 with V20 at 7.5 MHz, all files on ramdisk.

	notes...
		need to update above performance.
		should investigate recursive subdivision by factors other than 2.
		needs -m switch to specify a margin (for expanding the x & y range)
		can't handle boundary crossing edge of grid: referencing nonexistant
			grid points.
		needs option for cylindrical symmetry

*/
#define DESMET			/* deleting this only disables the timing function */

#include <stdio.h>
#include <math.h>

#define VERSION "0.2"
#define MAX_CROSS 300		/* max # times boundary can cross grid */
#define ENTRIES 300			/* maximum # points in boundary point curves */
#define MAXBOUNDARIES 50	/* maximum # separate curves in boundary */
#define PERMITTIVITY 8.85418e-12	/* permittivity of free space, farad/m */
#define SCALE 256
#define MAX(a,b) ((a)>(b)?(a):(b))

typedef struct
	{int grid;	/*	cell number */
	int direct;		/*	direction: 0, 1, 2, 3  for  +x, -x, +y, -y  */
	double dist;	/*	distance from cell center and boundary crossing point 
						(fraction of grid spacing)  */
	int volt;		/*	negative of boundary # (so it's an index into volt[]) */
	int dummy[4];	/*	...so sizeof(crossing) >= sizeof(boundary_rule)	*/
	} crossing;

typedef struct
	{int var[4];	/*	grid number on which this point depends	*/
	int coef[4];	/*	corresponding coefficient, scaled up by SCALE */
	} boundary_rule;

union 	/*	these can share space */
	{boundary_rule	b_rule[MAX_CROSS];
	crossing		cross_list[MAX_CROSS];
	} bo;

show_b(b) boundary_rule *b;
{	int i, sum;
	sum=0;
	for (i=0; i<4; i++) sum += b->coef[i];
	printf("rule %d: <%2d>*%d+<%2d>*%d+<%2d>*%d+<%2d>*%d  (ttl wt=%d)\n",
		b-bo.b_rule,
		b->var[0],b->coef[0],
		b->var[1],b->coef[1],
		b->var[2],b->coef[2],
		b->var[3],b->coef[3],
		sum);
	i=0;
}

char *c_label[4]={"right of","left of","above","below"};
show_c(c) crossing *c;
{	printf("boundary value %d is %f %s grid %d\n",
	c->volt, c->dist, c_label[c->direct], c->grid);
}


int lc;			/*	index of 1st unused entry in bo.cross_list	*/
int lb;			/*	index of 1st unused entry in bo.b_rule	*/

/*
	cell types are as follows...

			2 3 3 3 4
			9 1 1 1 5
			9 1 1 1 5
			9 1 1 1 5
			8 7 7 7 6

	...and type 0 cells have a specified potential (part of the
	boundary conditions) and are never updated.  There are special cell
	types are for boundaries where the normal gradient vanishes:

	         10  10  10  10
	       13              11
	       13              11
	       13              11
	         12  12  12  12
*/
int alpha=0, alphap1=1;
int *type; int *volt;
int debugging=0;
int show_resid=0;	/* nonzero if residuals are to be shown */
int grounded=0;		/* nonzero if edge of grid is to be grounded */
int show_boundaries=0;	/* nonzero if boundaries are to be shown */
int show_potentials=0;	/* nonzero if calculated potentials are to be shown */
int left_symmetric=0;	/* nonzero if x gradient vanishes on left edge of grid */
int lower_symmetric=0;	/* nonzero if y gradient vanishes on lower edge of grid */
int done=0;
long work=0;		/* number of cell updates so far */

double xmin=1.e30;
double xmax=-1.e30;
double ymin=1.e30;
double ymax=-1.e30;
double zmin=1.e30;
double zmax=-1.e30;
double margin=0.;		/* minimum amount to expand grid beyond given boundaries */
double *x;				/* pointer to data area for boundaries */
double *y;				/* pointer to data area for boundaries */
double xscale,yscale,zscale;
double edge;
double aspect,correct;
double p_voltage[MAXBOUNDARIES];
double residual();

int boundaries=0;		/* number of boundaries in input data */
int x_arguments=0;
int y_arguments=0;
int height=32;		/* default height of grid in cells */
int width=32;		/* default width of grid in cells */
int cycles=64;		/* default # relaxation cycles */
int n;				/* number of entries in x and y */
int index_array[MAXBOUNDARIES];	/* indexes into x and y */
int *p_data=index_array;
char outname[40];

FILE file;
FILE ofile;


main (argc,argv) int argc; char **argv;
{	int i,j,k,hw,col;
	double yval;
	char *s;
	read_data(argc,argv);
/*
	printf("there are %d points...\n",n);
	k=0;
	for (i=0; i<n; i++)
		{printf("%f %f\n",x[i],y[i]);
		if(i==p_data[k]) printf("       ...at voltage %f\n",p_voltage[k++]);
		}
*/
	if(s=index(argv[1],'.')) strncpy(outname,argv[1],s-argv[1]);
	else strcpy(outname,argv[1]);
	strcat(outname,".pot");
	ofile=fopen(outname,"w");
	if(ofile==0) {printf("can\'t create output file\n"); exit(1);}
	aspect=(ymax-ymin)/(xmax-xmin);
	if(height<width*aspect) height=width*aspect;
	else width=height/aspect;
	if(height<width*aspect)
		{correct=(ymax-ymin)*width/(double)height-(xmax-xmin);
		xmin -= correct/2.;
		xmax += correct/2.;
		}
	else
		{correct=(xmax-xmin)*height/(double)width-(ymax-ymin);
		ymin -= correct/2.;
		ymax += correct/2.;
		}
	hw=height*width;
	if(cycles<height) cycles=height;
	if(cycles<width) cycles=width;
	printf("x in (%g,%g) and y in (%g,%g)\n",xmin,xmax,ymin,ymax);
	type=malloc(hw*sizeof(int));
	volt=malloc((boundaries+hw)*sizeof(int)); 
	if(!type || !volt)
		{fprintf(stderr,"not enough workspace for %d by %d grid",height,width);
		exit(1);
		}
	volt += boundaries;		/* allow space for the boundary potentials */
	laplace(type,volt,width,height,cycles);
	if(show_potentials) show_volt(type,volt,width,height);
#ifdef FIXED
	charges(type,volt,width,height);
#endif
	printf("writing result to file %s\n",outname);
	if(-1==fprintf(ofile,"%10.4g%10.4g%5d    minimum and maximum x values, width\n",xmin,xmax,width)) exit(1);
	if(-1==fprintf(ofile,"%10.4g%10.4g%5d    minimum and maximum y values, height\n",ymin,ymax,height)) exit(1);
	yval=ymin;
	for (i=0; i<height; i++)
		{col=0;
		for (j=0; j<width; j++)
			{if(-1==fprintf(ofile,"%10.4g",volt[i*width+j]/zscale+zmin)) exit(1);
			if(++col>=7)
				{if(-1==fprintf(ofile,"\n")) exit(1);
				col=0;
				}
			}
		if(col)
			if(-1==fprintf(ofile,"\n")) exit(1);
		}
	fclose(ofile);
}

laplace (type,volt,width,height,cycles) int *type; int *volt,width,height,cycles;
{
#ifdef DESMET
	long start,tics();
#endif
	int i,j,k,hw;
	hw=height*width;
	if(height>width) i=height;
	else i=width;
	if(i>=8 && (width&1)==0 && (height&1)==0)
		{laplace(type,volt,width/2,height/2,cycles/2);
		i=hw; j=i/4;
		while(i)
			{i--; j--;
			volt[i]=volt[i-1]=volt[i-width]=volt[i-width-1]=volt[j];
			i--;
			if(i%width==0) i -= width;
			}
		}
	else
		{for (i=0; i<hw; i++) volt[i]=4096;
		}
	printf("height=%d, width=%d, cycles=%d\n",height,width,cycles);
#ifdef DESMET
	start=tics();
#endif
	fill(type,volt,width,height);
	k=0;
	if(show_boundaries)
		{for (i=height-1; i>=0; i--)
			{k=i*width;
			printf("%3d: ",i);
			for (j=0; j<width; j++)
				printf("%3d",type[k++]);
			printf("\n");
			}
		if(debugging) getchar();
		}
/*	show_volt(type,volt,width,height); */
	while(cycles>-1)
		{while(!csts() && cycles--)
			{if(!show_resid)
				printf("\015 %d cycles to go  ",cycles);
			relax(type,volt,width,height);
			work += hw;
			if(show_resid && ++done%5==0)
				printf("%D   %f\n",work,residual(type,volt,width,height));
			}
		if(cycles>-1)
			{printf("\nINTERRUPTED                 ");
			show_volt(type,volt,width,height);
			printf("continue calculation (Y/N)? ");
			if(toupper(getchar())!='Y')break;
			}
		}
#ifdef DESMET
	printf("\015 finished in %4.2f sec  \n",.01*(tics()-start));
#else
	printf("\015 finished               \n");
#endif
}

relax (type,v,w,h) int *type; int *v,w,h;
{	int i,wh;
	int hm1,wm1,j;
	boundary_rule *rule;
	long sum;

#ifdef SWITCH
	i=wh=w*h;
	while(i--)
		{if(*type>0)
			{switch(*type)
				{
				case  1: v[0]= (  v[-w] + v[-1] +     v[1] + v[w] + 2)/4; break;
				case  2: v[0]= (                      v[1] + v[w] + 1)/2; break;
				case  3: v[0]= (        2*v[-1] + 2*v[1]   + v[w] + 2)/5; break;
				case  4: v[0]= (          v[-1]            + v[w] + 1)/2; break;
				case  5: v[0]= (2*v[-w] + v[-1]          + 2*v[w] + 2)/5; break;
				case  6: v[0]= (  v[-w] + v[-1]                   + 1)/2; break;
				case  7: v[0]= (  v[-w] + 2*v[-1] + 2*v[1]        + 2)/5; break;
				case  8: v[0]= (  v[-w]         +     v[1]        + 1)/2; break;
				case  9: v[0]= (2*v[-w]         +   v[1] + 2*v[w] + 2)/5; break;
				case 10: v[0]= (2*v[ w]   + v[-1] + v[1]           + 2)/4; break;
				case 11: v[0]= (  v[-w] + 2*v[-1]          + v[ w] + 2)/4; break;
				case 12: v[0]= (            v[-1] + v[1] + 2*v[-w] + 2)/4; break;
				case 13: v[0]= (  v[-w]         + 2*v[1]   + v[ w] + 2)/4; break;
				}
			}
		else if (*type<0)
			{rule=&bo.b_rule[-*type];
			sum=0L;
#ifdef DEBUG
printf("\ngrid %d: ",wh-i-1); show_b(rule);
show_volt(type,volt,w,h);
#endif
			for (j=0; j<4; j++) sum +=	volt[rule.var[j]]*(long)rule.coef[j];
			v[0]=sum/SCALE;
#ifdef DEBUG
printf("setting potential to %5.2f\n",v[0]/zscale+zmin);
#endif
			}
		v++; type++;
		}
#else
	hm1=h-1;
	wm1=w-1;
	if(*type++)
		v[0]= (                      v[1] + v[w] + 1)/2;
	v++;
	for (j=1; j<wm1; j++)
		{if(*type++)
			v[0]= (           v[-1] + v[1] + v[w] + 1)/3;
		v++;
		}
	if(*type++)
		v[0]= (           v[-1]           + v[w] + 1)/2;
	v++;
	for (i=1; i<hm1; i++)
		{
		if(*type++)
			v[0]= (v[-w]            + v[1] + v[w] + 1)/3;
		v++;
		for (j=1; j<wm1; j++)
			{if(*type++)
				v[0]= (v[-w] + v[-1] + v[1] + v[w] + 2)/4;
			v++;
			}
		if(*type++)
			v[0]= (v[-w] + v[-1]           + v[w] + 1)/3;
		v++;
		}
	if(*type++)
		v[0]= (v[-w]            + v[1]           + 1)/2;
	v++;
	for (j=1; j<wm1; j++)
		{if(*type++)
			v[0]= (v[-w] + v[-1] + v[1]           + 1)/3;
		v++;
		}
	if(*type++)
		v[0]= (v[-w] + v[-1]                     + 1)/2;
	v++;

#endif
}

double residual (type,volt,w,h) int *type; int *volt,w,h;
{	int i,hw,n;
	int hm1,wm1,j;
	long r,d;

	n=0;
	i=hw=w*h;
	r=0L;
	while(i--)
		{switch(*type++)
			{case 0: n++; break;
			case 1: d=volt[0]-(volt[-w] + volt[-1] + volt[1] + volt[w] + 2)/4; r+=d*d; break;
			case 2: d=volt[0]-(                      volt[1] + volt[w] + 1)/2; r+=d*d; break;
			case 3: d=volt[0]-(           volt[-1] + volt[1] + volt[w] + 1)/3; r+=d*d; break;
			case 4: d=volt[0]-(           volt[-1]           + volt[w] + 1)/2; r+=d*d; break;
			case 5: d=volt[0]-(volt[-w] + volt[-1]           + volt[w] + 1)/3; r+=d*d; break;
			case 6: d=volt[0]-(volt[-w] + volt[-1]                     + 1)/2; r+=d*d; break;
			case 7: d=volt[0]-(volt[-w] + volt[-1] + volt[1]           + 1)/3; r+=d*d; break;
			case 8: d=volt[0]-(volt[-w]            + volt[1]           + 1)/2; r+=d*d; break;
			case 9: d=volt[0]-(volt[-w]            + volt[1] + volt[w] + 1)/3; r+=d*d; break;
			}
		volt++;
		}
	if (hw<=n)
		return 1.;
	return sqrt(((double)r)/(hw-n));
}

#ifdef FIXED
/*
	this no longer works...could be changed to solve linear system for
	the gradient at center of each grid point next to a boundary and
	integrate, but it's easier to let CONTOUR calculate the integral.
*/
charges (type,volt,width,height) int *type; int *volt,width,height;
{	int v,i,plate=0,out=0,hw;
	long charge;
	double C,Z,q[2],plate_voltage[2],sum=0.;
	hw=height*width;
	printf("\nComputing surface integral of grad V over each boundary\n");
	printf("(integral of -grad V dot N, where N is normal to the boundary)...\n");
	printf("\nboundary    potential    surface integral     charge\n");
	while(1)
		{for (i=0; i<hw; i++)
			if(type[i]==0)
				break;
		if(i==hw) break;
		v=volt[i];
		plate_voltage[plate&1]=v/zscale+zmin;
		charge=0.;
		if(i%width) charge += v-volt[i-1];
		if(i>=width) charge += v-volt[i-width];
		if(i<hw-width) charge += v-volt[i+width];
		if(i%width != width-1) charge += v-volt[i+1];
		type[i]=23;
		for (i++; i<hw; i++)
			{if(type[i]==0 && v==volt[i])
				{if(i%width) charge += v-volt[i-1];
				if(i>=width) charge += v-volt[i-width];
				if(i<hw-width) charge += v-volt[i+width];
				if(i%width != width-1) charge += v-volt[i+1];
				type[i]=23;		/* don't look at this cell again */
				}
			}
		sum += (q[plate&1]=(charge/zscale)*PERMITTIVITY);
		printf("%5d %14.2g %15.3g %15.3g\n",
			plate,v/zscale+zmin,q[plate&1]/PERMITTIVITY,q[plate&1]);
		plate++;
		}
	printf("    sum (should be zero)%12.3g %15.3g\n\n",sum/PERMITTIVITY,sum);
	if(plate==2)
		{q[0]=fabs(q[0]);
		q[1]=fabs(q[1]);
		C=MAX(q[1],q[0])/fabs(plate_voltage[1]-plate_voltage[0]);
		Z=1./C/2.9979e8;
		printf("characteristic impedance is %5.2f ohms\n",Z);
		printf("capacitance is %7.3g farad/m\n",C);
		printf("inductance is %7.3g henries/m\n",Z/2.9979e8);
		}
}
#endif
/*
	the equations relating a grid point and its four neighbors, in the
	absence of a boundary, is:
	
		( 0  1  1  1) (Ux   )   (Uright)
		( 0 -1  1  1) (Uy   ) = (Uleft )
 		( 1  0 -1  1) (Uxx  )   (Uup   )
		(-1  0 -1  1) (Uhere)   (Udown )

	where	Ux is h*(d/dx U)
			Uy is h*(d/dy U)
			Uxx is .5*h**2*(d^2/dx^2 U))
*/

double lhs[16]=
	{0., 1., 1., 1.,
	 0.,-1., 1., 1.,
	 1., 0.,-1., 1.,
	-1., 0.,-1., 1.};

#define LEFT_SIDE(g)	((g)%w==0)
#define RIGHT_SIDE(g)	((g)%w==w-1)
#define TOP_SIDE(g)		((g)<w)
#define LOWER_SIDE(g)	((g)>=hw-w)

	/*	fill type and voltage arrays based on user's boundary conditions */
fill (type,volt,w,h) int *type; int *volt,w,h;
{	boundary_rule *ru;
	crossing *this,*last;
	double invert(), decomp();
	double a[4][4], ai[4][4], cond, condp1, *p, f;
	double xt,yt,xt2,yt2;
	int i,j,g,weight,zero;
	int compare();
	int i,k,hw,zz;

	printf("establishing boundary conditions\n");
	hw=h*w;
	xscale=(w-1)/(xmax-xmin);
	yscale=(h-1)/(ymax-ymin);	/* same as xscale */
	edge=1./xscale;
/*	printf("xscale=%f, yscale=%f, edge=%f\n",xscale,yscale,edge); */
	zscale=32760./5./(zmax-zmin);

	for (i=0; i<hw; i++)	type[i]=1;		/* interior */
	if(grounded)
		{zero=floor((0.-zmin)*zscale+.499);
		for (i=w; i<hw; i+=w)	{type[i]=0; volt[i]=zero;}	/* left */
		for (i=0; i<w; i++)		{type[i]=0; volt[i]=zero;}	/* top */
		for (i=w-1; i<hw; i+=w)	{type[i]=0; volt[i]=zero;}	/* right */
		for (i=hw-w; i<hw; i++)	{type[i]=0; volt[i]=zero;}	/* bottom */
		}
	else
		{if(left_symmetric)
			for (i=w; i<hw; i+=w)	type[i]=13;		/* left */
		else
			for (i=w; i<hw; i+=w)	type[i]=9;		/* left */
		for (i=0; i<w; i++)		type[i]=3;			/* top */
		for (i=w-1; i<hw; i+=w)	type[i]=5;			/* right */
		if(lower_symmetric)
			for (i=hw-w; i<hw; i++)	type[i]=12;		/* bottom */
		else
			for (i=hw-w; i<hw; i++)	type[i]=7;		/* bottom */
		type[0]=2; type[w-1]=4; type[hw-w]=8; type[hw-1]=6;		/* corners */
		}
	i=k=0;
	lc=2;	/* bo.cross_list is empty */
	while(i<n)
		{xt=x[i]-xmin;
		yt=y[i]-ymin;
		zz=floor((p_voltage[k]-zmin)*zscale+.499);
		volt[-k-1]=zz;
		do	{xt2=x[i]-xmin;
			yt2=y[i]-ymin;
			mark(w,h,xt,yt,xt2,yt2,-k-1);
/*
			printf("%f,%f -> ",(x[i]-xmin)*xscale,(y[i]-ymin)*yscale);
			printf("mark(%d,%d, %d,%d, %d,%d, %d)\n",w,h,xx,yy,xx2,yy2,zz); 
*/
			xt=xt2; yt=yt2; 
			i++;
			} while(i<=p_data[k]);
		k++;
		}
/*	printf("the first of the %d crossings before sorting\n",lc-1);
	for (i=2; i<9; i++) show_c(bo.cross_list[i]); getchar();	*/
/*----------------  sort list  -------------*/
	hsort(&bo.cross_list[2],lc-2,sizeof(crossing),&compare);
/*----------------  solve each problem  -------------*/
#ifdef DEBUG
	printf("crossings after sorting\n");
	for (i=2; i<9; i++) show_c(bo.cross_list[i]); getchar();
#endif
	lb=1;	/* bo.b_rule is empty */
	this=&bo.cross_list[2];
	last=&bo.cross_list[lc];
	while(this<last)
		{p=a;
		for (i=0; i<16; i++) p[i]=lhs[i];
		g=this->grid;
#ifdef DEBUG
printf("rule %d, boundary conditions for grid point %d\n",lb,g);
#endif
		volt[g]=floor((p_voltage[-this->volt-1]-zmin)*zscale+.499);
		ru=&bo.b_rule[lb];
		ru->var[0]=g+1;
		ru->var[1]=g-1;
		ru->var[2]=g+w;
		ru->var[3]=g-w;
		do	{j=this->direct;
#ifdef DEBUG
printf("        "); show_c(this);
#endif
			f=this->dist;
			ru->var[j]=this->volt;
			a[j][0] *= f;
			a[j][1] *= f;
			a[j][2] *= f*f;
			this++;
			} while (this->grid == g);
/*		show_m("boundary condition matrix",a); getchar();	*/
		cond=invert(4,4,a,ai);
/*		printf("condition number is %f\n",cond);
		show_m("matrix inverse",ai);	*/
		condp1=cond+1.;
		if(cond==condp1)
			{fprintf(stderr,
				"can\'t solve linear system for boundary conditions\n");
			exit(1);
			}
		for (i=0; i<4; i++) ru->coef[i]=ai[3][i]*SCALE+.5;
		if(LOWER_SIDE(g))
			{if(lower_symmetric)
				ru->coef[2] += ru->coef[3];
			else
				{weight=ru->coef[2]=(ru->coef[2] + ru->coef[3])/4;
				weight = SCALE - weight + ru->coef[0]+ru->coef[1];
				ru->coef[0] += weight/2;
				ru->coef[1] += weight - weight/2;
				}
			ru->coef[3]=0;
			}
		if(LEFT_SIDE(g))
			{if(left_symmetric)
				ru->coef[0] += ru->coef[1];
			else
				{weight=ru->coef[0]=(ru->coef[0] + ru->coef[1])/4;
				weight = SCALE - weight + ru->coef[2] + ru->coef[3];
				ru->coef[2] += weight/2;
				ru->coef[3] += weight - weight/2;
				}
			ru->coef[1]=0;
			}
		if(TOP_SIDE(g))
			{weight=ru->coef[3]=(ru->coef[2] + ru->coef[3])/4;
			weight = SCALE - weight + ru->coef[0]+ru->coef[1];
			ru->coef[1] += weight/2;
			ru->coef[0] += weight - weight/2;
			ru->coef[2]=0;
			}
		if(RIGHT_SIDE(g))
			{weight=ru->coef[1]=(ru->coef[1] + ru->coef[0])/4;
			weight = SCALE - weight + ru->coef[2]+ru->coef[3];
			ru->coef[2] += weight/2;
			ru->coef[3] += weight - weight/2;
			ru->coef[0]=0;
			}
#ifdef DEBUG
printf("        result is ");show_b(ru);
#endif
		type[g]=-lb;
		lb++;
		}

}

compare(a,b) crossing *a,*b;
{	return (a->grid - b->grid);
}

mark (width,height,x1,y1,x2,y2,v) int width,height,v; double x1,y1,x2,y2;
{	int q,dx,dy,ax,ay,decide,up,over,i;
	int nx,nx2,ny,ny2;
	double del,xe,ye,f,t;
/*
	if(x1<0||x1>=width||x2<0||x2>=width||y1<0||y1>=height||y2<0||y2>=height)
		{fprintf(stderr,"calling sequence error: mark(%d, %d,%d, %d,%d, %d)\n",
			width,x1,y1,x2,y2,v);
		fprintf(stderr,"(point outside 0<=x<%d or 0<=y<%d)\n",width,height);
		exit(1);
		}
*/
#ifdef DEBUG
printf("line from  %5.2f,%5.2f  to  %5.2f,%5.2f at potential %5.2f\n",
x1,y1,x2,y2,p_voltage[-v-1]);
#endif
	if(x2<x1) {t=x1; x1=x2; x2=t; t=y1; y1=y2; y2=t;}
	nx2=x2*xscale;
	del=x2-x1;
	for (nx=x1*xscale+1; nx<=nx2; nx++)
		{xe=nx*edge;
		ye=(y1*(x2-xe) + y2*(xe-x1))/del;
		ny=ye*yscale;
		f=ye/edge-ny;
#ifdef DEBUG
printf("     crosses y grid line at  %f,%f  ->  %d,%d at distance %f\n",xe,ye,nx,ny,f);
#endif
		register_crossing(width,nx,ny,f,2,v);
		register_crossing(width,nx,ny+1,1.-f,3,v);
		}
	if(y2<y1) {t=x1; x1=x2; x2=t; t=y1; y1=y2; y2=t;}
	ny2=y2*yscale;
	del=y2-y1;
	for (ny=y1*yscale+1; ny<=ny2; ny++)
		{ye=ny*edge;
		xe=(x1*(y2-ye) + x2*(ye-y1))/del;
		nx=xe*xscale;
		f=xe/edge-nx;
#ifdef DEBUG
printf("     crosses x grid line at  %f,%f  ->  %d,%d at distance %f\n",xe,ye,nx,ny,f);
#endif
		register_crossing(width,nx,ny,f,0,v);
		register_crossing(width,nx+1,ny,1.-f,1,v);
		}
}

register_crossing(width,nx,ny,f,direction,v)
int width,nx,ny,direction,v;
double f;
{	crossing *new;

	if(lc>=MAX_CROSS) {fprintf(stderr,"too many boundary crosspoints\n"); exit(1);}
	new=&bo.cross_list[lc++];
	new->grid=nx+ny*width;
	new->direct=direction;
	new->dist=f;
	new->volt=v;
}

show_volt (type,volt,w,h) int *type; int *volt,w,h;
{	int i,j,k;
	k=0;
printf("boundaries...");
for (i=1; i<=boundaries; i++) printf("%d: %5.2f   ",i,volt[-i]/zscale+zmin);
printf("\n");
/*	printf("(press any key to abort printout)\n");
	for (i=0; i<h && !csts(); i++)
*/
	for (i=height-1; i>=0; i--)
		{k=i*w;
		printf("%3d: ",i);
		for (j=0; j<w; j++)
			{printf("%6.2f",volt[k++]/zscale+zmin);
			/* printf("%6d",volt[k++]); */
			}
		printf("\n");
		}
	for (i=height-1; i>=0; i--)
		{k=i*w+j;
		printf("%3d: ",i);
		for (j=0; j<width; j++)
			printf("%3d",type[k++]);
		printf("\n");
		}
	if(debugging) getchar();
}

read_data(argc,argv) int argc; char **argv;
{	int i,j,length;
	double xx,yy,d,*pd,sum;
	char *s,*t;
#define BUFSIZE 200
	static char buf[BUFSIZE];

	x=malloc(ENTRIES*sizeof(double));
	y=malloc(ENTRIES*sizeof(double));
	if(x==0 || y==0) {fprintf(stderr,"can\'t allocate buffer"); exit();}
	if(argc>1 && streq(argv[1],"?")) help();
	if(argc<=1 || *argv[1]=='-') file=stdin;
	else
		{if(argc>1)
			{file=fopen(argv[1],"r");
			if(file==0) {printf("file %s not found\n",argv[1]); exit();}
			argc--; argv++;
			}
		else help();
		}
	argc--; argv++;
	while(argc>0)
		{i=get_parameter(argc,argv);
		argv=argv+i; argc=argc-i;
		}
	p_data[0]=-1;
	if(height<2||width<2||cycles<2)
		{fprintf(stderr,"parameter too small: height=%d, width=%d, cycles=%d",
			height,width,cycles);
		exit(1);
		}
	i=0;
	while(i<ENTRIES)
		{if(fgets(buf,BUFSIZE,file)==0) break;
		t=buf;
		while(*t && isspace(*t)) t++;
		if(*t == 0) continue;		/* skip blank lines */
		buf[strlen(buf)-1]=0;		/* zero out the line feed */
		if(buf[0]==';') {printf("%s\n",buf); continue;}  /* skip comment */
		sscanf(buf,"%F %F",&x[i],&y[i]);
		if (x[i]<xmin) xmin=x[i];
		if (x[i]>xmax) xmax=x[i];
		if (y[i]<ymin) ymin=y[i];
		if (y[i]>ymax) ymax=y[i];
		s=buf;                      /* start looking for label */
		while(*s==' ')s++;			/* skip first number */
		while(*s && (*s!=' '))s++;
		while(*s==' ')s++;			/* skip second number */
		while(*s && (*s!=' '))s++;
		while(*s==' ')s++;
		if((length=strlen(s))&&(boundaries<MAXBOUNDARIES))
			{if(*s=='\"') s++;           /* label in quotes */
			sscanf(s,"%F",&p_voltage[boundaries]);
			if(p_voltage[boundaries]<zmin) zmin=p_voltage[boundaries];
			if(p_voltage[boundaries]>zmax) zmax=p_voltage[boundaries];
			p_data[boundaries++]=i;
			}
		i++;
		}
	if(grounded)
		{if(0.<zmin) zmin=0.;
		if(0.>zmax) zmax=0.;
		}
	n=i;
	if(!boundaries || p_data[boundaries-1]!=n-1)
		{p_data[boundaries]=i-1;
		p_voltage[boundaries++]=0.;
		}
	xmin -= margin;
	xmax += margin;
	ymin -= margin;
	ymax += margin;
}


/* get_parameter - process one command line option
		(return # parameters used) */
get_parameter(argc,argv) int argc; char **argv;
{	int i;
	double temp;

	if(streq(*argv,"-d")) {debugging=1; return 1;}
	else if(streq(*argv,"-h"))
		{if((argc>1) && numeric(argv[1])) height=atoi(argv[1]);
		return 2;
		}
	else if(streq(*argv,"-w"))
		{if((argc>1) && numeric(argv[1])) width=atoi(argv[1]);
		return 2;
		}
	else if(streq(*argv,"-g"))
		{grounded++;
		return 1;
		}
	else if(streq(*argv,"-m"))
		{i=get_double(argc,argv,1,&margin,&margin,&margin);
		if(margin<0.) margin=0.;
		return i;
		}
	else if(streq(*argv,"-n"))
		{if((argc>1) && numeric(argv[1])) cycles=atoi(argv[1]);
		return 2;
		}
	else if(streq(*argv,"-r"))
		{show_resid++;
		return 1;
		}
	else if(streq(*argv,"-b"))
		{show_boundaries++;
		return 1;
		}
	else if(streq(*argv,"-p"))
		{show_potentials++;
		return 1;
		}
/*	else if(streq(*argv,"-a"))
		{i=get_double(argc,argv,1,&temp,&temp,&temp);
		alpha=temp;
		alphap1=alpha+1;
		return i;
		}
*/
	else if(streq(*argv,"-x"))
		{i=get_double(argc,argv,2,&xmin,&xmax,&xmax);
		x_arguments=i-1;
		return i;
		}
	else if(streq(*argv,"-y"))
		{i=get_double(argc,argv,2,&ymin,&ymax,&ymax);
		y_arguments=i-1;
		return i;
		}
	else if(streq(*argv,"-xs"))
		{left_symmetric++;
		return 1;
		}
	else if(streq(*argv,"-ys"))
		{lower_symmetric++;
		return 1;
		}
	else gripe(argv);
}

get_double(argc,argv,permitted,a,b,c)
int argc,permitted; char **argv; double *a,*b,*c;
{	int i=1;
	if((permitted--)>0 && (argc>i) && numeric(argv[i])) *a=atof(argv[i++]);
	if((permitted--)>0 && (argc>i) && numeric(argv[i])) *b=atof(argv[i++]);
	if((permitted--)>0 && (argc>i) && numeric(argv[i])) *c=atof(argv[i++]);
	return i;
}

int streq(a,b) char *a,*b;
{	while(*a)
		{if(*a!=*b)return 0;
		a++; b++;
		}
	return 1;
}

gripe_arg(s) char *s;
{	fprintf(stderr,"argument missing for switch %s",s);
	help();
}

gripe(argv) char **argv;
{	fprintf(stderr,*argv); fprintf(stderr," isn\'t a legal argument \n\n");
	help();
}

help()
{	fprintf(stderr,"laplace   version %s",VERSION);
	fprintf(stderr,"\nusage: laplace [file][options]\n");
	fprintf(stderr,"options are:\n");
/*	fprintf(stderr,"  -a  num            acceleration factor (default %d) \n",alpha); */
	fprintf(stderr,"  -b               display boundaries\n");
	fprintf(stderr,"  -g               edge of grid is grounded (potential 0)\n");
	fprintf(stderr,"  -p               display calculated potentials\n");
	fprintf(stderr,"  -r               calculate and display residuals\n");
	fprintf(stderr,"  -m  margin       amount to expand grid beyond given boundaries (default 0) \n");
	fprintf(stderr,"  -n  num          number of relaxation cycles (default %d) \n",cycles);
	fprintf(stderr,"  -h  num          height of grid in cells (default %d)\n",height);
	fprintf(stderr,"  -w  num          width of grid in cells (default %d)\n",width);
	fprintf(stderr,"  -x  [min [max]]  specify x range \n");
	fprintf(stderr,"  -y  [min [max]]  specify y range \n");
	fprintf(stderr,"  -xs              symmetric across the line x=xmin \n");
	fprintf(stderr,"  -ys              symmetric across the line y=ymin \n");
	exit();
}

numeric(s) char *s;
{	char c;
	while(c=*s++)
		{if((c<='9' && c>='0') || c=='+' || c=='-' || c=='.') continue;
		return 0;
		}
	return 1;
}

#ifdef DESMET

/*
	tics - report hundreths of a second since midnight
*/

long tics()
{	int hr,min,sec,hun;

	_timer(&hr,&min,&sec,&hun);
	return (( (long) hr*60+min)*60+sec)*100+hun;
}

_timer()
{
#asm
	mov	ah,2ch
	int	21h
	mov	bx,[bp+4]
	mov	[bx],ch		;hours
	mov	byte [bx+1],0
	mov	bx,[bp+6]
	mov	[bx],cl		;minutes
	mov	byte [bx+1],0
	mov	bx,[bp+8]
	mov	[bx],dh		;seconds
	mov	byte [bx+1],0
	mov	bx,[bp+10]
	mov	[bx],dl		;hundreths
	mov	byte [bx+1],0
#endasm
}

#endif

show_m(s,a) char *s; double *a;
{	int i,j;
	printf("%s\n",s);
	for (i=0; i<4; i++)
		{for (j=0; j<4; j++) printf("%8.2f ",*a++);
		printf("\n");
		}
}