#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <ctype.h>

#define MAXLEN 8192

int verbose = 1;

/*----------------------------------------------------------------------*\
 | Computes PHI* shape function of Zahn and Roskies (1972) IEEE			|
 | Trans. C-21(3):269-281.  PHI* represents a closed plane curve "as	|
 | a function of arc length by the accumulated change in direction of	|
 | the curve since the starting point".  Phi is the net amount of		|
 | angular bend around its perimeter.  PHI* is normalized by removing	|
 | a circle from PHI, such that PHI*(0) and PHI*(2PI) = 0.  Thus,		|
 | "PHI* measures the way in which the shape in question differs from	|
 | a circular shape".													|
 |																		|
 | X and Y are the integer X,Y coordinates of NPT edge-points.			|
 | PHI* is the normalized net angular change as a function of ARC		|
 |																		|
 | G.P. Lohmann (WHOI, Woods Hole, MA 02543)  18 Oct 81					|
 | P.N. Schweitzer (USGS, Reston VA 22092)	 4 Nov 91 to C				|
 |																		|
 | After a program written in PL/1 by Zahn and Roskies.					|
\*----------------------------------------------------------------------*/

void phistar (double x[], double y[], double arc[], double phi[], int npt) {
	int i,j;
	double dx1,dx2,dy1,dy2,c,s,ssq,t,circle,delta;
	double p;

	/* If last point is not the same as the first, make it so. */

	if (x[npt-1] != x[0] || y[npt-1] != y[0]) {
		x[npt] = x[0];
		y[npt] = y[0];
		npt++;
		}

	/* Calculate arclength function */

	arc[0] = 0;
	phi[0] = 0;
	p = 0;
	for (i=1; i < npt; i++) {
		p += sqrt ((x[i] - x[i-1]) * (x[i] - x[i-1]) +
				   (y[i] - y[i-1]) * (y[i] - y[i-1]));
		arc[i] = p;
		}

	/* Store second point after the last (the last is first already */

	x[npt] = x[1];
	y[npt] = y[1];

	for (i=1; i < npt; i++) {
		dy1 = y[i] - y[i-1];
		dx1 = x[i] - x[i-1];
		dy2 = y[i+1] - y[i];
		dx2 = x[i+1] - x[i];
		ssq = (dx1*dx1 + dy1*dy1) * (dx2*dx2 + dy2*dy2);
		if (ssq == 0) {
			phi[i] = 0;
			x[i] = x[i-1];
			y[i] = y[i-1];
			}
		else {
			c = (dx1*dx2 + dy1*dy2)/sqrt(ssq);
			s = 0;
			if (c*c < 1) s = sqrt (1-c*c);

			/*----------------------------------------------------------*\
			 | Turbo C has a bug in its atan2 routine; it incorrectly	|
			 | considers -0 to be a negative number, which affects the	|
			 | value of phi if the outline takes a right angle.			|
			\*----------------------------------------------------------*/

			phi[i] = atan2 (s,c);
			if (dx1*dy2 - dx2*dy1 < 0) phi[i] *= -1;
			}
		}

	t = 0;
	for (i=0; i < npt; i++) {
		t += phi[i];
		phi[i] = t;
		}

	circle = phi[npt-1];
	if (verbose) {
		if (circle > 0) puts ("Outline was drawn counterclockwise");
		else			puts ("Outline was drawn clockwise");
		}

	delta = phi[npt-1]/arc[npt-1];
	for (i=0; i < npt; i++) phi[i] -= arc[i] * delta;

	if (circle > 0)
		for (i=0; i < npt/2; i++) {
			j = npt - i;
			t	  = phi[i];
			phi[i] = phi[j];
			phi[j] = t;
			}
	else {
		for (i=npt-2; i > 0; i--) phi[i+1] = phi[i];
		phi[0] = phi[npt-1];
		}
	npt--;
	}

/*----------------------------------------------------------------------*\
\*----------------------------------------------------------------------*/

main (int argc, char *argv[]) {
	int i,j,k,n,npt,idnum;
	FILE *in;
	char line [MAXLEN],*s;
	double x [MAXLEN];
	double y [MAXLEN];

	double arc [MAXLEN];
	double phi [MAXLEN];

	if (argc > 1) {
		if (in = fopen (argv[1],"r")) {

			/*----------------------------------------------------------*\
			\*----------------------------------------------------------*/

			idnum = 0;
			while (fgets (line,MAXLEN,in)) {
				if (s = strrchr (line,'\n')) *s = 0;

				/*--------------------------------------------------*\
				 | Specialized code to read BNA files				|
				\*--------------------------------------------------*/

				idnum++;
				for (s=line+strlen(line)-1; s > line && isdigit(*s); s--);
				npt = atoi (s);
				if (npt > 0) {
					if (verbose) printf ("shape %d has %d points\n",idnum,npt);
					i = 0;
					while (i < npt)
						if (fgets (line,MAXLEN,in))
							if (sscanf (line,"%lf%lf",&x[i],&y[i]) != 2) {
								puts ("Warning: expected <x> <y> in line below:");
								puts (line);
								}
							else {
								/* printf ("%4d: %12.6lf %12.6lf\n",i,x[i],y[i]); */
								i++;
								}
						else {
							puts ("Error: input file is shorter than expected");
							exit (1);
							}

					/*----------------------------------------------*\
					 | If the last point is not the same as the		|
					 | first, make it so.							|
					\*----------------------------------------------*/

					if (x[npt-1] != x[0] || y[npt-1] != y[0]) {
						if (verbose) printf ("(closing shape)\n");
						x[npt] = x[0];
						y[npt] = y[0];
						npt++;
						}

					/*----------------------------------------------*\
					 | Calculate phi-star function					|
					\*----------------------------------------------*/

					phistar (x,y,arc,phi,npt);

					}
				else {
					puts ("Warning: Expected <idnum> <npts> in line below:");
					puts (line);
					}
				}

			fclose (in);
			}
		else
			puts ("Error: Input file could not be opened");
		}
	else
		puts ("Usage: name input file on command line");
	}