/************************************************************************/
/*                                                                   	*/
/*	COPYRIGHT NOTICE:  This software was written and/or developed   */
/*      by John A. Hipp, Ph.D. at the Orthopaedic Biomechanics Lab      */
/*      of Beth Israel Hospital and Harvard Medical School.             */
/*                                                                      */
/*									*/
/*	Program:  t3m							*/
/*	Date:	May 2, 1995						*/
/*	Written by:	John A. Hipp, Ph.D.				*/
/*			Alan O. Jansujwicz, M.S.			*/
/*			Craig A. Simmons, M.S.				*/
/*									*/
/* 	Software to perform three-dimensional stereology.  The methods 	*/
/*	are based on a three-dimensional version of the directed secant	*/
/*	algorithm (Saltykov, 1958).  Standard morphology parameters are	*/
/*	calculated, based on the parallel plate model (Parfitt et al.,	*/
/*	JCI 72:1396-1409, 1983) and are described using the recommended	*/
/*	nomenclature (Parfitt et el., JBMR, 2:595-610, 1987).  3-D mean	*/
/*	intercept length vectors, MIL(theta,phi), are also calculated 	*/
/*	and fit to an ellipsoid from which the three principal MIL 	*/
/*	vectors are determined.						*/
/*									*/
/*	The program expects the following inputs:			*/
/*		1)  The image array is input through stdin, and is 	*/
/*			expected to be in byte form, where the 		*/
/*			background is represented by high voxel values	*/
/*			and the bone is	represented by voxel values 	*/
/*			less than the threshold. The voxels are assumed	*/
/*			to be cubiodal (equal side dimensions).  The 	*/
/*			format of the byte array is assumed to be 	*/
/*			(x, y, z).					*/
/*		2)  An ascii file called t3m.dat which contains the 	*/
/*			following parameters (in the this format):	*/
/*									*/
/*			xsize	ysize	zsize   (of image array)	*/
/*			xctr 	yctr 	zctr	(center coords of 	*/
/*							test sphere)	*/
/*			radius			(test sphere radius in 	*/
/*								voxels)	*/
/*			pixsize			(size in microns of one	*/
/*				    side of a voxel (an integer value))	*/
/*		3)  An ascii file called thresh.dat which contains the	*/
/*			threshold value which separates bone from 	*/
/*			background					*/
/*									*/
/*	For example, a typical t3m.dat file would be:			*/
/*									*/
/*			128 128 128					*/
/*			64 64 64					*/
/*			50						*/
/*			90						*/
/*									*/
/*	And a typical thresh.dat file would be:				*/
/*									*/
/*			100						*/
/*									*/
/*	For these examples, the program would read in a 128x128x128 	*/
/*	array with a voxel size of 90x90x90 microns cubed.  The 	*/
/*	analysis would be confined to a sphere with a 50 voxel radius, 	*/
/*	centered at the (64,64,64) position in the array.  Bone voxels 	*/
/*	would be defined as any voxel with a value below the threshold 	*/
/*	of 100.								*/
/*									*/
/*	The program produces the following outputs:			*/
/*		1)  Output to stderr containing the stereology analysis	*/
/*		2)  An ascii file (named through stdout) which contains */
/*			the stereology results in a space-delimited	*/
/*			format so that it can be read into various 	*/
/*			external applications. The order of the results	*/
/*			can be determined from the code.		*/
/*		3)  An ascii file called mil_3d.dat which contains a 	*/
/*			list of	the test grid rotations, and for each 	*/
/*			rotation, the number of intersections per unit 	*/
/*			test line length, the MIL and the number of	*/
/*			intersections.					*/
/*									*/
/*	We compile the program on a SUN using the command:		*/
/*									*/
/*	cc  t3m.c -o t3m -lm -O3 -Bdynamic 	*/
/*									*/
/*      You must have the Math++.h file in the include file search path */
/*									*/
/*	To run the program use:						*/
/*									*/
/*	t3m < array.byt >! data.out					*/
/*									*/
/************************************************************************/

#include <stdio.h>
#include <math.h>
#include <time.h>
#include <Math++.h>

#define EXIT_SUCCESS    0
#define EXIT_FAILURE    -1

#define BONE	0
#define BACKGROUND	255

#define MAXPTS 128  /*  The maximum # of MIL rotations	*/

#define LINE_POS sqrt((x - xctr) * (x - xctr) + (y - yctr) * (y - yctr) + (z - zctr) * (z - zctr))  /* Definition simply for convenience */

/*  Variables required to generate random numbers  */
static long state1[32] = {
	3,
	0x9a319039, 0x32d9c024, 0x9b663182, 0x5da1f342,
	0x7449e56b, 0xbeb1dbb0, 0xab5c5918, 0x946554fd,
	0x8c2e680f, 0xeb3d799f, 0xb11ee0b7, 0x2d436b86,
	0xda672e2a, 0x1588ca88, 0xe369735d, 0x904f35f7,
	0xd7158fd6, 0x6fa6f051, 0x616e6b96, 0xac94efdc,
	0xde3b81e0, 0xdf0a6fb5, 0xf103bc02, 0x48f340fb,
	0x36413f93, 0xc622c298, 0xf5a42ab8, 0x8a88d77b,
	0xf5ad9d0e, 0x8999220b, 0x27fb47b9
	};

long random();
char *setstate();
char *initstate();

extern char *calloc();

time_t  begin_time, end_time;
double	elapsed_time_min;

int	xsize, ysize, zsize;	/* size of array			*/
int 	xctr, yctr, zctr;	/* center of test sphere		*/
int	radius, thresh;		/* radius of test sphere and threshold	*/
double 	pixsize;		/* size of one side of voxel		*/

double	tline_sp;		/* spacing between test lines in mm	*/
double	tline_ln;		/* length of all test lines in mm	*/
int	num_tlines;		/* number of test lines in grid		*/
long	tintrsctn;		/* total # of intersections		*/
double 	rot1[MAXPTS];		/* random theta rotations		*/
double	rot2[MAXPTS];		/* random phi rotations			*/

double 	solidvol, bvtv;		/* solid volume in mm^3 and BV/TV	*/
double 	pl;			/* #intrsctns/unit test line length	*/
double	mil[MAXPTS];		/* magnitude of MIL for ith rotation	*/
double	sv, tb, tpd, tps;	/* BS/TV, Tb.Th, Tb.N and Tb.Sp		*/

double  ae,be,ce,de,ee,fe;	/* ellipsoid coefficients 		*/
double  corr,sqres;		/* ellipsoid correlation statistics	*/
double	stdev,variance;		/* ellipsoid variance statistics	*/

double  prin1,prin2,prin3;	/* principal MIL vectors		*/
double  orient1a,orient1b;	/* orientation of 1st principal		*/
double	orient2a,orient2b;	/* orientation of 2nd principal		*/
double	orient3a,orient3b;	/* orientation of 3rd principal		*/
double  deg1,deg2,deg3;		/* degrees of anisotropy		*/
double 	tval;			/* test statistic for MIL indpendence	*/
double	*ci,*test,*princi;	/* confidence interval and test stat	*/

FILE 	*fpi, *fpo, *fpt, *fpm;	/* file pointers			*/

long 	euler;			/* Euler connectivity value		*/


main() 
{
unsigned char *image;		/* image array				*/
unsigned tsize;			/* total size of image (i.e. XxYxZ)	*/

double	*eigvals,*eigvecs;	/* e-values and e-vectors for MIL tensor*/
double	sdev, adev, mean, minmil;/* statistical variables for MILs	*/
double	curt, skew, svar, s, p;	/* more stats for MILs			*/

long 	got;			/* temporary long variable		*/
double 	dtemp, dt;		/* temporary variables for calculations	*/
int 	i;			/* loop counter				*/

/*  The image file and the stereologic analysis file are defined through */
/*	directed input and output, respectively. 			 */

	fpi = stdin;
	fpo = stdout;

	eigvals    = (double *)calloc(9,sizeof(double));
	eigvecs    = (double *)calloc(9,sizeof(double));
	ci	   = (double *)calloc(3,sizeof(double));
	test	   = (double *)calloc(3,sizeof(double));
	princi     = (double *)calloc(6,sizeof(double));

/*  Parameter file t3m.dat contains the array size, the location and 	*/
/*	radius of the test sphere (in pixels) and the size of a pixel 	*/
/*	in microns  							*/

	if((fpt = fopen("t3m.dat","r")) == NULL) {
		fprintf(stderr,
		  "ERROR: could not find the t3m.dat parameter file\n\n");
		exit(EXIT_FAILURE);
	}

	fscanf(fpt,"%d %d %d", &xsize,&ysize,&zsize);
	fscanf(fpt,"%d %d %d", &xctr,&yctr,&zctr);
	fscanf(fpt,"%d", &radius);
	fscanf(fpt,"%d", &got);

	fclose(fpt);

	pixsize = (double)got / 1000.0;	/* convert the pixel size from 
						micrometers to mm */
 

/*   Parameter file thresh.dat contains the threshold value to 		*/
/*	differentiate bone from marrow.  A value below the threshold is */
/*	defined as bone 						*/

	if((fpt = fopen("thresh.dat","r")) == NULL) {
		fprintf(stderr,
		  "ERROR: could not find the thresh.dat parameter file\n\n");
		exit(EXIT_FAILURE);
	}

	fscanf(fpt,"%d", &thresh);

	fclose(fpt);

/*   Allocate memory for image array */
	
	tsize = xsize*ysize*zsize;
	image = (unsigned char *)calloc(tsize,sizeof(unsigned char));
	if(image == NULL) {
		fprintf(stderr,"ERROR: could not allocate memory for image\n");
		exit(EXIT_FAILURE);
	}

	time(&begin_time);  /* begin clock */

/*   Read in image array */

	got = fread(image, sizeof(unsigned char),tsize,fpi);
	if(got < tsize) {
		fprintf(stderr,"ERROR: only got %ld of %ld voxels\n",got,tsize);
		exit(EXIT_FAILURE);
	}
	fclose(fpi);

/*   Call routine to determine stereology parameters, based on random 	*/
/*	rotations of the test line grid 				*/

	got = gridrota(image);

/*   Determine descriptive statistics for MIL vectors  */

	mean = 0.0;
	for(i=0;i<MAXPTS;i++) mean += mil[i];
	mean /= (double)MAXPTS;				/* Mean MIL */
	adev = svar = skew = curt = 0.0;
	for(i=0;i<MAXPTS;i++) {
		adev += fabs(s=(mil[i] - mean));
		svar += (p = s*s);
		skew += (p *= s);
		curt += (p *= s);
	}
	adev /= (double)MAXPTS;
	svar /= ((double)MAXPTS - 1.000);
	sdev = sqrt(svar);				/* Std deviation */
	skew /= (MAXPTS * svar * sdev);			/* Skew */
	curt = curt / (MAXPTS * svar * svar) - 3.0;	/* Kurtosis */

	minmil = mil[0];
	for(i=1;i<MAXPTS;i++) 
		if(mil[i] < minmil) minmil = mil[i];	/* Minimum MIL */

/*   Call routine to determine the best fit ellipsoid equation for data	*/

	mil2_fit();

/*   Call routine to calculate the properties of the ellipsoid tensor  */

	anisotropy(eigvals,eigvecs);

/*   Call routine to determine the independence of the eigenvalues  */

	eigstats(eigvals,eigvecs);

/*   Call routine to sort the data  */

	sortdata(eigvals);

/*   Print the results to stderr  */ 

	fprintf(stderr, "\nSTEREOLOGY RESULTS\n");
	fprintf(stderr, "==================\n");
	fprintf(stderr, "ANALYSIS SPECIFICATIONS\n");
	fprintf(stderr, "-----------------------\n");
	fprintf(stderr, "Total grid size: %d x %d x %d\t", xsize, ysize, zsize);
	fprintf(stderr, "Sphere centered at: (%d, %d, %d)\n", xctr, yctr, zctr);
	fprintf(stderr, "Sphere radius (pixels): %d\t\t", radius);
	fprintf(stderr, "Threshold value: %d\n", thresh);
	fprintf(stderr, "Pixel size (mm^3/voxel): %6.4lf x %6.4lf x %6.4lf\n",
						pixsize, pixsize, pixsize);
	dt = radius * pixsize;
	dtemp = 4.000000 / 3.00000000 * M_PI * dt * dt * dt;
	fprintf(stderr, "VOLUME FRACTION\n");
	fprintf(stderr, "---------------\n");
	fprintf(stderr, "BV/TV: %.4lf\n",bvtv);
	fprintf(stderr, "Total volume(mm^3): %.2lf\t",dtemp);
	fprintf(stderr, "Solid volume(mm^3): %.2lf\n",solidvol);
	fprintf(stderr, "INTERSECTIONS\n");
	fprintf(stderr, "-------------\n");
	fprintf(stderr, "Test line spacing (mm): %.3lf\t\t",tline_sp);
	fprintf(stderr, "Number of test lines: %d\n",num_tlines);
	fprintf(stderr, "Avg tline length (mm): %.2lf\t\t",tline_ln);
	fprintf(stderr, "Total # intrsctns: %ld\n",tintrsctn);
	fprintf(stderr, "# of intrsctns/tline length(mm^-1): %.4lf\n", pl);
	fprintf(stderr, "MORPHOLOGY PARAMETERS\n");
	fprintf(stderr, "---------------------\n");
	fprintf(stderr, "BS/BV (mm^-1): %.4lf\t",sv);
	fprintf(stderr, "Tb.Th (mm): %.4lf\n",tb);
	fprintf(stderr, "Tb.N (mm^-1): %.4lf\t",tpd);
	fprintf(stderr, "Tb.Sp (mm): %.4lf\n",tps);
	fprintf(stderr, "ANISOTROPY PARAMETERS\n");
	fprintf(stderr, "---------------------\n");
	fprintf(stderr, "Principal MILs (mm):\t%.4lf\t%.4lf\t%.4lf\n",
							prin1, prin2, prin3); 
	fprintf(stderr, "\nMean MIL (mm): %.4lf\t", mean);
	fprintf(stderr, "Stdev of MILs (mm): %.4lf\n",sdev);
	fprintf(stderr, "Kurtosis: %.4lf\t",curt);
	fprintf(stderr, "Skew of MILs: %.4lf\n",skew);
	fprintf(stderr, "\nOrientation 1 (theta,phi): (%.2lf, %.2lf)\n",
						orient1a, orient1b);
	fprintf(stderr, "Orientation 2 (theta,phi): (%.2lf, %.2lf)\n",
						orient2a, orient2b);
	fprintf(stderr, "Orientation 3 (theta,phi): (%.2lf, %.2lf)\n",
						orient3a, orient3b);
	fprintf(stderr, "\nDegree of anisotropy MIL 1: %.4lf\n", deg1);
	fprintf(stderr, "Degree of anisotropy MIL 2: %.4lf\n", deg2);
	fprintf(stderr, "Degree of anisotropy MIL 3: %.4lf\n", deg3);
	fprintf(stderr, "95%% Confidence intervals for principals:\n");
	fprintf(stderr, "\t%6.4lf < MIL 1 < %6.4lf \n",*(princi+0),*(princi+1));
	fprintf(stderr, "\t%6.4lf < MIL 2 < %6.4lf \n",*(princi+2),*(princi+3));
	fprintf(stderr, "\t%6.4lf < MIL 3 < %6.4lf \n",*(princi+4),*(princi+5));
	fprintf(stderr, "ELLIPSOID FIT PARAMETERS\n");
	fprintf(stderr, "------------------------\n");
	fprintf(stderr, "Ellipsoid coefficients:\n");
	fprintf(stderr, "\tae = %.4lf\tbe = %.4lf\tce = %.4lf\n", ae, be, ce);
	fprintf(stderr, "\tde = %.4lf\tee = %.4lf\tfe = %.4lf\n", de, ee, fe);
	fprintf(stderr, "\nSum of squares: %.4lf\t", sqres);
	fprintf(stderr, "Correlation coeff. (r): %.4lf\n", corr);
	fprintf(stderr, "EIGENVALUE/VECTOR PARAMETERS\n");
	fprintf(stderr, "----------------------------\n");
	fprintf(stderr, "Eigenvalues:");
	fprintf(stderr, "\t%.4lf\t%.4lf\t%.4lf\n", 
				*(eigvals+0), *(eigvals+4), *(eigvals+8));
	fprintf(stderr, "Eigenvectors:\n");
	fprintf(stderr, "\tE-vector 1: (%.4lf,%.4lf,%.4lf)\n", *(eigvecs+0),
		 			*(eigvecs+3),*(eigvecs+6));
	fprintf(stderr, "\tE-vector 2: (%.4lf,%.4lf,%.4lf)\n", *(eigvecs+1), 
					*(eigvecs+4),*(eigvecs+7));
	fprintf(stderr, "\tE-vector 3: (%.4lf,%.4lf,%.4lf)\n", *(eigvecs+2),
					*(eigvecs+5),*(eigvecs+8));
	fprintf(stderr,"Eigenvalue Independence:\n");
	fprintf(stderr,"\ttvalue   teststat1   teststat2   teststat3\n");
	fprintf(stderr,"\t------   ---------   ---------   ---------\n");
	fprintf(stderr,"\t%6.4lf     %6.4lf      %6.4lf      %6.4lf \n",
					tval,*(test+0),*(test+1),*(test+2));

/*   Write the a subset of the results to a space delimited file.  The 	*/
/*	order of the columns can be determined from the comments in the	*/
/*	following code.							*/

	fprintf(fpo,"%6.4lf ",pixsize);		/* pixel size in mm */
	fprintf(fpo,"%8.4lf ",solidvol);	/* # pixels solid volume */
	fprintf(fpo,"%8.4lf ",bvtv);		/* BV/TV */
	fprintf(fpo,"%8.4lf ",tline_sp);	/* test line spacing */
	fprintf(fpo,"%8.4lf ",tline_ln);	/* total length of test line */
	fprintf(fpo,"%ld ",tintrsctn);		/* total # of intersections */
	fprintf(fpo,"%8.4lf ",pl);		/* #intrsct/test line length */
	fprintf(fpo,"%8.4lf ",sv);		/* BS/BV */
	fprintf(fpo,"%8.4lf ",tb);		/* Tb.Th */
	fprintf(fpo,"%8.4lf ",tpd); 		/* Tb.N */
	fprintf(fpo,"%8.4lf ",tps); 		/* Tb.Sp */
	fprintf(fpo,"%lf %lf %lf ", prin1, prin2, prin3);      /* Princ MILs */
	fprintf(fpo,"%8.4lf ",mean);   		/* mean MIL */
	fprintf(fpo,"%8.4lf ",sdev);  		/* stdev of MIL */
	fprintf(fpo,"%8.4lf ",minmil); 		/* minimum MIL */
	fprintf(fpo,"%lf %lf ", orient1a, orient1b);	    /* Princ1 angles */
	fprintf(fpo,"%lf  %lf ", orient2a, orient2b);	    /* Princ2 angles */
	fprintf(fpo,"%lf %lf ", orient3a, orient3b);	    /* Princ3 angles */	
	fprintf(fpo,"%lf %lf %lf ", deg1, deg2, deg3);	 /* degrees of aniso */
	fprintf(fpo,"%6.4lf ", *(princi+0));	      /* Lower CI for Princ1 */
	fprintf(fpo,"%6.4lf ", *(princi+1));	      /* Upper CI for Princ1 */
	fprintf(fpo,"%6.4lf ", *(princi+2));	      /* Lower CI for Princ2 */
	fprintf(fpo,"%6.4lf ", *(princi+3));	      /* Upper CI for Princ2 */
	fprintf(fpo,"%6.4lf ", *(princi+4));	      /* Lower CI for Princ3 */
	fprintf(fpo,"%6.4lf ", *(princi+5));	      /* Upper CI for Princ3 */
	fprintf(fpo,"%lf  %lf  %lf  %lf  %lf  %lf ", ae,be,ce,de,ee,fe);
						   /* ellipsoid coefficients */
	fprintf(fpo,"%lf %lf ", sqres, corr);	       /* ellipsoid SS and r */
	fprintf(fpo,"%6.4lf ", *(test+0) - tval);      /* stat for eval1 ind */
	fprintf(fpo,"%6.4lf ", *(test+1) - tval);      /* stat for eval2 ind */
	fprintf(fpo,"%6.4lf ", *(test+2) - tval);      /* stat for eval3 ind */
	fclose(fpo);

	time(&end_time);	/* stop the clock */
	elapsed_time_min = (end_time -  begin_time) / 60.;

	fprintf(stderr,"\nAll DONE!!! Elapsed time = %3.2f min\n",
							elapsed_time_min);
	exit(0);
}


/************************************************************************/
/*   Function:  gridrota (array)					*/
/*									*/
/*	Function to scan the image array using a three-dimensional 	*/
/*	version of the directed secant method.  The image is scanned by */
/*	a 3-D test grid at randomly determined orientations.  The 	*/
/*	orientations are defined by spherical angles (theta, phi).  	*/
/*	Theta is the rotation about the x-axis and phi is the rotation 	*/
/*	about the z-axis.  Based on the threshold value, the image is 	*/
/*	converted to a binary format and the bone volume fraction is 	*/
/*	determined.  The image is then systematically scanned and 	*/
/*	intersections are recorded when the binary value of the current	*/
/*	voxel differs from the binary value of the previous voxel.	*/
/*	Standard morphology parameters are calculated based on the 	*/
/*	parallel plate model:						*/
/*		BV/TV = # bone voxels/total # voxels in test sphere	*/
/*		Tb.N = total # intersections/total length of test lines	*/
/*		Tb.Th = BV/TV / Tb.N					*/
/*		Tb.Sp = (1-BV/TV) / Tb.N				*/
/*		BS/BV = 2 * Tb.N / BV/TV				*/
/*	Three-dimensional mean intercept length vectors are calculated 	*/
/*	for each rotation as:						*/
/*		MIL(theta, phi) = (BV/TV) / ((1/2) * #intersections	*/
/*					per unit test line length)	*/
/*									*/
/************************************************************************/

gridrota(array)
unsigned char *array;
{

double	maxrand;			/* value of maximum random num.	*/
unsigned char seed;			/* seed for random number gen.	*/
int	n;				/* for random number gen.	*/

double	dradius;			/* test sphere radius (double)	*/
int	lmarg, rmarg;			/* x limits of analysis region	*/
int	tlay, blay;			/* y limits of analysis region	*/
int 	top, bottom;			/* z limits of analysis region	*/
double 	tvol;				/* total number of voxels	*/

double	theta, phi;			/* orientation angles		*/
double	sin_theta, cos_theta;		/* trig calculations		*/
double	sin_phi, cos_phi;		/* more trig calculations	*/
double	ctcp, ctsp;			/* and more trig calculations	*/
double 	r[3], norm_r[3];		/* test line direction vectors	*/
double 	x, y, z;			/* position on test line	*/
double 	x_off, y_off, z_off;		/* offset to get to start pos.	*/
double 	dist;				/* dist. from ctr to grid pt.	*/
double 	dx, dy, dz;			/* grid position (double)	*/
int	ix, iy, iz;			/* integer grid position	*/
int	iline_sp;			/* spacing of test lines (int)	*/
long 	intrsctn;			/* number of intersections	*/
long 	cur_pt;				/* current location in array	*/
unsigned char preval;			/* previous voxel value		*/

int	irot;				/* rotation loop counter	*/
int	i;				/* loop counter			*/

	fprintf (stderr, "\nStarting analysis\n");
	
	dradius = (double)radius;

	maxrand = 2147483647;   /* The max value generated by the random*/
				/* number routine is (2**31 - 1) 	*/

/*   Generate the random rotations */

	seed = 1;
	n = 128;
	initstate(seed, (char *) state1, n);
	setstate(state1);

	for(irot=0;irot < MAXPTS;irot++) {
		rot1[irot] =  random() * M_PI / (double)maxrand;
		rot2[irot] =  random() * M_PI / (double)maxrand;
	}
	
/*  Calculate the test line spacing in pixels in x-y plane using a 200	*/ 
/*	micron spacing							*/

	iline_sp = 0.2000 / pixsize;  

	tline_sp = (double)iline_sp * pixsize; 

/*  Open the file for that will contain the MIL calculated at each rotation */

	if((fpm = fopen("mil_3d.dat","w")) == NULL) {
		fprintf(stderr,"ERROR: could not write output to mil_3d.dat\n");
		return(1);
	}

	tintrsctn = 0;
	num_tlines = 0;

/*  Determine the limits of the test region, in case the requested sphere is	*/
/*	not within the array dimensions.					*/

	lmarg = xctr - radius;
	if (lmarg < 0) lmarg = 0;
	rmarg = xctr + radius;
	if (rmarg > xsize) lmarg = xsize;

	tlay = yctr - radius;
	if (tlay < 0) tlay = 0;
	blay = yctr + radius;
	if (blay > ysize) blay = ysize;

	top = zctr - radius;
	if (top < 0) top = 0;
	bottom = zctr + radius;
	if (bottom > zsize) bottom = zsize;

/*  Calculate the total test line length (same for all rotations).  	*/
/*	For any (x,y) point on the test grid, the test line length is 	*/
/*	twice the perpendicular distance from the x-y plane to the edge */
/*	of the test sphere.		  				*/	

	tline_ln = 0.0000000000000;
	for(ix = lmarg; ix <= rmarg; ix += iline_sp) {
		dx = (double)ix - (double)xctr;
		for(iy = tlay; iy <= blay; iy += iline_sp) {
			dy = (double)iy - (double)yctr;
			dist = sqrt(dx * dx + dy * dy);
			if(dist <= dradius) {
				dz = sqrt(dradius * dradius - dist*dist);
				tline_ln += 2.0000 * dz;
			}
		}
	}

	tline_ln *= pixsize;  /* Convert from pixels -> mm */

/*  Calculate the solid volume (#voxels) of the test region and		*/
/*	threshold the image to a bone/background representation. 	*/

	solidvol = tvol = 0.00000000000000;

	for(iz = top; iz < bottom; iz++) {
		for (iy = tlay; iy < blay; iy++) {
			for (ix = lmarg; ix < rmarg; ix++) {
				dx = (double)ix - (double)xctr;
				dy = (double)iy - (double)yctr;
				dz = (double)iz - (double)zctr;

				cur_pt = (iz * ysize + iy) * xsize + ix;

			/* Threshold image to BONE or BACKGROUND  */

				if (*(array + cur_pt) < thresh)
					*(array + cur_pt) = BONE;
				else *(array + cur_pt) = BACKGROUND;

			/* Determine # of BONE pixels */

				dist = sqrt(dx * dx + dy * dy + dz * dz);
				if(dist <= radius) {
					tvol++;
					if (*(array + cur_pt) == BONE) 
								solidvol++;
					else *(array + cur_pt) = BACKGROUND;
				}
			}
		}
	}

	bvtv = solidvol / tvol;	 /* Volume fraction, BV/TV */

	solidvol *= pixsize * pixsize * pixsize;  /* Scale voxels -> mm^3 */

	fprintf (stderr, "\nThe volume fraction is: %lf\n", bvtv);

	fprintf (stderr, "\nRotation counter (total of %d):\n", MAXPTS);
				
/*   Scan the image for each of the randomly selected rotations  */
	
	for( irot = 0; irot < MAXPTS; irot++) {   	
 
		theta = rot1[irot];
		phi = rot2[irot];

		sin_phi = sin(phi);
		cos_phi = cos(phi);
		sin_theta = sin(theta);
		cos_theta = cos(theta);
		ctcp = cos_theta*cos_phi;
		ctsp = cos_theta*sin_phi;

/*   Define direction vector for this rotation 	*/

		r[0] = cos_phi*sin_theta;
		r[1] = sin_phi*sin_theta;
		r[2] = cos_theta;

/*   Normalize direction vector to +1 for the maximum component */

	if ((fabs(r[0]) >= fabs(r[1]))&&(fabs(r[0]) >= fabs(r[2]))) {
		norm_r[0] = 1.0;
		norm_r[1] = r[1] / r[0];
		norm_r[2] = r[2] / r[0];
	}
	else if ((fabs(r[1]) >= fabs(r[0]))&&(fabs(r[1]) >= fabs(r[2]))) {
		norm_r[1] = 1.0;
		norm_r[0] = r[0] / r[1];
		norm_r[2] = r[2] / r[1];
	}
	else {
		norm_r[2] = 1.0;
		norm_r[0] = r[0] / r[2];
		norm_r[1] = r[1] / r[2];
	}

/*   Make sure signs are correct */

	for (i=0; i<3; i++) {
		if (r[i] > 0) norm_r[i] = fabs(norm_r[i]);
		else norm_r[i] = (-fabs(norm_r[i]));
	}

	intrsctn = 0;

/*  Loop over the test grid which is equally spaced in x and y */

	for(ix = lmarg; ix <= rmarg; ix += iline_sp) {
		for(iy = tlay; iy <= blay; iy += iline_sp) {

			dx = (double)ix - (double)xctr;
			dy = (double)iy - (double)yctr;
			dist = sqrt(dx * dx + dy * dy);
	
			if(dist < dradius) {
				dz = sqrt(dradius * dradius - dist * dist);
				x_off = dz * r[0];
				y_off = dz * r[1];
				z_off = dz * r[2];

				x = dx*ctcp - dy*sin_phi + (double)xctr;
				y = dx*ctsp + dy*cos_phi + (double)yctr;
				z = (-dx)*sin_theta + (double)zctr;

			/* Starting coordinate for this rotation */
				x -= x_off;
				y -= y_off;
				z -= z_off;

			/* First determine whether the line starts in 	*/
			/* bone or background		 		*/
				cur_pt = (anint(z)*ysize + anint(y))*xsize +
								 anint(x);
				preval = *(array + cur_pt);

			/* Increment current position along the test line */
				x += norm_r[0];
				y += norm_r[1];
				z += norm_r[2];

				while (LINE_POS <= dradius) {
					cur_pt = (anint(z) * ysize + 
						anint(y)) * xsize + anint(x);
			
				/* Count an intersection if the current	*/
				/* voxel value differs from the previous*/
				/* value.				*/
					if(*(array + cur_pt) != preval) {
						intrsctn += 1;
						preval = *(array + cur_pt);
						}

					x += norm_r[0];
					y += norm_r[1];
					z += norm_r[2];
				} /* end while on line loop */

				num_tlines++;

				} /* end loop over z axis */
	
			} /* end loop over y axis */

		} /* end loop over x axis */
						
/* Calculate the number of intersections per test line length and the 	*/
/*	mean intercept length (MIL) for this rotation.			*/

	pl = (double)intrsctn / tline_ln;
	mil[irot] = (2.0 * bvtv) / pl; 
	tintrsctn += intrsctn;
	if ((irot % 10) == 0) fprintf(stderr, "%d", irot);
	else fprintf(stderr,".");

/* Write the angles, #intrsctns/test line length, MIL and #intrsctns 	*/
/*	to the file mil_3d.dat.						*/ 

	fprintf(fpm,"%lf %lf %lf %lf %ld\n",
		rot1[irot],rot2[irot], pl, mil[irot], intrsctn);

	}   /* end loop over MAXPTS rotations */


	fprintf(stderr,"\nDONE GRIDROTA\n\n");

/* Calculate morphology parameters: BS/BV, Tb.Th, Tb.N, Tb.Sp, using 	*/
/*	parallel plate model (Parfitt, et al. JCI 72:1396-1409, 1983)	*/

	pl = (double)tintrsctn / (tline_ln * MAXPTS);

	sv = 2 * pl / bvtv;	/* BS/BV, [mm^2 mm^-3] */
	tb = bvtv / pl;		/* Tb.Th, [mm] */
	tpd = pl;		/* Tb.N, [mm^-1] */
	tps = (1.0-bvtv)/pl;	/* Tb.Sp, [mm] */
	num_tlines /= MAXPTS;	/* # of test lines for a single rotation */
	
	fclose(fpm);

return(EXIT_SUCCESS);
}


/************************************************************************/
/* Function:  mil2_fit()						*/
/*									*/
/*	This function fits the 3-D MIL data to an ellipsoid, and 	*/
/*	determines the goodness of the fit.  These methods are based on */
/*	2-D methods of Whitehouse (JMicroscopy 101:153-168, 1974) and	*/
/*	Harrigan and Mann (JMatSci 19:761-767, 1984).  The approach is 	*/
/*	to plot the locus of the end points of the MIL vectors issuing 	*/
/*	from a common center and fit them to an ellipsoid of general 	*/
/*	formula:							*/
/*									*/
/*	A*n1^2 + B*n2^2 + C*n3^2 + D*n1*n2 + E*n1*n3 + F*n2*n3 = 1/L^2	*/
/*									*/
/*	where L is the length of the MIL, ni are the direction cosines	*/
/*	between L and the base vectors in an arbitrary coordinate 	*/
/*	system and A...F are the ellipsoid coefficients.		*/
/*									*/
/*	For this code, a multivariable linear least squares fitting 	*/
/*	technique is used to fit the data by solving the linear system 	*/
/*	A * x = b, where A contains the projection data for each 	*/
/*	rotation, x is a column vector of the ellipsoid coefficients 	*/
/*	and b is a column vector of 1/L^2, where L is the magnitude of 	*/
/*	the MIL vector.  The solution is formulated as:			*/
/*									*/
/*			x = (At * A)^-1 * At * b			*/
/*									*/	/*	The sum of the squares of the residuals and the correlation 	*/
/*	coefficient of the fit are also calculated.			*/
/*									*/
/************************************************************************/

mil2_fit() 
{

double	x, y, z;			/* direction vector projections	*/
double	a[MAXPTS][6];			/* A matrix 			*/
double	at[6][MAXPTS];			/* A transpose 	(At)		*/
double	ata[6][6];			/* A transpose * A  (AtA)	*/
double	atai[6][6];			/* inverse of AtA  (AtA)^-1	*/
double	c[6][MAXPTS];			/* (AtA)^-1 * At		*/

double	xvector[6];			/* vector of ellipsoid coeff.	*/

double	b[MAXPTS];			/* vector of 1/(MIL)^2		*/
double	bmn, bfitmn;			/* mean of measured and fit b	*/
double	bfit[MAXPTS];			/* b vector based on fit	*/
double	expdif, expfit, expsqr;		/* stats for measured data	*/
double	fitsqr, fitdif;			/* stats for fit data		*/

int	i,j,k;				/* loop variables 		*/
int	degree;				/* degree of vectors (6)	*/

	degree = 6;


/* Initialize the matrices */

	for (i=0;i<degree;i++) {
		for (j=0;j<degree;j++) {
			ata[i][j]  = 0.0;
			atai[i][j] = 0.0;
		}
	}

	for (i=0;i<degree;i++) {
		for (k=0;k<MAXPTS;k++) {
			c[i][k] = 0.0;
		}
	}

	for (i=0;i<degree;i++) {
		xvector[i] = 0.0;
	}


/* Determine A and b matrices of Ax = b */
		
	for (i=0;i<MAXPTS;i++) {

		x = cos(rot2[i])*sin(rot1[i]);
		y = sin(rot2[i])*sin(rot1[i]);
		z = cos(rot1[i]);

		a[i][0] = x * x;
		a[i][1] = y * y;
		a[i][2] = z * z;
		a[i][3] = x * y;
		a[i][4] = y * z;
		a[i][5] = x * z;

		b[i] = 1.0 / (*(mil+i) * *(mil+i));

	}

/* Determine At of AtAx = Atb */

	for (i=0;i<MAXPTS;i++) {
		at[0][i] = a[i][0];
		at[1][i] = a[i][1];
		at[2][i] = a[i][2];
		at[3][i] = a[i][3];
		at[4][i] = a[i][4];
		at[5][i] = a[i][5];
	}

/* Multiply At * A */

	for (i=0;i<degree;i++) {
		for (j=0;j<degree;j++) {
			for (k=0;k<MAXPTS;k++) {
				ata[i][j] += at[i][k]*a[k][j];
			}
		}
	}


/* Find inverse of AtA */

	linInverseMatrix (degree,degree,ata,atai);


/* Multiply (AtA)^-1 * At */

	for (i=0;i<degree;i++) {
		for (k=0;k<MAXPTS;k++) {
			for (j=0;j<degree;j++) {
				c[i][k] += atai[i][j]*at[j][k];
			}
		}
	}


/* Multiply (AtA)^-1 * At * b */

	for (i=0;i<degree;i++) {
		for (k=0;k<MAXPTS;k++) {
			xvector[i] += c[i][k]*b[k];
		}
	}


/* Determine the goodness of fit.  Begin by calculating the sum of the	*/
/* 	square of the residuals. 					*/

	sqres = 0.0;

	for(i=0;i<MAXPTS;i++) {
		bfit[i] = xvector[0]*a[i][0] + xvector[1]*a[i][1] + 
			  xvector[2]*a[i][2] + xvector[3]*a[i][3] + 
			  xvector[4]*a[i][4] + xvector[5]*a[i][5];

		sqres += pow((b[i] - bfit[i]), 2.0);
	}

/*  Find the mean of b[] and bfit[] */

	bmn = 0.0;
	bfitmn = 0.0;

	for(i=0; i<MAXPTS; i++) {
		bmn += b[i];
		bfitmn += bfit[i];
	}

	bmn /= (double)MAXPTS;
	bfitmn /= (double)MAXPTS;

/* Find the variance of the error in the ellipsoid fit */

	variance = sqres / (double)(MAXPTS - 6 - 1);

	stdev	 = sqrt(variance);

/*  Calculate the correlation coefficients */

	expfit=0.0;
	expsqr=0.0;
	fitsqr=0.0;

	for(i=0; i< MAXPTS; i++) {
		expdif = b[i] - bmn;
		fitdif = bfit[i] - bfitmn;
		expfit += expdif * fitdif;
		expsqr += expdif * expdif;
		fitsqr += fitdif * fitdif;
	}

	if(( expsqr < 0.0) || (fitsqr < 0.0) ) corr = 1.0;  /* no zero div*/
	else	corr = expfit / sqrt(expsqr*fitsqr);
	
/* Return ellipse coefficients through external variables */

	ae = xvector[0];
	be = xvector[1];
	ce = xvector[2];
	de = xvector[3];
	ee = xvector[4];
	fe = xvector[5];

return (EXIT_SUCCESS);
}


/************************************************************************/
/* Function: anisotropy(eigvals, eigvecs)				*/
/*									*/
/*	This function defines the MIL tensor based on the ellipsoid fit	*/
/*	and determines the principal eigenvectors and eigenvalues, 	*/
/*	which are used to define the principal MIL vector orientations 	*/
/*	and magnitudes, respectively.  The MIL tensor assumes material	*/
/*	orthotropy and is defined as (Harrigan and Mann, JMatSci 	*/
/*	19:761-767, 1984):						*/
/*									*/
/*				| ae	  de/2	  ee/2 |		*/
/*			    M = | de/2	  be	  fe/2 |		*/
/*				| ee/2	  fe/2	  ce   |		*/
/*		where ae, be, ..,fe are the ellipsoid coefficients	*/
/*									*/
/*	The orientations of the principal MILs are defined by spherical */
/*	angles,	a and b, where a and b correspond to theta and phi,	*/
/*	respectively.  The magnitudes of the principal MILs are 	*/
/*	calculated as the inverse of the square root of the absolute 	*/
/*	values of the eigenvalues.					*/
/*	The three degrees of anisotropy are defined as the relative 	*/
/*	differences between the three prinicpal MILs and the mean of 	*/
/*	the principals.							*/
/*									*/
/*	The eigenvalues and eigenvectors are returned as arguments.	*/
/*									*/
/************************************************************************/

anisotropy(eigvals,eigvecs)

double	*eigvals,		/* Used for MIL tensor, then evalues	*/
	*eigvecs;		/* eigenvectors 			*/

{
double	eps;			/* convergence criterion 		*/
double	prinmn;			/* mean principal value of MIL		*/

int	dim,			/* dimension of evalue and evector array*/
	order,			/* order of matrices 			*/
	numiters,		/* number of iterations required	*/
	maxiters,		/* maximum nuber of iterations allowed	*/
	status;			/* status returned by eigJacobi		*/

	dim		= 3;
	order		= 3;
	numiters	= 0;
	maxiters	= 50;
	eps		= 1.0e-6;

/* Create the MIL tensor */

	*(eigvals+0) = ae;
	*(eigvals+1) = de/2.0;
	*(eigvals+2) = fe/2.0;
	*(eigvals+3) = de/2.0;
	*(eigvals+4) = be;
	*(eigvals+5) = ee/2.0;
	*(eigvals+6) = fe/2.0;
	*(eigvals+7) = ee/2.0;
	*(eigvals+8) = ce;

/* Calculate the eigenvalues and eigenvectors of the matrix/tensor */

	status = eigJacobi(dim,eigvals,eigvecs,order,eps,&numiters,maxiters);

	if (status == 1) {
		fprintf(stderr,"\nCould not converge in eigJacobi");
	}

/* Determine the angles of orientation for the structure */

	orient1a = atan2(*(eigvecs+6), sqrt(pow(*(eigvecs+0),2.0) + 
						pow(*(eigvecs+3),2.0)));
	orient1b = atan2(*(eigvecs+3),*(eigvecs+0));

	orient2a = atan2(*(eigvecs+7), sqrt(pow(*(eigvecs+1),2.0) +
						pow(*(eigvecs+4),2.0)));
	orient2b = atan2(*(eigvecs+4),*(eigvecs+1));

	orient3a = atan2(*(eigvecs+8), sqrt(pow(*(eigvecs+2),2.0) +
						pow(*(eigvecs+5),2.0)));
	orient3b = atan2(*(eigvecs+5),*(eigvecs+2));	

/* Convert the angles from radians -> degrees */

	orient1a *= (180.0/M_PI);
	if(orient1a < 0.0) orient1a += 180.0;
	orient1b *= (180.0/M_PI);
	if(orient1b < 0.0) orient1b += 180.0;
	orient2a *= (180.0/M_PI);
	if(orient2a < 0.0) orient2a += 180.0;
	orient2b *= (180.0/M_PI);
	if(orient2b < 0.0) orient2b += 180.0;
	orient3a *= (180.0/M_PI);
	if(orient3a < 0.0) orient3a += 180.0;
	orient3b *= (180.0/M_PI);
	if(orient3b < 0.0) orient3b += 180.0;

/* Determine the principal values */

	prin1 = 1.000 / sqrt(fabs(*(eigvals+0)));  
	prin2 = 1.000 / sqrt(fabs(*(eigvals+4)));  
	prin3 = 1.000 / sqrt(fabs(*(eigvals+8)));  


/* Determine the degrees of anisotropy */

	prinmn = (prin1 + prin2 + prin3)/3.0;

	deg1 = (prin1 - prinmn)/prinmn;
	deg2 = (prin2 - prinmn)/prinmn;
	deg3 = (prin3 - prinmn)/prinmn;

return (EXIT_SUCCESS);
}


/************************************************************************/
/* Function: eigstats(eigvals, eigvecs)					*/
/*									*/
/*	This function determines the statistical independence of the 	*/
/*	eigenvalues and the confidence intervals for the principal MILs.*/
/*									*/
/*	The variance in each of the eigenvalues is determined from the	*/
/*	the variance of the estimate of the fit ellipsoid (Se) weighted */
/*	by the inverse of the appropriate Gaussian multiplier term, Cii.*/
/*	The inverse Gaussian multipliers are coefficients from the 	*/
/*	inverse	of the sum-of-squares matrix.  The test statistic is 	*/
/*	defined as:							*/
/*									*/
/*		t = |evali - evalj| / [(Cii+Cjj-2*Cij)*Se]^0.5		*/
/*									*/
/*	The confidence interval for each eigenvalue is given by:	*/
/*									*/
/*		-/+ t(0.05) * (Cii*Se^2)^0.5				*/
/*									*/
/*	and the 95% confidence intervals are calculated for the 	*/
/*	principcal MILs.						*/
/*									*/
/*	The statistical independence of the eigenvalues is tested at a	*/
/*	0.05 confidence level (alpha), 	assuming 128 MIL vectors 	*/
/*	(rotations).  The critical t-value for these conditions is 	*/
/*	2.429.  Eigenvalues are statistically independent if the test 	*/
/*	statistic is greater than the critical t-value.  The test 	*/
/*	statistics are calculated for each of eigenvalue 1 versus 2, 	*/
/*	eigenvalue 2 versus 3, and eigenvalue 1 versus 3.  If all three	*/
/*	eigenvalues are different, the structure is orthotropic.  If	*/
/*	two eigenvalues are statistically equivalent, the structure is	*/
/*	transversely isotropic.  If all three are equivalent, the 	*/
/*	structure is isotropic.	Note that this classification will	*/
/*      depend on the significance level you use to determine 		*/
/*	statistical independence 					*/
/*									*/
/************************************************************************/

eigstats(eigvals,eigvecs)
double	*eigvals,*eigvecs;		/* eigenvalues and eigenvectors	*/

{
double	x[3];				/* unit direction vector	*/
double	xprin[3];			/* e-vector dot direct vector	*/

double	a[MAXPTS][3];			/* A matrix in princ coords	*/
double	at[3][MAXPTS];			/* A transpose (At)		*/
double	ata[3][3];			/* At * A (AtA)			*/
double	c[3][3];			/* inverse of AtA ((AtA)^-1)	*/
double	b[MAXPTS];			/* vector of 1/MIL^2		*/

int	i,j,k;				/* loop variables		*/
int	degree;				/* degree of matrices (3)	*/

	degree	= 3;

/* Initialize matrices */

	for (i=0;i<MAXPTS;i++) {
		for (j=0;j<degree;j++) {
			a[i][j]  = 0.0;
			at[j][i] = 0.0;
		}
	}

	for (i=0;i<degree;i++) {
		for (j=0;j<degree;j++) {
			ata[i][j] = 0.0;
			c[i][j]   = 0.0;
		}
	}


/* Determine A matrix in principal coordinates */

	for (i=0;i<MAXPTS;i++) {

		x[0] = cos(rot2[i])*sin(rot1[i]);
		x[1] = sin(rot2[i])*sin(rot1[i]);
		x[2] = cos(rot1[i]);

		for (j=0;j<degree;j++)
			xprin[j] = (*(eigvecs+0+j) * x[0]) + 
			(*(eigvecs+3+j) * x[1]) + (*(eigvecs+6+j) * x[2]);

		a[i][0] = xprin[0] * xprin[0];
		a[i][1] = xprin[1] * xprin[1];
		a[i][2] = xprin[2] * xprin[2];

		b[i] = 1.0/(*(mil+i)**(mil+i));
	}


/* Determine At of AtAx = Atb */

	for (i=0;i<MAXPTS;i++) {
		at[0][i] = a[i][0];
		at[1][i] = a[i][1];
		at[2][i] = a[i][2];
	}


/* Determine the sum-of-squares matrix*/

	for (i=0;i<degree;i++) {
		for (j=0;j<degree;j++) {
			for (k=0;k<MAXPTS;k++) {
				ata[i][j] += at[i][k]*a[k][j];
			}
		}
	}


/* Determine the inverse of the sum-of-squares matrix, c */

	linInverseMatrix (degree,degree,ata,c);


/* t value calculation  t[1 - alpha/2m; n - k - 1]			*/
/*	alpha	= confidence level	= 0.05				*/
/*	m	= # of tests		= 3				*/
/*	n	= # of observations	= 128				*/
/*	k	= # of variables	= 6  				*/
/*									*/
/*	t[0.99167; 121] = 2.429						*/

	tval = 2.429;


/* Determine test statistics for eigenvalue independence */

	*(test+0) = (fabs(*(eigvals+0) - *(eigvals+4))) / 
			(sqrt(variance*(c[0][0] + c[1][1] - 2*c[0][1])));

	*(test+1) = (fabs(*(eigvals+4) - *(eigvals+8))) / 
			(sqrt(variance*(c[1][1] + c[2][2] - 2*c[1][2])));

	*(test+2) = (fabs(*(eigvals+0) - *(eigvals+8))) / 
			(sqrt(variance*(c[0][0] + c[2][2] - 2*c[0][2])));


/* Determine confidence intervals for individual eigenvalues */

	for (i=0;i<degree;i++) {
		*(ci+i) = tval * sqrt(c[i][i] * variance);
	}

/* Determine the confidence intervals for the principal MILs */

	*(princi+0) = copysign((1.0 / sqrt(fabs(*(eigvals+0) +
 					*(ci+0)))),(*(eigvals+0) + *(ci+0)));
	*(princi+1) = copysign((1.0 / sqrt(fabs(*(eigvals+0) -
					*(ci+0)))),(*(eigvals+0) - *(ci+0)));
	*(princi+2) = copysign((1.0 / sqrt(fabs(*(eigvals+4) +
					*(ci+1)))),(*(eigvals+4) + *(ci+1)));
	*(princi+3) = copysign((1.0 / sqrt(fabs(*(eigvals+4) -
					*(ci+1)))),(*(eigvals+4) - *(ci+1)));
	*(princi+4) = copysign((1.0 / sqrt(fabs(*(eigvals+8) +
					*(ci+2)))),(*(eigvals+8) + *(ci+2)));
	*(princi+5) = copysign((1.0 / sqrt(fabs(*(eigvals+8) -
					*(ci+2)))),(*(eigvals+8) - *(ci+2)));

return (EXIT_SUCCESS);
}


/************************************************************************/
/* Function: sortdata(eigvals)						*/
/*									*/
/*	This function sorts the principal MIL data so that the results 	*/
/*	output to stdout and stderr are ordered from longest to		*/
/*	shortest principal MIL.						*/
/*									*/
/************************************************************************/

sortdata(eigvals)
double	*eigvals;
{

double	ptemp1,ptemp2,ptemp3;
double	otemp1a,otemp1b,otemp2a,otemp2b,otemp3a,otemp3b;
double	dtemp1,dtemp2,dtemp3;

double	*ttemp,*pctemp;

int	i;

	ttemp  = (double *)calloc(3,sizeof(double));
	pctemp = (double *)calloc(6,sizeof(double));

	ptemp1 = prin1;
	ptemp2 = prin2;
	ptemp3 = prin3;

	otemp1a = orient1a;
	otemp1b = orient1b;
	otemp2a = orient2a;
	otemp2b = orient2b;
	otemp3a = orient3a;
	otemp3b = orient3b;

	dtemp1 = deg1;
	dtemp2 = deg2;
	dtemp3 = deg3;

	for(i=0;i<3;i++) {
		*(ttemp+i) = *(test+i);
	}

	for(i=0;i<6;i++) {
		*(pctemp+i) = *(princi+i);
	}


/* Order the data from longest to shortest principal axis */

	if (ptemp1>ptemp2 && ptemp1>ptemp3) {
		if (ptemp2>ptemp3) {
			prin1 = copysign(ptemp1,*(eigvals+0));
			prin2 = copysign(ptemp2,*(eigvals+4));
			prin3 = copysign(ptemp3,*(eigvals+8));

			orient1a = otemp1a;
			orient1b = otemp1b;
			orient2a = otemp2a;
			orient2b = otemp2b;
			orient3a = otemp3a;
			orient3b = otemp3b;

			deg1 = dtemp1;
			deg2 = dtemp2;
			deg3 = dtemp3;

			*(princi+0) = *(pctemp+0);
			*(princi+1) = *(pctemp+1);
			*(princi+2) = *(pctemp+2);
			*(princi+3) = *(pctemp+3);
			*(princi+4) = *(pctemp+4);	
			*(princi+5) = *(pctemp+5);

			*(test+0) = *(ttemp+0);
			*(test+1) = *(ttemp+1);
			*(test+2) = *(ttemp+2);

		}
		else {
			prin1 = copysign(ptemp1,*(eigvals+0));
			prin2 = copysign(ptemp3,*(eigvals+8));
			prin3 = copysign(ptemp2,*(eigvals+4));

			orient1a = otemp1a;
			orient1b = otemp1b;
			orient2a = otemp3a;
			orient2b = otemp3b;
			orient3a = otemp2a;
			orient3b = otemp2b;

			deg1 = dtemp1;
			deg2 = dtemp3;
			deg3 = dtemp2;

			*(princi+0) = *(pctemp+0);
			*(princi+1) = *(pctemp+1);
			*(princi+2) = *(pctemp+4);
			*(princi+3) = *(pctemp+5);
			*(princi+4) = *(pctemp+2);
			*(princi+5) = *(pctemp+3);

			*(test+0) = *(ttemp+2);
			*(test+1) = *(ttemp+1);
			*(test+2) = *(ttemp+0);
		}
	}

	if (ptemp2>ptemp1 && ptemp2>ptemp3) {
		if (ptemp1>ptemp3) {
			prin1 = copysign(ptemp2,*(eigvals+4));
			prin2 = copysign(ptemp1,*(eigvals+0));
			prin3 = copysign(ptemp3,*(eigvals+8));

			orient1a = otemp2a;
			orient1b = otemp2b;
			orient2a = otemp1a;
			orient2b = otemp1b;
			orient3a = otemp3a;
			orient3b = otemp3b;

			deg1 = dtemp2;
			deg2 = dtemp1;
			deg3 = dtemp3;

			*(princi+0) = *(pctemp+2);
			*(princi+1) = *(pctemp+3);
			*(princi+2) = *(pctemp+0);
			*(princi+3) = *(pctemp+1);
			*(princi+4) = *(pctemp+4);
			*(princi+5) = *(pctemp+5);

			*(test+0) = *(ttemp+0);
			*(test+1) = *(ttemp+2);
			*(test+2) = *(ttemp+1);
		}
		else {
			prin1 = copysign(ptemp2,*(eigvals+4));
			prin2 = copysign(ptemp3,*(eigvals+8));
			prin3 = copysign(ptemp1,*(eigvals+0));

			orient1a = otemp2a;
			orient1b = otemp2b;
			orient2a = otemp3a;
			orient2b = otemp3b;
			orient3a = otemp1a;
			orient3b = otemp1b;

			deg1 = dtemp2;
			deg2 = dtemp3;
			deg3 = dtemp1;

			*(princi+0) = *(pctemp+2);
			*(princi+1) = *(pctemp+3);
			*(princi+2) = *(pctemp+4);
			*(princi+3) = *(pctemp+5);
			*(princi+4) = *(pctemp+0);
			*(princi+5) = *(pctemp+1);

			*(test+0) = *(ttemp+1);
			*(test+1) = *(ttemp+2);
			*(test+2) = *(ttemp+0);
		}
	}

	if (ptemp3>ptemp1 && ptemp3>ptemp2) {
		if (ptemp1>ptemp2) {
			prin1 = copysign(ptemp3,*(eigvals+8));
			prin2 = copysign(ptemp1,*(eigvals+0));
			prin3 = copysign(ptemp2,*(eigvals+4));

			orient1a = otemp3a;
			orient1b = otemp3b;
			orient2a = otemp1a;
			orient2b = otemp1b;
			orient3a = otemp2a;
			orient3b = otemp2b;

			deg1 = dtemp3;
			deg2 = dtemp1;
			deg3 = dtemp2;

			*(princi+0) = *(pctemp+4);
			*(princi+1) = *(pctemp+5);
			*(princi+2) = *(pctemp+0);
			*(princi+3) = *(pctemp+1);
			*(princi+4) = *(pctemp+2);
			*(princi+5) = *(pctemp+3);

			*(test+0) = *(ttemp+2);
			*(test+1) = *(ttemp+0);
			*(test+2) = *(ttemp+1);
		}
		else {
			prin1 = copysign(ptemp3,*(eigvals+8));
			prin2 = copysign(ptemp2,*(eigvals+4));
			prin3 = copysign(ptemp1,*(eigvals+0));

	
			orient1a = otemp3a;
			orient1b = otemp3b;
			orient2a = otemp2a;
			orient2b = otemp2b;
			orient3a = otemp1a;
			orient3b = otemp1b;

			deg1 = dtemp3;
			deg2 = dtemp2;
			deg3 = dtemp1;

			*(princi+0) = *(pctemp+4);
			*(princi+1) = *(pctemp+5);
			*(princi+2) = *(pctemp+2);
			*(princi+3) = *(pctemp+3);
			*(princi+4) = *(pctemp+0);
			*(princi+5) = *(pctemp+1);

			*(test+0) = *(ttemp+1);
			*(test+1) = *(ttemp+0);
			*(test+2) = *(ttemp+2);
		}
	}

return (EXIT_SUCCESS);
}


/*  The next four subroutines are used with permission of Triaxis, Inc,
Los Alamos, NM 87544.  They are from  Version 2.0 of their Math++ library */
/**********************************************************************
 *		Copyright (c) Triakis, Inc. 1989		      *
 *								      *
 *  Function: linSystemSolve                                          *
 *                                                                    *
 *  Version:  2.0                                                     *
 *                                                                    *
 *  Purpose:  Solution of linear system: Ax = b.                      *
 *            Used in conjunction with function linLuDecomp.          *
 *            Not used if linLuDecomp has detected a singularity.     *
 *                                                                    *
 *            linSystemSolve does the forward elimination             *
 *            and back substitution.  Splitting the process into      *
 *            two parts allows solutions for multiple b vectors       *
 *            with a minimum of recomputation.                        *
 *                                                                    *
 *  Returns:  Status ( 0 = success and 1 = failure ).                 *
 *                                                                    *
 *  Call:     linSystemSolve( dim, n, a, b, pivot );                  *
 *                                                                    *
 *      input:                                                        *
 *                                                                    *
 *          dim = Maximum dimension of A in main.                     *
 *            n = order of matrix.                                    *
 *            a = Triangularized matrix obtained from linLuDecomp.    *
 *            b = Right hand side vector.                             *
 *        pivot = Pivot vector obtained from linLuDecomp.             *
 *                                                                    *
 *      output:                                                       *
 *                                                                    *
 *            b = solution vector, x.                                 *
 *                                                                    *
 *  External functions:                                               *
 *            none.                                                   *
 *                                                                    *
 * Algorithm: Forsythe, Malcolm, and Moler from the text              *
 *            "Computer Methods for Mathematical Computations",       *
 *            1977, Prentice Hall.                                    *
 *                                                                    *
 *            Coded in structured 'C' by Tim Julian and Eric Blommer. *
 *                                                                    *
 **********************************************************************/

						/* Make this function capable of any size of matrix 'A'. */
#undef  A
#define A(X,Y)  *(a + (X) * dim + (Y))
						/* Define program constants. */
#undef  GOOD
#define GOOD  0
#undef  BAD
#define BAD   1

int 
linSystemSolve (dim, order, a, b, pivot)

    int		dim,
		order;
    double	*a,
		b[];
    int		pivot[];

{

    int             i,
                    k,
                    m;
    double          temp;

						/* Forward elimination, trivial case. */
    if (order == 1) {
	if (A (0, 0) == 0.0) {
	    return BAD;
	}
	else {
	    b[0] /= A (0, 0);
	    return GOOD;
	}
    }
						/* Forward elimination, non-trivial case. */
    for (k = 0; k < order - 1; ++k) {
	m = pivot[k];
	temp = b[m];
	b[m] = b[k];
	b[k] = temp;
	for (i = k + 1; i < order; ++i) {
	    b[i] += A (i, k) * temp;
	}
    }
						/* Back substitution */
    for (k = order - 1; k >= 0; --k) {
	if (A (k, k) == 0) {
	    return BAD;
	}
	b[k] /= A (k, k);
	temp = -b[k];
	for (i = 0; i < k; ++i) {
	    b[i] += A (i, k) * temp;
	}
    }

    return GOOD;

} /* End linSystemSolve */

/**********************************************************************
 *		Copyright (c) Triakis, Inc. 1989		      *
 *								      *
 *  Function: linLuDecomp                                             *
 *                                                                    *
 *  Version:  2.0                                                     *
 *                                                                    *
 *  Purpose:  Decomposes a real general matrix by LU decomposition    *
 *            and estimates the condition of the matrix.              *
 *                                                                    *
 *            Function linSystemSolve computes solution to linear     *
 *	      systems after decomposition by this function.           *
 *                                                                    *
 *  Returns:  Condnum = An estimate of the condition number of "A".   *
 *            A good rule of thumb in using condnum is that by taking *
 *            the log10 of the condnum and subtracting that value     *
 *            from the number of double precision significant digits  *
 *            availible, the remainder is the approximate number of   *
 *            significant digits in your answer.                      *
 *            Condnum is set to 3.33e+33 if exact singularity         *
 *            is detected.                                            *
 *                                                                    *
 *  Call:     status = linLuDecomp( dim, order, a, pivot, &condnum ); *
 *                                                                    *
 *       input:                                                       *
 *          dim  = Maximum dimensioned columns of "A".                *
 *         order = Order of the matrix  (actual size of "A").         *
 *            a  = Square matrix to be decomposed.                    *
 *                                                                    *
 *      output:                                                       *
 *            a = Contains an upper triangular matrix u and a         *
 *                lower triangular matrix of multipliers.             *
 *      condnum = An estimate of the condition number of A.           *
 *                For the linear system A*x = b, changes in A and b   *
 *                may cause changes condnum times as large in x.      *
 *                Condnum is set to 3.33e+33 if exact singularity     *
 *                is detected.                                        *
 *        pivot = The pivot vector.                                   *
 *                pivot[i] is the index of the ith pivot row.         *
 *                pivot[order-1] is (-1) ** (number of interchanges)  *
 *                                                                    *
 *  External functions:                                               *
 *            linSystemSolve(),                                       *
 *            fabs(),                                                 *
 *            calloc(),                                               *
 *            free().                                                 *
 *                                                                    *
 * Algorithm: Forsythe, Malcolm, and Moler from the text              *
 *            "Computer Methods for Mathematical Computations",       *
 *            1977, Prentice Hall( pages 30 - 62 ).                   *
 *            Also see the text "Linpack Users' Guide" by Dongarra,   *
 *            Moler, Bunch, and Stewart, 1979 by Siam.                *
 *            ( pages 1.1 - 1.34 and source listings in back )        *
 *            Note the routines dgeco, dgefa, dgedi, dgesl,           *
 *            ddot, idamax, dscal, daxpy, and dasum.                  *
 *            Corresponding single precision routines are sgeco,      *
 *            sgefa, sgedi, sgesl, sdot, isamax, sscal, saxpy,        *
 *            and sasum.                                              *
 *                                                                    *
 *  Misc:    The determinant of A can be obtained on output by        *
 *           det_a = pivot[order-1] * A(0,0) * A(1,1) *               *
 *                                ... * A(order-1,order-1).           *
 *           See the book Computer Methods for Mathematical           *
 *           Computations by Forsythe, Malcolm and Moler              *
 *           pages 30 - 62 for further explanation.                   *
 *                                                                    *
 *           Coded in structured 'C' by Tim Julian and Eric Blommer.  *
 *                                                                    *
 **********************************************************************/

						/* The following lines allow the function to be separately compiled */
#ifndef NULL
#include <stdio.h>
#include <math.h>
#endif
						/* Allow the function to accept any order for matrix 'A' */
#undef  A
#define A(X,Y)  *(a + (X)*dim + (Y))
						/* Define program constants. */
#undef  GOOD
#define GOOD  0
#undef  BAD
#define BAD   1
#undef  ABEND
#define ABEND 1

int 
linLuDecomp (dim, order, a, pivot, condnum)

    int             dim,
                    order,
                    pivot[];
    double         *a,
                   *condnum;

{
 						/* Local variable declarations. */
    double         *work,
                    ek,
                    temp,
                    anorm,
                    ynorm,
                    znorm;
    int             i,
                    j,
                    k,
                    kplus1,
                    m,
                    flag;
						/* Function declarations. */
    char           *calloc ();
    int		   linSystemSolve ();

						/* Check that "a" matrix order not larger than maximum size */
    if (order > dim) {
	(void) fprintf (stderr, " In linLuDecomp: order > Max. dim.\n");
	(void) exit (ABEND);
    }
						/* Check that "a" matrix order not smaller than one. */
    if (order < 1) {
	(void) fprintf (stderr, " In linLuDecomp: order < 1\n");
	(void) exit (ABEND);
    }

						/* Dynamically allocate work vector appropriate for problem. */
    if ((work = (double *) calloc (order, sizeof (double))) == NULL) {
	(void) fprintf (stderr, " In linLuDecomp: can't allocate work vector.\n");
	(void) exit (ABEND);
    }

    pivot[order - 1] = 1;
						/* Treat the trivial case. */
    if (order == 1) {
	if (A (0, 0) != 0.0) {
	    *condnum = 1.0;
	    free (work);
	    return GOOD;
	}
	else {
						/* Exact singularity detected for trivial case */
	    *condnum = 3.33e33;
	    free (work);
	    return BAD;
	}
    }
						/* Compute the 1-norm of matrix 'a'.                    */
						/* The 1-norm of matrix 'a' is defined as the largest   */
						/* sum of the absolute values of each of the elements   */
						/* of a single column vector across all column vectors. */
    anorm = 0.0;
    for (j = 0; j < order; ++j) {
	temp = 0.0;
	for (i = 0; i < order; ++i) {
	    temp += fabs (A (i, j));
	}
	if (temp > anorm) {
	    anorm = temp;
	}
    }
						/* Gaussian elimination with partial pivoting. */
						/* See Linpack function dgefa or sgefa         */
    for (k = 0; k < order - 1; ++k) {
	kplus1 = k + 1;
						/* Find pivot. */
	m = k;
	for (i = kplus1; i < order; ++i) {
	    if (fabs (A (i, k)) > fabs (A (m, k))) {
		m = i;
	    }
	}
	pivot[k] = m;
						/* See Linpack function dgedi or sgedi */
	if (m != k) {
	    pivot[order - 1] = -pivot[order - 1];
	}
						/* Skip step if pivot is zero. */
	if (A (m, k) != 0.0) {

	    if (m != k) {
		temp = A (m, k);
		A (m, k) = A (k, k);
		A (k, k) = temp;
	    }
						/* Compute multipliers. */
	    temp = -1.0 / A (k, k);

	    for (i = kplus1; i < order; ++i) {
		A (i, k) *= temp;
	    }
						/* Row elimination with column indexing */
	    for (j = kplus1; j < order; ++j) {
		temp = A (m, j);
		A (m, j) = A (k, j);
		A (k, j) = temp;
		if (temp != 0.0) {
		    for (i = kplus1; i < order; ++i) {
			A (i, j) += A (i, k) * temp;
		    }
		}
	    }
	}
    } /* End k loop. */
						/* Calculate condition number of matrix.                         */
						/* condnum = (1-norm of a)*(an estimate of 1-norm of a-inverse)  */
						/* Estimate obtained by one step of inverse iteration for the    */
						/* small singular vector.  This involves solving two systems     */
						/* of equations, (a-transpose)*y = e and a*z = y where e         */
						/* is a vector of +1 or -1 chosen to cause growth in y.          */
						/* Estimate = (1-norm of z)/(1-norm of y)                        */
						/* First, solve (a-transpose)*y = e                              */
    for (k = 0; k < order; ++k) {
	temp = 0.0;
	if (k != 0) {
	    for (i = 0; i < k - 1; ++i) {
		temp += A (i, k) * work[i];
	    }
	}
	if (temp < 0.0) {
	    ek = -1.0;
	}
	else {
	    ek = 1.0;
	}
	if (A (k, k) == 0.0) {
						/* Exact singularity detected for non-trivial case */
	    *condnum = 3.33e33;
	    free (work);
	    return BAD;
	}
	work[k] = -(ek + temp) / A (k, k);
    }
    for (k = order - 2; k >= 0; --k) {
	temp = 0.0;
	for (i = k + 1; i < order; ++i) {
	    temp += A (i, k) * work[k];
	}
	work[k] = temp;
	m = pivot[k];
	if (m != k) {
	    temp = work[m];
	    work[m] = work[k];
	    work[k] = temp;
	}
    }

    ynorm = 0.0;
    for (i = 0; i < order; ++i) {
	ynorm += fabs (work[i]);
    }
						/* Solve a*z = y. */
    flag = linSystemSolve (dim, order, a, work, pivot);
    if (flag == 1) {
	(void) fprintf (stderr, "In linLuDecomp: error detected in function linSystemSolve.\n");
	(void) exit (ABEND);
    }

    znorm = 0.0;
    for (i = 0; i < order; ++i) {
	znorm += fabs (work[i]);
    }
						/* Estimate condition number. */
    *condnum = anorm * znorm / ynorm;
    if (*condnum < 1.0) {
	*condnum = 1.0;
    }

    free (work);

    return GOOD;

} /* End linLuDecomp */

/**********************************************************************
 *		Copyright (c) Triakis, Inc. 1989		      *
 *								      *
 *  Function: linInverseMatrix                                        *
 *                                                                    *
 *  Version:  2.0                                                     *
 *                                                                    *
 *  Purpose:  Determines the inverse of a matrix via LU decomposition.*
 *                                                                    *
 *  Returns:  Nothing of value.                                       *
 *                                                                    *
 *  Call:     linInverseMatrix( dim, order, a, in );                  *
 *                                                                    *
 *       input:                                                       *
 *          dim = maximum dimensioned columns of "A".                 *
 *         order = order of the matrix  (actual size of "A").         *
 *            a = square matrix to be decomposed.                     *
 *                                                                    *
 *         output:                                                    *
 *           in = Inverse of matrix "a".                              *
 *                                                                    *
 *  External functions:                                               *
 *            linSystemSolve(),					      *
 *            linLuDecomp(),                                          *
 *            calloc(),                                               *
 *            free().                                                 *
 *                                                                    *
 * Algorithm: Tim Julian and Eric Blommer.                            *
 *                                                                    *
 **********************************************************************/

						/* The following lines allow the function to be separately compiled. */
#ifndef NULL
#include <stdio.h>
#include <math.h>
#endif
						/* Cause the function to accept any size for arrays "A" and "IN" */
#undef  A
#define A(X,Y)  *(a + (X)*dim + (Y))
#undef  IN
#define IN(X,Y)  *(in + (X)*dim + (Y))
#define BIGNUM  1.0E15

#undef  ABEND
#define ABEND 1

void 
linInverseMatrix (dim, order, a, in)

    int             dim,
                    order;
    double         *a,
                   *in;
{
    int             i,
                    j;				/* work variables */
    int             status;			/* return flag to check functions */
    int            *pivot;			/* pivoting from function linLuDecomp */
    int             linLuDecomp ();		/* decompose matrix a into LU matrix */
    int             linSystemSolve ();		/* solve system of simultaneous eqns. */
    double         *b;				/* B vector */
    double          condnum;			/* matrix condition number */
    char           *calloc ();

    if (order > dim) {
	(void) fprintf (stderr, "Inverse: order is greater than dim!\n");
	(void) exit (ABEND);
    }

						/* allocate work vectors */

    if ((b = (double *) calloc (order, sizeof (double))) == NULL) {
	(void) fprintf (stderr, " In linInverseMatrix, can't allocate b vector.\n");
	(void) exit (ABEND);
    }

    if ((pivot = (int *) calloc (order, sizeof (int))) == NULL) {
	(void) fprintf (stderr, " In linInverseMatrix, can't allocate pivot vector.\n");
	(void) exit (ABEND);
    }

						/* Initialize data values */
    condnum = 0.0;

						/* decompose matrix */
    status = linLuDecomp (dim, order, a, pivot, &condnum);

    if (status == 1) {
	(void) fprintf (stderr, "In Inverse, error detected in function linLuDecomp.\n");
	(void) exit (ABEND);
    }

    if (condnum > BIGNUM) {
	(void) fprintf (stderr, "Inverse: Matrix singular to double precision.\n");
	(void) fprintf (stderr, "Condition number : %e\n", condnum);
	(void) exit (ABEND);
    }
						/* backsolve matrix */
    for (i = 0; i < order; ++i) {
	for (j = 0; j < order; ++j) {
	    if (j == i) {
		b[j] = 1.0;
	    }
	    else {
		b[j] = 0.0;
	    }
	}
	status = linSystemSolve (dim, order, a, b, pivot);

	if (status == 1) {
	    (void) fprintf (stderr, "In Inverse, error detected in function linSystemSolve.\n");
	    (void) exit (ABEND);
	}
	for (j = 0; j < order; ++j) {
	    IN (j, i) = b[j];
	}
    }

    free (b);
    free (pivot);

} /* end linInverseMatrix */

/**********************************************************************
 *		Copyright (c) Triakis, Inc. 1989		      *
 *                                                                    *
 *  Function: eigJacobi                                               *
 *                                                                    *
 *  Version:  2.0                                                     *
 *                                                                    *
 *  Purpose:  Computes the eigenvalues and eigenvectors of a real     *
 *            symmetric matrix via the Jacobi method.                 *
 *                                                                    *
 *  Returns:  status = (0 = success, 1 = failure)                     *
 *                                                                    *
 *  Call:     status =  eigJacobi(  d, e1, e2, o, e, it, mit ), where:*
 *                                                                    *
 *       input:                                                       *
 *              d = maximum dimension of matrix eigenvalue matrix     *
 *             e1 = eigenvalue matrix.                                *
 *             e2 = eigenvector matrix.                               *
 *              o = order of the matrices                             *
 *              e = convergence criterion                             *
 *            mit = maximum number of iterations allowed              *
 *                                                                    *
 *      output:                                                       *
 *             e1 = eigenvalue matrix.                                *
 *             e2 = eigenvector matrix.                               *
 *             it = actual number of iterations                       *
 *                                                                    *
 *  External functions:                                               *
 *            fabs(),                                                 *
 *            calloc(),                                               *
 *            free(),                                                 *
 *            sqrt().                                                 *
 *                                                                    *
 *  Author(s): Tim Julian and Eric Blommer                            *
 *                                                                    *
 *  Misc:    None.                                                    *
 *                                                                    *
 *           Coded in structured 'C' by Tim Julian and Eric Blommer.  *
 *                                                                    *
 **********************************************************************/

						/* The following four lines allow this function */
 						/* to be compiled separately                    */
#ifndef NULL
#include <stdio.h>
#include <math.h>
#endif

#undef  BAD
#define BAD   1
#undef  GOOD
#define GOOD  0
#undef  ABEND
#define ABEND 1
						/* the following statements allow any size of eigenvalues or vectors */
						/* to be passed to this function                                     */
#define EIGVALS(i,j) *(eigvals+(i)*dim+(j))
#define EIGVECS(i,j) *(eigvecs+(i)*dim+(j))

int 
eigJacobi (dim, eigvals, eigvecs, order, eps, numiters, maxiters)

    int             dim;
    double         *eigvals,
                   *eigvecs;
    double          eps;
    int            *numiters,
                    maxiters,
                    order;
{

    double          vo,
                    f,
                    u,
                    uf,
                    alpha,
                    beta,
                    c,
                    s;
    double         *temp1,
                   *temp2,
                   *temp3;
    int             i,
                    j,
                    l,
                    m,
                    p,
                    q;
    char           *calloc ();

						/* dynamically allocate temporary vectors */
    if ((temp1 = (double *) calloc (dim, sizeof (double))) == NULL) {
	(void) fprintf (stderr, "Can't allocate vector temp1\n");
	(void) exit (ABEND);
    }

    if ((temp2 = (double *) calloc (dim, sizeof (double))) == NULL) {
	(void) fprintf (stderr, "Can't allocate vector temp2\n");
	(void) exit (ABEND);
    }

    if ((temp3 = (double *) calloc (dim, sizeof (double))) == NULL) {
	(void) fprintf (stderr, "Can't allocate vector temp3\n");
	(void) exit (ABEND);
    }
						/* make identity matrix */
    for (i = 0; i < order; ++i) {
	for (j = 0; j <= i; ++j) {
	    if (i == j) {
		EIGVECS (i, j) = 1.0;
	    }
	    else {
		EIGVECS (i, j) = EIGVECS (j, i) = 0.0;
	    }
	}
    }

    vo = 0.0;
    for (i = 0; i < order; ++i) {
	for (j = 0; j < order; ++j) {
	    if (i != j) {
		vo += EIGVALS (i, j) * EIGVALS (i, j);
	    }
	}
    }

    u = sqrt (vo) / (double) order;
    uf = eps * u;

    for (uf = eps * u; uf < u; u /= order) {
	for (l = 0; l < order - 1; ++l) {
	    for (m = l + 1; m < order; ++m) {
		if (fabs (EIGVALS (l, m)) >= u) {
		    p = l;
		    q = m;
		    for (i = 0; i < order; ++i) {
			temp1[i] = EIGVALS (p, i);
			temp2[i] = EIGVALS (i, p);
			temp3[i] = EIGVECS (i, p);
		    }
		    f = EIGVALS (p, p);
		    alpha = 0.5 * (EIGVALS (p, p) - EIGVALS (q, q));
		    beta = sqrt (EIGVALS (p, q) * EIGVALS (p, q) + alpha * alpha);
		    c = sqrt (0.5 + fabs (alpha) / (2.0 * beta));
		    s = alpha * (-EIGVALS (p, q)) / (2.0 * beta * fabs (alpha) * c);
		    for (j = 0; j < order; ++j) {
			if (j != p) {
			    if (j != q) {
				EIGVALS (p, j) = c * EIGVALS (p, j) - s * EIGVALS (q, j);
				EIGVALS (q, j) = s * temp1[j] + c * EIGVALS (q, j);
				EIGVALS (j, p) = c * EIGVALS (j, p) - s * EIGVALS (j, q);
				EIGVALS (j, q) = s * temp2[j] + c * EIGVALS (j, q);
			    }
			}
			EIGVECS (j, p) = c * EIGVECS (j, p) - s * EIGVECS (j, q);
			EIGVECS (j, q) = s * temp3[j] + c * EIGVECS (j, q);
		    }

		    EIGVALS (p, p) = c * c * EIGVALS (p, p) + s * s *
			    EIGVALS (q, q) - 2.0 * c * s * EIGVALS (p, q);
		    EIGVALS (q, q) = s * s * f + c * c * EIGVALS (q, q) +
			    2.0 * c * s * EIGVALS (p, q);
		    EIGVALS (p, q) = 0.0;
		    EIGVALS (q, p) = 0.0;

		    if (++*numiters > maxiters) {
			(void) fprintf (stderr,
					"number of iterations exceeded %d", maxiters);
			return BAD;
		    }
		}
	    }
	}
    }

    free (temp1);
    free (temp2);
    free (temp3);

    return GOOD;

} /* end eigJacobi */