/*	
	Mesh Generator for Neurons v2.9.2
   	Rogene M. Eichler West

	Copyright @ 1995-7 Rogene M. Eichler West
	All Rights Reserved

   	
   	Questions can be directed to the author at

   	     < rogene@med.umn.edu >
	     < rogene@bbf.uia.ac.be >
	     < rogene@bbb.caltech.edu >

	Documentation in accompanying README file

***** to do *****
	input formats
		eutectic v 4 and 5
		neurolucida
	output formats
		check neuron form
		seperate into smaller routines
	check banding
	soma slicing
*/

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

/* These defines are specific for NTS data files */
#define STARTLINE 7
#define MAXNO 100

typedef struct ptinf{
		int ptnum;
		char pttype[3];
		double pos[3];
		double rad;
		}PTINF;
typedef struct slicepts{
		int compno;
		double slicerad;
		double slicelength;
		double slicepos[3];
		int parent;
		int lchild;
		int rchild;
		int sib;
		}SLICEPTS;
struct ptinf *ptdat;
struct slicepts *slice;
double Rm, Ri, RadCrit,LCrit; 
double midpt[3], rangesoma[3];
double NEtmp[4];
int eutpoints;
int somapts = 0;
int priden = 0;
int endno = 0;
int branchno = 0;
int *somachild;
int track[3];
int parent = 1;
int sibling = 0;
int *hines;
int outtype[6];
int GENopt;
int CSA = 1;
int lastPT;
int firstPT;
int netmp;

void getdata(void);
void scaler(void);
void errorcheck(void);
void checker(void);
void reorder(void);
void dupcheck(void);
void balance(void);
void slicer(void);
void somaslicer(void);
double distance(double,double,double,double,double,double);
void writepoint(double[3],double,double);
void conserveSA(void);
void radiusslicer(void);
void lengthslicer(void);
int adjerror(int);
void topology(void);
void hinesnum(void);
void eichwestband(void);
void output(void);
void custom_out(double *, double *);
void genesis_out(double *, double *);

main()
{
	getdata();
	scaler();
	errorcheck(); 
	slicer();
	topology();
	hinesnum();
	eichwestband();
	output(); 
	printf("\nSuccessful Mesh Generation!\n");
}

/* ****************************************************************** */
void getdata(void)
{
int i;
char tmpchar[100], tmpchr[100];
char filename[40];
FILE *eutectdata;

/* 	Welcome, read data, and query for user specs.	*/	

	printf("Mesh Generator for Neurons v2.9.3\nby Rogene M. Eichler West\n\n\n");
	printf("Please enter the file name of the Eutectic Data\n");
	
	scanf("%s",filename);
	eutectdata = fopen(filename,"r");	
	if (eutectdata == NULL)
		{
		printf("File could not be found\n");
		exit(1);
		}		
		
	fscanf(eutectdata,"%s%s%d", tmpchar,tmpchar, &eutpoints);
	printf("\nYour input data has %d datapoints according to\n",eutpoints);
	printf("the header.\n");
	
/* Additional User defined criteria */
		printf("Enter Rm (Ohm*cm2)\t");
		scanf("%lf",&Rm);
		printf("Enter Ri (Ohm*cm)\t");
		scanf("%lf",&Ri);
		printf("Enter the maximum intracompartment radius range (um)\t");
		scanf("%lf",&RadCrit);
		RadCrit /= 10000.0;
		printf("Enter the maximum electrotonic length per compartment\t");
		scanf("%lf",&LCrit);
	if (eutpoints > 0){		
	ptdat = (struct ptinf *) malloc(sizeof(struct ptinf) * eutpoints);
	if (ptdat == NULL)
		{
		puts("Memory allocation error 1");
		exit(1);
	}
	}
		
		for (i=0; i<STARTLINE;i++) {
		fgets(tmpchar,MAXNO,eutectdata);
		}
	
		for (i=0; i<eutpoints; i++) {
		
		fgets(tmpchar,MAXNO,eutectdata);
		
		sscanf(tmpchar,"%d %s %s %lf %lf %lf %lf %s %s",
		&ptdat[i].ptnum, ptdat[i].pttype,tmpchr,&ptdat[i].pos[0],
		&ptdat[i].pos[1],&ptdat[i].pos[2],&ptdat[i].rad,tmpchr,tmpchr);

		/* in case the datafile has been altered */
                ptdat[i].ptnum = i + 1;

		/* ignore outline points */
		if ((strncmp(ptdat[i].pttype,"OS",2) == 0) ||
			(strncmp(ptdat[i].pttype,"OCP",3) == 0) ||
			(strncmp(ptdat[i].pttype,"OE",2) == 0) ||
			(strncmp(ptdat[i].pttype,"DS",2) == 0) ||
			(strncmp(ptdat[i].pttype,"DCP",3) == 0) ||
			(strncmp(ptdat[i].pttype,"DE",2) == 0)) {
			i--;
			eutpoints--;
			ptdat[i].ptnum = i - 1;
		}

		}

	fclose(eutectdata);

	printf("\nChoose from the following output formats:\n");
	printf("1\tCustom code\n");
	printf("2\tGenesis scripts\n");
	printf("3\tNeuron scripts\n");	
	scanf("%d", &outtype[0]);
		if (outtype[0] == 2) {
			printf("\n\tGenesis:\n");
			printf("\t 1\tAbsolute Coordinates\n");
			printf("\t 2\tRelative Coordinates\n");
			scanf("%d", &GENopt);
		}
	printf("\nChoose from the following options:\n");
	printf("Yes = 1, No = 0\n");
	printf("SGI objects?");
	scanf("%d",&outtype[1]);
	printf("Individual electrotonic lengths?");
	scanf("%d",&outtype[2]);
	printf("Electrotonic length from soma?");
	scanf("%d",&outtype[3]);
	printf("Distance from soma?");
	scanf("%d",&outtype[4]);
	printf("\n");
/* 
	Make optimization the default 
	printf("Optimization option:\n");
	printf("Implement Eichler West renumbering?");
	scanf("%d",&outtype[5]);
*/
	outtype[5] = 1;
 
}

/* ****************************************************************** */
void scaler(void)
{
int i;
int flg1 = 0;

/* convert from microns to centimeters
   divide diameter to obtain radius
*/
	for (i=0; i<eutpoints; i++) {
		ptdat[i].pos[0] /= 10000.0;
		ptdat[i].pos[1] /= 10000.0;
		ptdat[i].pos[2] /= 10000.0;
		ptdat[i].rad /= 20000.0;
		if (ptdat[i].rad == 0.0) {
			if (flg1 == 0) {
				printf("This file has radii of zero");
				printf("Changing to 0.5 micron.\n");	
				flg1 = 1;
			}
			printf("\t\t%d\n",ptdat[i].ptnum);
			ptdat[i].rad = .00005;
		}
	}		
}

/* ****************************************************************** */
void errorcheck(void)
{
	checker();
	reorder(); 
	dupcheck();
	balance();
}

/* ******************************************************************** */

void checker(void)
{
int i;
/* Convert nonrecognized point types.			
			Change ES to NE
			Change SB, FS to CP
			Convert all other points to CP
*/

	for (i=0; i<eutpoints; i++) {
		if ((strncmp(ptdat[i].pttype,"SOS",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"SCP",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"SOE",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"MTO",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"BTO",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"TTO",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"BAE",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"MAE",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"TAE",3) == 0) ||
	    	(strncmp(ptdat[i].pttype,"NE",2) == 0) ||
	    	(strncmp(ptdat[i].pttype,"BP",2) == 0) ||
	    	(strncmp(ptdat[i].pttype,"CP",2) == 0)) {
		    }				 
		else if (strncmp(ptdat[i].pttype,"ES",2) == 0) {
		    		strcpy(ptdat[i].pttype,"NE");
		    		}
		else if (strncmp(ptdat[i].pttype,"FS",2) == 0) {
		    		strcpy(ptdat[i].pttype,"CP");
		    		}
		else if (strncmp(ptdat[i].pttype,"SB",2) == 0) {
		    		strcpy(ptdat[i].pttype,"CP");
		    		}
		else {
		printf("\nWarning! The point type %s is not", ptdat[i].pttype);
		printf("recognized by this program.\n");
		printf("Please delete these point types from the input file");
		printf("and try again.\n");
		printf("Please contact the author regarding the failed\n");
		printf("point type so that a fix can be included in future\n");
		printf("versions of the	mesher.\n");
		exit(1);
		}
		}
}
/* ******************************************************************** */
void reorder(void)
{
int i,j,k,numbuff,lastne,nerequire,test;
double minx,maxx,miny,maxy,minz,maxz;
int mto = 0;
int bto = 0;
int tto = 0;
int mae = 0;
int bae = 0;
int tae = 0;
int nosoma = 0;
int *mtopt, *btopt, *ttopt, *maept, *baept, *taept;
int *list, *tmplist, *bufflist;
double y, prox,d;
typedef struct mtopairs{
		int to;
		int ae;
		double prox;
		}MTOPAIRS;
struct mtopairs *pmtop;
typedef struct btopairs{
		int to;
		int ae;
		double prox;
		}BTOPAIRS;
struct btopairs *pbtop;
typedef struct ttopairs{
		int to;
		int ae;
		double prox;
		}TTOPAIRS;
struct ttopairs *pttop;
typedef struct shortterm{
		int ptnum;
		char pttype[3];
		double pos[3];
		double rad;
		}SHORTTERM;
struct shortterm *starr;

/*  Examine the TTO/TAE, BTO/BAE, MTO/MAE points.
	Some are continuations in a new section of tissue.
	NTS has a function to realign the coordinates. However, if it
	was not used, then this reorder fix has no chance.
	
	ALL *TO points must be accounted for by either proximity
	to the soma or by proximity to a *AE point.
	
	No BTO or TTO on return - either MTO or CP
	No BAE, TAE, or MAE on return - all NE
	
	The assumption is that *AE out of the plane will
	always match *TO, with the same * (* = B,T, or M).
*/

/* 	Find the max radius and center of the soma.
	This will be used to determine which *TO points
	are primary dendrites off the soma. Save this info 
	for later, as it will be needed for the graphics file.
*/
		for (i=0; i<eutpoints; i++) {
		if ((strncmp(ptdat[i].pttype,"SOS",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"SCP",3) == 0) ||
		    (strncmp(ptdat[i].pttype,"SOE",3) == 0)) {
		    	somapts ++;
		    }
		if (strncmp(ptdat[i].pttype,"SOE",3) == 0) {
			nosoma++;
		}
		}
		if (nosoma > 1) {
			printf("The file contains multiple soma.\n");
			printf("\t\tSorry.\n");
                        exit(1);
                }
		if (somapts == 0) {
			printf("This file does not contain any soma points.\n");
			printf("\t\tSorry.\n");
			exit(1);
		}
		minx = ptdat[0].pos[0];
		maxx = ptdat[0].pos[0];
		miny = ptdat[0].pos[1];
		maxy = ptdat[0].pos[1];
		minz = ptdat[0].pos[2];
		maxz = ptdat[0].pos[2];
		for (i=1; i<somapts; i++) {
			if (minx < ptdat[i].pos[0]) {
				minx = ptdat[i].pos[0];
				}
			if (maxx > ptdat[i].pos[0]) {
				maxx = ptdat[i].pos[0];
				}
			if (miny < ptdat[i].pos[1]) {
				miny = ptdat[i].pos[1];
				}
			if (maxy > ptdat[i].pos[1]) {
				maxy = ptdat[i].pos[1];
				}
			if (minz < ptdat[i].pos[2]) {
				minz = ptdat[i].pos[2];
				}
			if (maxz > ptdat[i].pos[2]) {
				maxz = ptdat[i].pos[2];
				}
		}								
			midpt[0] = (minx + maxx)/2.0;
			midpt[1] = (miny + maxy)/2.0;
			midpt[2] = (minz + maxz)/2.0;
			rangesoma[0] = fabs(minx - maxx);
			rangesoma[1] = fabs(miny - maxy);
			rangesoma[2] = fabs(minz - maxz);
		
/* Count the *TO and *AE points */
		
		for (i=0; i<eutpoints; i++) {
		if (strncmp(ptdat[i].pttype,"MTO",3) == 0) {
		    		mto ++;
		    	}
		else if (strncmp(ptdat[i].pttype,"BTO",3) == 0) {
		    		bto ++;
		    	}
		else if (strncmp(ptdat[i].pttype,"TTO",3) == 0) {
					tto ++;
				}
		else if (strncmp(ptdat[i].pttype,"MAE",3) == 0) {
					mae ++;
				}
		else if (strncmp(ptdat[i].pttype,"BAE",3) == 0) {
					bae ++;
				}
		else if (strncmp(ptdat[i].pttype,"TAE",3) == 0) {
					tae ++;
				}
		}
/* Allocate memory based on the count */
if (mto > 0) {
	mtopt = (int *) malloc(sizeof(int) * mto);
	pmtop = (struct mtopairs *) malloc(sizeof(struct mtopairs) * mto);
	if ((mtopt == NULL) ||
		(pmtop == NULL)) {
		puts("Memory allocation error 2");
                exit(1);
                }
	}
if (bto > 0) {
        btopt = (int *) malloc(sizeof(int) * bto);
	pbtop = (struct btopairs *) malloc(sizeof(struct btopairs) * bto);	
        if ((btopt == NULL) ||
		(pbtop == NULL)) {
                puts("Memory allocation error 3");
                exit(1);
                }
        }
if (tto > 0) {
        ttopt = (int *) malloc(sizeof(int) * tto);
	pttop = (struct ttopairs *) malloc(sizeof(struct ttopairs) * tto);	
        if ((ttopt == NULL) ||
		(pttop == NULL)) {
                puts("Memory allocation error 4");
                exit(1);
                }
        }
if (mae > 0) {
        maept = (int *) malloc(sizeof(int) * mae);
        if (maept == NULL) {
                puts("Memory allocation error 5");
                exit(1);
                }
        }
if (bae > 0) {
        baept = (int *) malloc(sizeof(int) * bae);
        if (baept == NULL) {
                puts("Memory allocation error 6");
                exit(1);
                }
        }
if (tae > 0) {
        taept = (int *) malloc(sizeof(int) * tae);
        if (taept == NULL) {
                puts("Memory allocation error 7");
                exit(1);
                }
        }

/*	Place the address of each type in an array */

		mto = 0;
		bto = 0;
		tto = 0;
		mae = 0;
		bae = 0;
		tae = 0;		
		for (i=0; i<eutpoints; i++) {
		if (strncmp(ptdat[i].pttype,"MTO",3) == 0) {
			mtopt[mto] = ptdat[i].ptnum;
			mto ++;
		    }
		else if (strncmp(ptdat[i].pttype,"BTO",3) == 0) {
				btopt[bto] = (ptdat[i].ptnum);
		    	bto ++;
		    	}
		else if (strncmp(ptdat[i].pttype,"TTO",3) == 0) {
				ttopt[tto] = (ptdat[i].ptnum);
				tto ++;
				}
		else if (strncmp(ptdat[i].pttype,"MAE",3) == 0) {
				maept[mae] = (ptdat[i].ptnum);
				mae ++;
				}
		else if (strncmp(ptdat[i].pttype,"BAE",3) == 0) {
				baept[bae] = (ptdat[i].ptnum);
				bae ++;
				}
		else if (strncmp(ptdat[i].pttype,"TAE",3) == 0) {
				taept[tae] = (ptdat[i].ptnum);
				tae ++;
				}
		}

/*  It is assumed that any *TO point contiguous with the soma
	can be found within (midpt +/- rangesoma).
	If a *TO point is within the soma space, the point becomes
	a MTO and is taken off the list for the next section.
	Changing the value to -99 means that it is a primary dendrite.
*/
	d = .001;	
	for (i=0; i<mto; i++) {
		if ((ptdat[mtopt[i]].pos[0] < (midpt[0]+rangesoma[0]+d)) &&
		(ptdat[mtopt[i]-1].pos[0] > (midpt[0]-rangesoma[0]-d)) &&
		(ptdat[mtopt[i]-1].pos[1] < (midpt[1]+rangesoma[1]+d)) &&
		(ptdat[mtopt[i]-1].pos[1] > (midpt[1]-rangesoma[1]-d)) &&
		(ptdat[mtopt[i]-1].pos[2] < (midpt[2]+rangesoma[2]+d)) &&
		(ptdat[mtopt[i]-1].pos[2] > (midpt[2]-rangesoma[2]-d))) { 
			strcpy(ptdat[mtopt[i]-1].pttype,"MTO");
			mtopt[i] = -99;
		}
	}
	for (i=0; i<bto; i++) {
		if ((ptdat[btopt[i]].pos[0] < (midpt[0]+rangesoma[0]+d)) &&
		(ptdat[btopt[i]-1].pos[0] > (midpt[0]-rangesoma[0]-d)) &&
		(ptdat[btopt[i]-1].pos[1] < (midpt[1]+rangesoma[1]+d)) &&
		(ptdat[btopt[i]-1].pos[1] > (midpt[1]-rangesoma[1]-d)) &&
		(ptdat[btopt[i]-1].pos[2] < (midpt[2]+rangesoma[2]+d)) &&
		(ptdat[btopt[i]-1].pos[2] > (midpt[2]-rangesoma[2]-d))) {
			strcpy(ptdat[btopt[i]-1].pttype,"MTO");
			btopt[i] = -99;
			}
	}	
	for (i=0; i<tto; i++) {
		if ((ptdat[ttopt[i]].pos[0] < (midpt[0]+rangesoma[0]+d)) &&
		(ptdat[ttopt[i]-1].pos[0] > (midpt[0]-rangesoma[0]-d)) &&
		(ptdat[ttopt[i]-1].pos[1] < (midpt[1]+rangesoma[1]+d)) &&
		(ptdat[ttopt[i]-1].pos[1] > (midpt[1]-rangesoma[1]-d)) &&
		(ptdat[ttopt[i]-1].pos[2] < (midpt[2]+rangesoma[2]+d)) &&
		(ptdat[ttopt[i]-1].pos[2] > (midpt[2]-rangesoma[2]-d))) {
			strcpy(ptdat[ttopt[i]-1].pttype,"MTO");
			ttopt[i] = -99;
			}
	}	
	
/* initialize new arrays */
	for (i=0; i<mto; i++) {
		pmtop[i].to = -99;
		pmtop[i].ae = -99;
		pmtop[i].prox = -99;
	}
	for (i=0; i<bto; i++) {
		pbtop[i].to = -99;
		pbtop[i].ae = -99;
		pbtop[i].prox = -99;
	}
	for (i=0; i<tto; i++) {
		pttop[i].to = -99;
		pttop[i].ae = -99;
		pttop[i].prox = -99;
	}
	
/* 	For these remaining points, match up *TO with *AE	*/
		for (i=0; i<mto; i++) {
			prox = 1000.0;
			if (mtopt[i] != -99) {
				if (mae == 0) {
	printf("\nThis dataset contains a branch not originating\n");
	printf("at the soma, yet not continuing from a previous\n");
	printf("MAE ending. This code will not deal with such\n");
	printf("structures. See datapoint %d.\n",mtopt[i]);
	printf("\n\t\t\tSorry.");
	exit(1);
	}
		for (j=0; j<mae; j++) {
	y = distance(ptdat[mtopt[i]-1].pos[0],ptdat[maept[j]-1].pos[0],
		ptdat[mtopt[i]-1].pos[1], ptdat[maept[j]-1].pos[1],
		ptdat[mtopt[i]-1].pos[2], ptdat[maept[j]-1].pos[2]);
			if (y < prox) {
				prox = y;
				pmtop[i].to = (ptdat[mtopt[i]-1].ptnum);
				pmtop[i].ae = (ptdat[maept[j]-1].ptnum);
				pmtop[i].prox = prox;
				}
			}
			if (prox > .001) {
	printf("\nThe closest match for the branch originating at\n");
	printf("point %d is over 10 microns away. I doubt\n",mtopt[i]);
	printf("these points actually match and conclude there is\n");
	printf("no match for this structure. This code will not deal\n");
	printf("with such structures.\n");
	printf("Did you use the \"merge\" Eutectic function to join\n");
	printf("multiple tissue slide segments?\n");
	printf("\n\t\t\tSorry.\n");
	exit(1);
	}	
}		 
}
	for (i=0; i<mto; i++) {
		for (j=i+1; j<mto; j++) {
		if ((mtopt[i] != -99) || (mtopt[j] != -99)) {
				if (pmtop[i].ae == pmtop[j].ae) {
printf("\nTwo MTO points, %d and %d\n",mtopt[i],mtopt[j]);
printf("have chosen the same MAE point. I have not\n");
printf("has a chance to code a fix to this possibility.\n");
printf("Send mail to the author for a fix.\n");
printf("\n\t\t\tSorry.");
exit(1);
		}
	}
	}
}				
	for (i=0; i<bto; i++) {
		prox = 1000.0;
		if (btopt[i] != -99) {
			if (bae == 0) {
	printf("\nThis dataset contains a branch not originating\n");
	printf("at the soma, yet not continuing from a previous\n");
	printf("BAE ending. This code will not deal with such\n");
	printf("structures. See datapoint %d.\n", btopt[i]);
	printf("\n\t\t\tSorry.");
	exit(1);
		}
	for (j=0; j<bae; j++) {
	y = distance(ptdat[btopt[i]-1].pos[0],ptdat[baept[j]-1].pos[0],
		ptdat[btopt[i]-1].pos[1], ptdat[baept[j]-1].pos[1],
		ptdat[btopt[i]-1].pos[2], ptdat[baept[j]-1].pos[2]);
		if (y < prox) {
			prox = y;
			pbtop[i].to = (ptdat[btopt[i]-1].ptnum);
			pbtop[i].ae = (ptdat[baept[j]-1].ptnum);
			pbtop[i].prox = prox;
			}
			}
		if (prox > .001) {
	printf("\nThe closest match for the branch originating at\n");
	printf("point %d is over 10 microns away. I doubt\n",btopt[i]);
	printf("these points actually match and conclude there is\n");
	printf("no match for this structure. This code will not deal\n");
	printf("with such structures.\n");
        printf("Did you use the \"merge\" Eutectic function to join\n");
        printf("multiple tissue slide segments?\n");
        printf("\n\t\t\tSorry.\n");
	exit(1);
			}	
		}
	}	
	for (i=0; i<bto; i++) {
		for (j=i+1; j<bto; j++) {
			if ((btopt[i] != -99) || (btopt[j] != -99)) {
				if (pbtop[i].ae == pbtop[j].ae) {
	printf("\nTwo BTO points, %d and %d\n",btopt[i],btopt[j]);
	printf("have chosen the same BAE point. I have not\n");
	printf("has a chance to code a fix to this possibility.\n");
	printf("Send mail to the author for a fix.\n");
	printf("\n\t\t\tSorry.");
	exit(1);
			}
		}
	}
}			
	for (i=0; i<tto; i++) {
		prox = 1000.0;
		if (ttopt[i] != -99) {
			if (tae == 0) {
	printf("\nThis dataset contains a branch not originating\n");
	printf("at the soma, yet not continuing from a previous\n");
	printf("TAE ending. This code will not deal with such\n");
	printf("structures. See datapoint %d.\n", ttopt[i]);
	printf("\n\t\t\tSorry.");
	exit(1);
			}
	for (j=0; j<tae; j++) {
	y = distance(ptdat[ttopt[i]-1].pos[0],ptdat[taept[j]-1].pos[0],
		ptdat[ttopt[i]-1].pos[1], ptdat[taept[j]-1].pos[1],
		ptdat[ttopt[i]-1].pos[2], ptdat[taept[j]-1].pos[2]);
		if (y < prox) {
			prox = y;
			pttop[i].to = (ptdat[ttopt[i]-1].ptnum);
			pttop[i].ae = (ptdat[taept[j]-1].ptnum);
			pttop[i].prox = prox;
			}
		}
		if (prox > .001) {
	printf("\nThe closest match for the branch originating at\n");
	printf("point %d is over 10 microns away. I doubt\n",ttopt[i]);
	printf("these points actually match and conclude there is\n");
	printf("no match for this structure. This code will not deal\n");
	printf("with such structures.\n");
        printf("Did you use the \"merge\" Eutectic function to join\n");
        printf("multiple tissue slide segments?\n");
        printf("\n\t\t\tSorry.\n");
	exit(1);
		}
	}	
}			
	for (i=0; i<tto; i++) {
		for (j=i+1; j<tto; j++) {
			if ((ttopt[i] != -99) || (ttopt[j] != -99)) {
				if (pttop[i].ae == pttop[j].ae) {
	printf("\nTwo TTO points, %d and %d\n",ttopt[i],ttopt[j]);
	printf("have chosen the same TAE point. I have not\n");
	printf("has a chance to code a fix to this possibility.\n");
	printf("Send mail to the author for a fix.\n");
	printf("\n\t\t\tSorry.");
	exit(1);
				}
			}
		}
	}	
/* No longer needed, now that pairs have been processed */
	if (mto > 0) {
		free(mtopt);
	}
	if (bto > 0) {
		free(btopt);
	}
	if (tto > 0) {
		free(ttopt);
	}
	if (mae > 0) {
		free(maept);
	}
	if (bae > 0) {
		free(baept);
	}
	if (tae > 0) {
		free(taept);
	}
/* set the list for reordering */
	list = (int *) malloc(sizeof(int) * eutpoints);
	tmplist = (int *) malloc(sizeof(int) * eutpoints);
	bufflist = (int *) malloc(sizeof(int) * eutpoints);
	if ((list == NULL) ||
		(tmplist == NULL) ||
		(bufflist == NULL)) {
		puts("Memory allocation error 8");
		exit(1);
		}
	for (i=0; i<eutpoints; i++) {
		list[i] = ptdat[i].ptnum;
		tmplist[i] = 0;
		bufflist[i] = 0;
		}
		
/* Reorder the numbering */	
/* Start the MTO */
for (i=0; i<mto; i++) {
	if (pmtop[i].to != -99) {
	
/* verify that the *ae has a lower value that the *to */
	if (pmtop[i].to < pmtop[i].ae) {
	printf("\nThe MTO point should have a higher point number\n");
	printf("than the MAE point. The proximity test has failed.\n");
	printf("\n\t\t\tSorry.");
	exit(1);
	}		
/* copy the numbering list into the templist up to and including the ae */
	test = 0;
	j = 0;
	while (test == 0) {
		tmplist[j] = list[j];
		if (list[j] == pmtop[i].ae) {
			test = 1;
		}
		j ++;
	}
/* take the values after ae and up to but not including *to, put in buffer */
	test = 0;
	numbuff = 0;
	j = (pmtop[i].ae);
	while (test == 0) {
		if (list[j] == pmtop[i].to) {
			test = 1;
		}
		else {
			bufflist[j] = list[j];
			numbuff ++;
			j ++;
		}
	}		
/* figure out how long the tree goes after the to. */
	nerequire = 1;
	test = 0;
	j = (pmtop[i].to - 1);
	k = pmtop[i].ae;
	while (test == 0) {
		if (strncmp(ptdat[j].pttype,"BP",2) == 0) {
			nerequire ++;
			}
		else if ((strncmp(ptdat[j].pttype,"NE",2) == 0) ||
			 (strncmp(ptdat[j].pttype,"MAE",3) == 0) ||
			 (strncmp(ptdat[j].pttype,"BAE",3) == 0) ||
			 (strncmp(ptdat[j].pttype,"TAE",3) == 0)){
			nerequire --;
			}
		tmplist[k] = list[j];
		if (nerequire == 0) {
			test = 1;
			lastne = list[j];
			}	
		k ++;
		j ++;		
		}
/* write the buffer values into the templist */
	j = pmtop[i].ae;
	while (j < pmtop[i].to -1 ){
		tmplist[k] = bufflist[j];
		j ++;
		k ++;
		}
/* write the remaining values after the to's ne into the templist */
	for (j = (lastne + 1); j < eutpoints; j++){
		tmplist[j] = list[j];
		}		
/* write 'em all to list */
	for ( i = 0; i < eutpoints; i ++) {
		list[i] = tmplist[i];
	}
}
}
/* start the BTO */
for (i=0; i<bto; i++) {
	if (pbtop[i].to != -99) {
/* verify that the *ae has a lower value that the *to */
	if (pbtop[i].to < pbtop[i].ae) {
	printf("\nThe BTO point should have a higher point number\n");
	printf("than the BAE point. The proximity test has failed.\n");
	printf("\n\t\t\tSorry.");
	exit(1);
	}		
/* copy the numbering list into the templist up to and including the ae */
	test = 0;
	j = 0;
	while (test == 0) {
		tmplist[j] = list[j];
		if (list[j] == pbtop[i].ae) {
			test = 1;
		}
		j ++;
	}
/* take the values after ae and up to but not including *to, put in buffer */
	test = 0;
	numbuff = 0;
	j = (pbtop[i].ae);
	while (test == 0) {
		if (list[j] == pbtop[i].to) {
			test = 1;
			}
		else {
			bufflist[j] = list[j];
			numbuff ++;
			j ++;
		}
	}		
/* figure out how long the tree goes after the to. */
	nerequire = 1;
	test = 0;
	j = (pbtop[i].to - 1);
	k = pbtop[i].ae;
	while (test == 0) {
		if (strncmp(ptdat[j].pttype,"BP",2) == 0) {
			nerequire ++;
			}
		else if ((strncmp(ptdat[j].pttype,"NE",2) == 0) ||
			 (strncmp(ptdat[j].pttype,"MAE",3) == 0) ||
			 (strncmp(ptdat[j].pttype,"BAE",3) == 0) ||
			 (strncmp(ptdat[j].pttype,"TAE",3) == 0)){
			nerequire --;
			}
		tmplist[k] = list[j];
		if (nerequire == 0) {
			test = 1;
			lastne = list[j];
			}	
		k ++;
		j ++;		
		}
/* write the buffer values into the templist */
	j = pbtop[i].ae;
	while (j < pbtop[i].to -1){
		tmplist[k] = bufflist[j];
		j ++;
		k ++;
		}
/* write the remaining values after the to's ne into the templist */
	for (j = (lastne + 1); j < eutpoints; j++){
		tmplist[j] = list[j];
		}		
/* write 'em all to list */
	for ( i = 0; i < eutpoints; i ++) {
		list[i] = tmplist[i];
	}
}
}			

/* start the TTO */
for (i=0; i<tto; i++) {
	if (pttop[i].to != -99) {
/* verify that the *ae has a lower value that the *to */
	if (pttop[i].to < pttop[i].ae) {
	printf("\nThe TTO point should have a higher point number\n");
	printf("than the TAE point. The proximity test has failed.\n");
	printf("\n\t\t\tSorry.");
	exit(1);
	}		
/* copy the numbering list into the templist up to and including the ae */
	test = 0;
	j = 0;
	while (test == 0) {
		tmplist[j] = list[j];
		if (list[j] == pttop[i].ae) {
			test = 1;
		}
			j ++;
	}
/* take the values after ae and up to but not including *to, put in buffer */
	test = 0;
	numbuff = 0;
	j = (pttop[i].ae);
	while (test == 0) {
		if (list[j] == pttop[i].to) {
			test = 1;
		}
		else {
			bufflist[j] = list[j];
			numbuff ++;
			j ++;
		}
	}		
/* figure out how long the tree goes after the to. */
	nerequire = 1;
	test = 0;
	j = (pttop[i].to - 1);
	k = pttop[i].ae;
	while (test == 0) {
		if (strncmp(ptdat[j].pttype,"BP",2) == 0) {
			nerequire ++;
			}
		else if ((strncmp(ptdat[j].pttype,"NE",2) == 0) ||
			 (strncmp(ptdat[j].pttype,"MAE",3) == 0) ||
			 (strncmp(ptdat[j].pttype,"BAE",3) == 0) ||
			 (strncmp(ptdat[j].pttype,"TAE",3) == 0)){
			nerequire --;
			}
		tmplist[k] = list[j];
		if (nerequire == 0) {
			test = 1;
			lastne = list[j];
			}	
		k ++;
		j ++;		
		}
/* write the buffer values into the templist */
	j = pttop[i].ae;
	while (j < pttop[i].to -1){
		tmplist[k] = bufflist[j];
		j ++;
		k ++;
		}
/* write the remaining values after the to's ne into the templist */
	for (j = (lastne + 1); j < eutpoints; j++){
		tmplist[j] = list[j];
		}		
/* write 'em all to list */
	for ( i = 0; i < eutpoints; i ++) {
		list[i] = tmplist[i];
	}
}
}			

/* Free up space not needed */	
	if (mto > 0) {
	free(pmtop);
	}
        if (bto > 0) {
	free(pbtop);
	}
	if (tto > 0) {
	free(pttop);
	}
	free(tmplist);
	free(bufflist);

/* write it all back to the original data structure */
starr = (struct shortterm *) malloc(sizeof(struct shortterm) * eutpoints);
	if (starr == NULL)
		{
		puts("Memory allocation error 9");
		exit(1);
		}
	for (i = 0; i < eutpoints; i++) {
		starr[i].ptnum = ptdat[list[i]-1].ptnum;
		strcpy(starr[i].pttype,ptdat[list[i]-1].pttype);
		starr[i].pos[0] = ptdat[list[i]-1].pos[0];
		starr[i].pos[1] = ptdat[list[i]-1].pos[1];
		starr[i].pos[2] = ptdat[list[i]-1].pos[2];
		starr[i].rad = ptdat[list[i]-1].rad;
	}
	for (i = 0; i < eutpoints; i++) {
		ptdat[i].ptnum = starr[i].ptnum;
		strcpy(ptdat[i].pttype,starr[i].pttype);
		ptdat[i].pos[0] = starr[i].pos[0];
		ptdat[i].pos[1] = starr[i].pos[1];
		ptdat[i].pos[2] = starr[i].pos[2];
		ptdat[i].rad = starr[i].rad;
	}

/* Free up space not needed */	
	free(list);	
	free(starr);
	
/* 	Make adjacents CP */
	for (i = 0; i < eutpoints-1; i++) {
		if ((strncmp(ptdat[i].pttype,"MAE",3) == 0) &&
			(strncmp(ptdat[i+1].pttype,"MTO",3) == 0)) {
			strcpy(ptdat[i].pttype,"CP");
			strcpy(ptdat[i+1].pttype,"CP");
		}
		if ((strncmp(ptdat[i].pttype,"BAE",3) == 0) &&
			(strncmp(ptdat[i+1].pttype,"BTO",3) == 0)) {
			strcpy(ptdat[i].pttype,"CP");
			strcpy(ptdat[i+1].pttype,"CP");
		}
		if ((strncmp(ptdat[i].pttype,"TAE",3) == 0) &&
			(strncmp(ptdat[i+1].pttype,"TTO",3) == 0)) {
			strcpy(ptdat[i].pttype,"CP");
			strcpy(ptdat[i+1].pttype,"CP");
		}
	}
		
/* change all remaining AE points to NE */
	for (i = 0; i < eutpoints; i++) {
		if ((strncmp(ptdat[i].pttype,"MAE",3) == 0) ||
			(strncmp(ptdat[i].pttype,"BAE",3) == 0) ||
			(strncmp(ptdat[i].pttype,"TAE",3) == 0)) {
			strcpy(ptdat[i].pttype,"NE");
			}
		}
}

/* ****************************************************************** */

void dupcheck(void)
{
int i, j, dup, *duppt;
double EPS = 1.0e-6;
double displace = 1.0e-6;

/* Get rid of duplicate points */
/* see if any exist */
	dup = 0;
	for (i=0; i<eutpoints-1; i++) {
		if (((fabs(ptdat[i].pos[0] - ptdat[i+1].pos[0])) < EPS) &&
		    ((fabs(ptdat[i].pos[1] - ptdat[i+1].pos[1])) < EPS) &&
		    ((fabs(ptdat[i].pos[2] - ptdat[i+1].pos[2])) < EPS)) {
		    dup ++;
		}
	}

/* if they exist, record the point number for alteration */ 
	if (dup > 0) {
		duppt = (int *) malloc(sizeof(int) * dup);
		if (duppt == NULL){
			puts("Memory allocation error 10");
			exit(1);
			} 
		j = 0;
		for (i=0; i<eutpoints-1; i++) {
		if (((fabs(ptdat[i].pos[0] - ptdat[i+1].pos[0])) < EPS) &&
		    ((fabs(ptdat[i].pos[1] - ptdat[i+1].pos[1])) < EPS) &&
		    ((fabs(ptdat[i].pos[2] - ptdat[i+1].pos[2])) < EPS)) {
		    duppt[j] = ptdat[i].ptnum;
		    j ++;
		    }
		 }
/* Warn the end user */
printf("\nPoints with duplicate spatial coordinates have been \n"); 
printf("entered in the dataset, presumably during the original data \n");
printf("collection. These are points:\n");
	for (i=0; i<dup; i++){
		printf("\t\t%d\n",duppt[i]);
		}
printf("These points will be altered by this program to compensate.\n");
printf("The first point will be displaced by 0.01 micron in all three\n");
printf("coordinates. You may want to examine the file and make your\n");
printf("own alterations to the original datafile.\n");
	}


/* Make the alterations */
	for (i=0; i<dup; i++){
		ptdat[duppt[i]-1].pos[0] += displace;
		ptdat[duppt[i]-1].pos[1] += displace;
		ptdat[duppt[i]-1].pos[2] += displace;
		}
	if (dup > 0) { 
		free(duppt);
		 }
}

/* ******************************************************************** */

void balance(void)
{
int i;


/* 	Count number of soma points
	Count number of primary dendrites (MTO)
	Count number of terminal endings (ES)
	Count branch points (BP) 
*/		

	for (i=0; i<eutpoints; i++) {
		if (strncmp(ptdat[i].pttype,"MTO",3) == 0) {
		    		priden ++;
		    	}
		else if (strncmp(ptdat[i].pttype,"NE",2) == 0) {
		    		 endno ++;
		    	}
		else if (strncmp(ptdat[i].pttype,"BP",2) == 0) {
					branchno ++;
				}
	}
		
/* Check for balancing */
	if (branchno+priden != endno) {
		printf("\nThere is an unrecoverable tree balance error.\n");
		printf("This datafile has been either been altered, or the\n");
		printf("collection of data was flawed. The total number of\n");
		printf("branch points plus the number of primary dendrites\n");
		printf("should equal the number of endpoints.\n");
		printf("That is not the case for this tree.\n\n");
		printf("Branches: %d\n", branchno);
		printf("Primary Dendrites: %d\n", priden);
		printf("Terminal endings: %d\n", endno);
		printf("\n\n\t\tSorry.");
		exit(1);
	}

}

/* ****************************************************************** */

void slicer(void)
{

int j = somapts;
int wbp = 0;
typedef struct bpinfo{
	int bppt;
	int parent;
	int sib;
	} BPINFO;
struct bpinfo *wbppt;	
int branchdone;
int segdone;
double y,tmp1;
double initrad;
int schildct = 0;
int radsplflg = 0;
int bpflag = 0;
int neflag = 0;
int i, reset;
int prompt = 1000;

double x;

/* 
Begin creating compartments. 

A branch is the complete traversal of an MTO. 

A segment is a region between mto/bp, mto/ne, bp/ne, or radius/L slice 
conditions.
	
Note, the ptdat points go away forever in this routine. Look for the slice.

track[0] = iteration constant for the slice points
track[1] = iteration constant for ptdat points
track[2] = interation constant for compartment number
	
The soma is compartment 1, so we start with track[2] = 2.
*/

slice = (struct slicepts *) malloc((sizeof(struct slicepts) * (1)));
        if (slice == NULL){
                puts("Memory allocation error 11");
                exit(1);
        }
	if (branchno > 1) {
	wbppt = (struct bpinfo *) malloc((sizeof(struct bpinfo) * (branchno)));
	}
	else {
	wbppt = (struct bpinfo *) malloc((sizeof(struct bpinfo) * (1)));
	}
	
	if (wbppt == NULL) {
		puts("Memory allocation error 12");
		exit(1);
	}

	somaslicer();
	
/* Initialize */
	wbppt[0].bppt = 0;
	wbppt[0].parent = 0;
	wbppt[0].sib = 0;
		
	track[0] = 1;
	track[1] = j;
	track[2] = 2;
	firstPT = track[1]; 
	lastPT = track[1];

	printf("\nProgress in slicer\n");
	while (j < eutpoints-1) {
		somachild[schildct] = track[2];
		schildct++;
		parent = 1;
		sibling = 0;
		branchdone = 1;
		wbp = 0;
		while (branchdone > 0) {
			segdone = 1;
			y = 0.0;
			initrad = ptdat[j].rad;
			writepoint(ptdat[j].pos,initrad,y);
			track[0]++;
			while(segdone > 0) {
				tmp1 = fabs(ptdat[j+1].rad - initrad);
				if (tmp1 > RadCrit) { 
					if (y > 0.0) {
						lengthslicer();
					}
					else {
						track[0]--;
					}
					radiusslicer();
					radsplflg = 1;
					segdone = -99;
				}
				else { 
			  y += distance(ptdat[j].pos[0],ptdat[j+1].pos[0], 
		 	          	ptdat[j].pos[1],ptdat[j+1].pos[1],
		 		      	ptdat[j].pos[2],ptdat[j+1].pos[2]);
					writepoint(ptdat[j+1].pos,initrad,y); 
				 	track[0]++;
				}
				if (strncmp(ptdat[j+1].pttype,"BP",2) == 0) {
					segdone = -99;
					branchdone++;
					wbp++;
					wbppt[wbp].bppt = j+1;
					bpflag = 1;
				}
				if (strncmp(ptdat[j+1].pttype,"NE",2) == 0) {
					wbp--;
					if (wbp >= 0) {
			NEtmp[0] = ptdat[j+1].pos[0];
			NEtmp[1] = ptdat[j+1].pos[1];
			NEtmp[2] = ptdat[j+1].pos[2];
			NEtmp[3] = ptdat[j+1].rad;
			ptdat[j+1].pos[2] = ptdat[wbppt[wbp+1].bppt].pos[2];
			ptdat[j+1].pos[1] = ptdat[wbppt[wbp+1].bppt].pos[1];
			ptdat[j+1].pos[0] = ptdat[wbppt[wbp+1].bppt].pos[0];
			ptdat[j+1].rad = ptdat[wbppt[wbp+1].bppt].rad;
						neflag = 1;
						netmp = neflag;
					}
					segdone = -99;
					branchdone--;
				}
				if ((segdone == -99) &&
					(radsplflg != 1)) {
					lengthslicer();
				}
			    	if (bpflag == 1){
					wbppt[wbp].parent = track[2]-1;
					wbppt[wbp].sib = track[2];
				} 
				if (neflag == 1) {
					parent = wbppt[wbp+1].parent;
					sibling = wbppt[wbp+1].sib;
				}
				if (strncmp(ptdat[j+2].pttype,"MTO",3) == 0){ 
					j++;
					track[1]++;
					lastPT = track[1];
					branchdone = 0;
				}
				j++;
				track[1]++;
				lastPT = track[1];
				if (segdone == -99) {firstPT = track[1];}
				radsplflg = 0;
				bpflag = 0;
				neflag = 0;
				netmp = neflag;
				if (track[0] > prompt) {
					printf("%d\t", track[0]);  
					prompt += 1000;
				}
			}
		}
	}	
	printf("\n\n");
	free(wbppt);
	free(ptdat);


}

/* ****************************************************************** */
void somaslicer(void)
{
	slice[0].compno = 1;
	slice[0].slicerad = (rangesoma[0]+rangesoma[1]+rangesoma[2])/3.0;
	slice[0].slicelength = (rangesoma[0]+rangesoma[1]+rangesoma[2])/3.0;
	slice[0].slicepos[0] = midpt[0];
	slice[0].slicepos[1] = midpt[1];
	slice[0].slicepos[2] = midpt[2];
	slice[0].parent = 0;
	slice[0].sib = 0;
	slice[0].rchild = 0;
	slice[0].lchild = 0;

	somachild = (int *) malloc(sizeof(int) * priden);
}

/* ****************************************************************** */
double distance(double x1,double x2,double y1,double y2,double z1,double z2)
{
/* The distance formula in 3-D */
	double dist;
	dist = ((x1-x2)*(x1-x2))+((y1-y2)*(y1-y2))+((z1-z2)*(z1-z2));
	return (sqrt(dist));
}

/* ****************************************************************** */
void writepoint(double pos[3], double rad, double y)
{
int i = track[0];
int k = track[2];

        slice = realloc(slice,(sizeof(struct slicepts) * (i+1)));
        if (slice == NULL){
                puts("Memory allocation error 13");
                exit(1);
        }

	slice[i].compno = k;
	slice[i].slicerad = rad;
	slice[i].slicelength = y;
	slice[i].slicepos[0] = pos[0];
	slice[i].slicepos[1] = pos[1];
	slice[i].slicepos[2] = pos[2];
	slice[i].parent = parent;
				
/* Initialize here. These are handled later */ 
	slice[i].lchild = 0;
	slice[i].rchild = 0;
/* half handled here, the other later */
	slice[i].sib = sibling;

}

/* ******************************************************************* */
void conserveSA(void)
{
int i, compno;
double sa, eutd, dist, b, m, d;
double x1, x2, y1, y2, radius;
double firstpt[3], lastpt[3];
double origdist; 

if (CSA == 1) {

	i = track[0]-1;

        lastpt[0] = slice[i].slicepos[0];
        lastpt[1] = slice[i].slicepos[1];
        lastpt[2] = slice[i].slicepos[2];

	compno = slice[i].compno;

        dist = 0.0;
        while (slice[i].compno == slice[i-1].compno) {
        	dist += distance(slice[i].slicepos[0],slice[i-1].slicepos[0],
                             slice[i].slicepos[1],slice[i-1].slicepos[1],
                             slice[i].slicepos[2],slice[i-1].slicepos[2]);
        	i--;
        }
	
	firstpt[0] = slice[i].slicepos[0];
        firstpt[1] = slice[i].slicepos[1];
        firstpt[2] = slice[i].slicepos[2];

	origdist = dist;

	sa = 0.0;

	while (dist > 0.0) {
	
	if (firstPT == lastPT) {

       	if (netmp == 0) {      
               	eutd = distance(ptdat[firstPT].pos[0],ptdat[firstPT+1].pos[0], 
                            ptdat[firstPT].pos[1],ptdat[firstPT+1].pos[1],    
                            ptdat[firstPT].pos[2],ptdat[firstPT+1].pos[2]);   
               	m = (ptdat[firstPT+1].rad - ptdat[firstPT].rad)/eutd;      
		y2 = ptdat[firstPT+1].rad;
       	}
       	if (netmp == 1) {      
               	eutd = distance(ptdat[firstPT].pos[0],NEtmp[0],
                       	ptdat[firstPT].pos[1], NEtmp[1], ptdat[firstPT].pos[2],
                       	NEtmp[2]);     
               	m = (NEtmp[3] - ptdat[firstPT].rad)/eutd; 
		y2 = NEtmp[3];
       	}
	}

	if (firstPT < lastPT) {
                eutd = distance(ptdat[firstPT].pos[0],ptdat[firstPT+1].pos[0],
                            ptdat[firstPT].pos[1],ptdat[firstPT+1].pos[1],
                            ptdat[firstPT].pos[2],ptdat[firstPT+1].pos[2]);
                m = (ptdat[firstPT+1].rad - ptdat[firstPT].rad)/eutd;
		y2 = ptdat[firstPT+1].rad;
	}

	if (eutd != 0.0) {

	b = ptdat[firstPT].rad;

        x1 = distance(ptdat[firstPT].pos[0],firstpt[0], ptdat[firstPT].pos[1],
                        firstpt[1], ptdat[firstPT].pos[2], firstpt[2]);
        
	y1 = (m*x1) + b;

	if (dist < (eutd - x1)) {	/* then ends in this segment */
                x2 = x1 + dist;
                y2 = (m*x2) + b;
        	radius = (y1 + y2)/2.0;
               	sa += (radius * dist);
		dist -= dist;
	}
	if (dist >= (eutd - x1)) { 	/* then cut and repeat until done */	
		radius = (y1 + y2)/2.0;
		sa += (radius * (eutd - x1)); 
		dist -= (eutd - x1);
		firstPT++;
		firstpt[0] = ptdat[firstPT].pos[0];
		firstpt[1] = ptdat[firstPT].pos[1];
		firstpt[2] = ptdat[firstPT].pos[2];
	}
	}
	
	if ((eutd == 0.0) && (dist >= (eutd - x1))) {
	        firstPT++;
                firstpt[0] = ptdat[firstPT].pos[0];
                firstpt[1] = ptdat[firstPT].pos[1];
                firstpt[2] = ptdat[firstPT].pos[2];
	}
	}

/* calculate radius to make compartment SA match */

	i = track[0]-1;

	while (slice[i].compno == slice[i-1].compno) {
		slice[i].slicerad = sa/origdist; 
		slice[i].slicelength = origdist;
		i--;
	}
	slice[i].slicerad = sa/origdist;
	slice[i].slicelength = origdist;
}
}

/* ******************************************************************* */
void radiusslicer(void)
{
double x, y[3], z;
int divmax = 0;
double fdivmax = 0.0;
int divflg = 0;
int h = 0;
typedef struct radsltmp{
		double tmprad;
		double tmppos[3];
		}RADSLTMP;
struct radsltmp *radpt;
int j = track[1];
double eps;
	
/* 
The length between point j and j+1 needs to be divided so that the
radii do not vary by more than a user specified criterium in micron.
	
Note the radius in "slicer" is recorded from the reference radius. (The
radius at the start of the segment.) In this routine, however, the radius
is from the second point of the subdivided segment. It should not matter- 
the whole point of the criterium is to choose a value for which the radius
can vary without numerical consequences. 
*/

/* 	Determine how many subcompartments are needed	*/
	x = (fabs(ptdat[j].rad - ptdat[j+1].rad))/RadCrit;
	eps = (1.0e-6)*x;
	
	while (divflg == 0){
		divmax++;
		fdivmax += 1.0;
		if (fdivmax >= (x+eps)){
			divflg = 1;
		}
	}

/* 	Create a temporary array to store the points in. */
radpt = (struct radsltmp *) malloc(sizeof(struct radsltmp) * (divmax+1));
	if (radpt == NULL) {
		puts("Memory allocation error 14");
		exit(1);
	}

	
/* 	Determine increment of new points	*/
	y[0] = (ptdat[j+1].pos[0] - ptdat[j].pos[0])/(fdivmax);
	y[1] = (ptdat[j+1].pos[1] - ptdat[j].pos[1])/(fdivmax);
	y[2] = (ptdat[j+1].pos[2] - ptdat[j].pos[2])/(fdivmax);	
	x = (ptdat[j+1].rad - ptdat[j].rad)/(fdivmax);

/* dupcheck should have taken care of this.... hmmm. check later */
if ((y[0] == 0.0) && (y[1] == 0.0) && (y[2] == 0.0)) {
	ptdat[j+1].pos[0] += 1.0e-6;
	ptdat[j+1].pos[1] += 1.0e-6;
	ptdat[j+1].pos[2] += 1.0e-6;
	y[0] = (ptdat[j+1].pos[0] - ptdat[j].pos[0])/(fdivmax);
        y[1] = (ptdat[j+1].pos[1] - ptdat[j].pos[1])/(fdivmax);
        y[2] = (ptdat[j+1].pos[2] - ptdat[j].pos[2])/(fdivmax); 
}

/* 	Write new points to temparray */
	radpt[h].tmprad	= ptdat[j].rad;
	radpt[h].tmppos[0] = ptdat[j].pos[0];
	radpt[h].tmppos[1] = ptdat[j].pos[1];
	radpt[h].tmppos[2] = ptdat[j].pos[2];	
	
	for (h = 1; h < divmax+1; h++) {
		radpt[h].tmprad	= radpt[h-1].tmprad + x;
		radpt[h].tmppos[0] = radpt[h-1].tmppos[0] + y[0];
		radpt[h].tmppos[1] = radpt[h-1].tmppos[1] + y[1];
		radpt[h].tmppos[2] = radpt[h-1].tmppos[2] + y[2];
	}
/* write each pair to slice */
	for (h = 0; h < divmax; h++) {
		writepoint(radpt[h].tmppos,radpt[h+1].tmprad,0.0); 
		track[0]++;	
		z = distance(radpt[h].tmppos[0], radpt[h+1].tmppos[0],
				 radpt[h].tmppos[1], radpt[h+1].tmppos[1],
				 radpt[h].tmppos[2], radpt[h+1].tmppos[2]);
		writepoint(radpt[h+1].tmppos,radpt[h+1].tmprad,z); 
		track[0]++;	
		lengthslicer();
	}
	free(radpt);
}

/* ******************************************************************* */

void lengthslicer(void)
{
int loop = track[0]-1;
int segmentlen = 1;
double lambda, Ltest, testloop, newarraylength;
double diffdist, fardist, bipos[3];
double patchdist, inloop, crawl, newdist;
double x, y, z;
double eps;
int i, j, testinc, ininc, last, nullpt;
typedef struct lengthtmp{
		double tmprad;
		double tmplength;
		double tmppos[3];
		}LENGTHTMP;
struct lengthtmp *lengthpt;
double tmpdist;
int wcond = 0;

nullpt = adjerror(loop);

if (nullpt == 0) {

/* How many points are in this segment? */
	while (slice[loop].compno == slice[loop-1].compno) {
		segmentlen++;
		loop--;
	}
	
/* How many points will be needed? */
	lambda = sqrt((Rm*slice[track[0]-1].slicerad)/(2.0*Ri));
	Ltest = slice[track[0]-1].slicelength/lambda;
	testloop = 1.0;
	testinc = 1;
	while (Ltest > LCrit) {
		testloop += 1.0;
		testinc++;
		Ltest = slice[track[0]-1].slicelength/(testloop*lambda);
	}
	newarraylength = slice[track[0]-1].slicelength/testloop;

/* Begin slicing if needed */

	if (testinc > 1) {

/* Allocate space for the temporary array */
lengthpt=(struct lengthtmp *) malloc(sizeof(struct lengthtmp) * (segmentlen));
		if (lengthpt == NULL) {
			puts("Memory allocation error 15");
			exit(1);
		}
		
/* write to the temp array. */
		for (i=0; i< segmentlen; i++) {
			lengthpt[i].tmprad = slice[loop].slicerad;
			lengthpt[i].tmplength = slice[loop].slicelength;
			lengthpt[i].tmppos[0] = slice[loop].slicepos[0];
			lengthpt[i].tmppos[1] = slice[loop].slicepos[1];
			lengthpt[i].tmppos[2] = slice[loop].slicepos[2];
			loop++;
		}

		track[0] -= segmentlen;

/* write the first point. This is a reference point. */
writepoint(lengthpt[0].tmppos,lengthpt[0].tmprad,lengthpt[0].tmplength);
		track[0]++;	

/* Begin the pointwise traversal */
	for (i=1; i< segmentlen; i++) {
		eps = (1.0e-4)*newarraylength;

		if (wcond == 7) {
writepoint(lengthpt[i-1].tmppos,lengthpt[i-1].tmprad,lengthpt[i-1].tmplength);
		track[0]++;
		wcond = 0;
		}
		if (lengthpt[i].tmplength < (newarraylength+eps)) {
writepoint(lengthpt[i].tmppos,lengthpt[i].tmprad,lengthpt[i].tmplength);
		track[0]++; 
		if (i == (segmentlen-1)) {
			parent = track[2];
			track[2]++;
		        if (CSA == 1) { conserveSA(); }
			sibling = 0;
		}
		}
		else {
/* calculate position of intermediary point and write */
		diffdist = newarraylength - lengthpt[i-1].tmplength;
		fardist = lengthpt[i].tmplength - lengthpt[i-1].tmplength;
				bipos[0] = (lengthpt[i-1].tmppos[0] + 
((lengthpt[i].tmppos[0] - lengthpt[i-1].tmppos[0])*(fabs(diffdist/fardist))));
				bipos[1] = (lengthpt[i-1].tmppos[1] + 
((lengthpt[i].tmppos[1] - lengthpt[i-1].tmppos[1])*(fabs(diffdist/fardist))));
				bipos[2] = (lengthpt[i-1].tmppos[2] + 
((lengthpt[i].tmppos[2] - lengthpt[i-1].tmppos[2])*(fabs(diffdist/fardist))));
		writepoint(bipos,lengthpt[i].tmprad,newarraylength);
				track[0]++;
				parent = track[2];
				track[2]++;	
				if (CSA == 1) { conserveSA(); }
				sibling = 0;
				writepoint(bipos,lengthpt[i].tmprad,0.0);
				track[0]++;	
/* displace distance by newarraylength */
	tmpdist = distance(lengthpt[i-1].tmppos[0], bipos[0],
		lengthpt[i-1].tmppos[1],bipos[1],
		lengthpt[i-1].tmppos[2],bipos[2]);
	tmpdist += lengthpt[i-1].tmplength;
				for (j = i; j < segmentlen; j++) {
				lengthpt[j].tmplength -= tmpdist;
				}

/* Recursion until the length at point i is less than the criterion */
			if (lengthpt[i].tmplength < (newarraylength+eps)) {
writepoint(lengthpt[i].tmppos,lengthpt[i].tmprad,lengthpt[i].tmplength);
				track[0]++;
			if (i ==(segmentlen-1)) {
				parent = track[2];
				track[2]++;
				if (CSA == 1) { conserveSA(); }
				sibling = 0;
			}
			}
			else {
/* how far apart are they? */
			patchdist = lengthpt[i].tmplength/newarraylength;
			eps = (1.0e-4)*patchdist;
					inloop = 1.0;
					ininc = 1;
					while ((patchdist+eps) > inloop) {
						inloop += 1.0;
						ininc++;
					}
			x = ((lengthpt[i].tmppos[0] - bipos[0])/inloop);
			y = ((lengthpt[i].tmppos[1] - bipos[1])/inloop);
			z = ((lengthpt[i].tmppos[2] - bipos[2])/inloop);

			newdist = lengthpt[i].tmplength/inloop;

                tmpdist = distance(lengthpt[i].tmppos[0], bipos[0],
                lengthpt[i].tmppos[1],bipos[1],
                lengthpt[i].tmppos[2],bipos[2]);

				for (last = 1; last < ininc + 1; last++) {
						bipos[0] += x;
						bipos[1] += y;
						bipos[2] += z;
				writepoint(bipos,lengthpt[i].tmprad,newdist);
						track[0]++;
						parent = track[2];
						track[2]++;	
						if (CSA == 1) { conserveSA(); }
						sibling = 0;
						wcond = 7;
					if (last != (ininc)) {
				writepoint(bipos,lengthpt[i].tmprad,0.0);
						track[0]++;	
						}
					}

/* reset all the distances */
					for (j = i; j < segmentlen; j++) {
				lengthpt[j].tmplength -= tmpdist;
					}
				}
			}
		}
	}
	
/* if subdivision not needed */
	else {
		parent = track[2];
		track[2]++;
		if (CSA == 1) { conserveSA(); }
		sibling = 0;
		}
if (testinc > 1) {
       	free(lengthpt);
}
}
}

/* ****************************************************************** */
int adjerror(int loop)
{
int i,j, ptct, nullpt;
double dist;

double tmp2;
int tmp;

/* dupcheck can't get those dups at all branch points */

/* first check the slice points just created */

	j = 0;
	ptct = 0;
	nullpt = 0;
	tmp = loop;
        while (slice[loop].compno == slice[loop-1].compno) {
		ptct++;
                dist = distance(slice[loop].slicepos[0],slice[loop-1].slicepos[0],
                             slice[loop].slicepos[1],slice[loop-1].slicepos[1],
                             slice[loop].slicepos[2],slice[loop-1].slicepos[2]);
		if (slice[loop-1].compno != slice[loop-2].compno) {
			j = -99;
		}	
		if ((dist == 0.0) && (j == 0)) {
			for (i = loop-1; i < tmp; i++) {
        			slice[i].compno = slice[i+1].compno;
        			slice[i].slicerad = slice[i+1].slicerad; 
        			slice[i].slicelength = slice[i+1].slicelength;
        			slice[i].slicepos[0] = slice[i+1].slicepos[0]; 
        			slice[i].slicepos[1] = slice[i+1].slicepos[1]; 
        			slice[i].slicepos[2] = slice[i+1].slicepos[2]; 
        			slice[i].parent = slice[i+1].parent; 
        			slice[i].lchild = slice[i+1].lchild; 
        			slice[i].rchild = slice[i+1].rchild;
        			slice[i].sib = slice[i+1].sib; 
			}
			track[0]--;
		}
		if ((dist == 0.0) && (j == -99) && (ptct == 1)) {
			track[0]-2;
			track[2]--;
			nullpt = 1;
                }
		loop--;
	}

	return(nullpt);
}

/* ****************************************************************** */
void topology(void)
{
int prompt, i, reset;

/* Final Topology */
	printf("\nTopology Setting\n\t\t(this may take a bit...)\n");
	prompt = 1000;
	for (i = 0; i < track[0]; i++) {
    	if (slice[i].parent != 0) {
        	for (reset = 1; reset < track[0]; reset++) {
            	if ((slice[reset].compno == slice[i].parent) &&
                    (slice[reset].lchild == 0)) {
                	slice[reset].lchild = slice[i].compno;
                }
                else if ((slice[reset].compno == slice[i].parent) &&
                		 (slice[reset].lchild != 0) &&
                         (slice[i].compno != slice[reset].lchild)) {
                    slice[reset].rchild = slice[i].compno;
                }
            }
        }
    	if (slice[i].sib != 0) {
        	for (reset = 1; reset < track[0]; reset++) {
            	if ((slice[reset].compno == slice[i].sib) &&
                    (slice[reset].sib == 0))    {
                	slice[reset].sib = slice[i].compno;
                }
          	}
  		}
  		if (i > prompt) {
					printf("%d\t", i);
					prompt += 1000;
				}
	} 

	slice[0].parent = 0;
	slice[0].lchild = 0;
	slice[0].rchild = 0;
	slice[0].sib = 0; 
	printf("\n\n");

}

/* ****************************************************************** */
void hinesnum(void)
{
int i;
int hinesno = 1;
int capflg = 0;

/* Allocate Hines memory */

	hines = (int *) malloc(sizeof(int) * (track[2]+1));
	if (hines == NULL) {
		puts("Memory allocation error 16");
		exit(1);
	}

/* Initialize Hines Array */

	for (i = 0; i < track[2]+1; i++) { 
		hines[i] = -99;
	}

	printf("Begin Hines Renumbering\n");
	
/* Start with the NE points */
	for (i = track[0]-1; i > 0; i--) {
		if (( slice[i].rchild == 0) &&
			(slice[i].lchild == 0)){
			hines[slice[i].compno] = hinesno;
			hinesno++;					
			while (slice[i].compno == slice[i-1].compno) {
				i--;
			}
			capflg = 0;
			if (i == 1) {
				capflg = 1;
			}
/* do until quit condition reached */
			while (capflg == 0) {
				if (((slice[i].sib != 0) && 
					(hines[slice[i].sib] != -99)) ||
					(slice[i].sib == 0)) {
					i--;
					hines[slice[i].compno] = hinesno;
					hinesno++;
				while (slice[i].compno == slice[i-1].compno) {
						i--;
					}
					if ((slice[i].parent == 1) ||
					    (i == 0)) {
						capflg = 1;
					}
				}
				else {
					capflg = 1;
				}
			}	
		}
	}	
	
	hines[1] = hinesno;
	hines[0] = 0;
	
/* Renumber everything */
	for (i = 0; i<track[0]; i++) {
		slice[i].compno = hines[slice[i].compno];
		slice[i].parent = hines[slice[i].parent];
		slice[i].sib = hines[slice[i].sib];
		slice[i].lchild = hines[slice[i].lchild];
		slice[i].rchild = hines[slice[i].rchild];
	}
	for (i = 0; i < priden; i++) {
		somachild[i] = hines[somachild[i]];
	}

printf("\nThe total number of compartments is %d. \n",hinesno);
	track[2] = hinesno + 1;
}


/* ****************************************************************** */
void eichwestband(void)
{
int maxindex, minindex;
int startcount, endcount;
int currentwidth, minwidth, maxwidth;
int summer, Rb, low, high, savings;
int *renum, *renumtmp, *bandchild;
int topsub, botsub, diffsub;
int i,j,k,looper1,looper2;
int backup, Eflag, shift, ewnum;
int flg1, flg2, flg3, flg4, flg5;
int st1, st2;
float timesave, persave;
struct slicepts *ew;

if (outtype[5] == 1) {

printf("\nBegin Eichler West & Wilcox renumbering\n");

/* Manuscripts describing this renumbering method can be
obtained from the author. Mailing information in the
intro comments above. */

renum = (int *) malloc(sizeof(int) * (track[2]+1));
if (renum == NULL) {
	puts("Memory allocation error 17");
	exit(1);
}
renumtmp = (int *) malloc(sizeof(int) * (track[2]+1));
if (renumtmp == NULL) {
        puts("Memory allocation error 18");
        exit(1);
}
bandchild = (int *) malloc(sizeof(int) * (priden+1));
if (bandchild == NULL) {
	puts("Memory allocation error 19");
	exit(1);
}

looper1 = 0;
looper2 = 0;
Eflag = 0;
for (i = 0; i <= track[2]; i++) {
	renum[i] = i;
	renumtmp[i] = i;
}

maxindex = 0;
minindex = 0;
bandchild[maxindex] = track[2] - 2;
bandchild[minindex] = track[2] - 2;
currentwidth = 1;
maxwidth = 1;
minwidth = 1;
endcount = track[2] - 2;
startcount = 1;

if (priden > 0) {

	bandchild[maxindex] = 0;

	for (i = 0; i < priden; i++) {
		if (i != priden-1) {
			bandchild[i] = somachild[i] - somachild[i+1];
		}
		else {
			bandchild[i] = somachild[i];
		}
		if (bandchild[i] >= bandchild[maxindex]) {
			maxindex = i;
		}
		if (bandchild[i] <= bandchild[minindex]) {
			minindex = i;
		}
	}	

	currentwidth = track[2] - 1 - bandchild[priden-1];
	maxwidth = track[2] - 1 - bandchild[minindex];
	minwidth = track[2] - 1 - bandchild[maxindex];

	startcount = somachild[maxindex+1] + 1;
	endcount = somachild[maxindex];

	if (maxindex != priden-1) {
		looper2 = 1;	
		j = 1;
		for (i = startcount; i < endcount + 1; i++) {
			renum[j] = i;
			j++;
		}
		for (i = 1; i < startcount; i++) {
			renum[j] = i;
			j++;
		}
		for (i = endcount+1; i <= track[2]; i++) {
			renum[j] = i;
			j++;
		}
	}
}

if ((minwidth < (bandchild[maxindex]-1)) || (priden == 0)) {
	i = 0;
	while ((slice[i].compno != endcount) && (i < track[0])) {
		i++;
	}
	
	while (((slice[i].lchild == 0) || (slice[i].rchild == 0)) && 
	    (Eflag == 0) && (i < track[0])) {
		i++;
		if ((slice[i+1].parent == track[2] - 1) ||
			(i+1 >= track[0])){
				 Eflag = 1;
		}
	}
	st1 = slice[i].compno;

	if (Eflag == 0) {
		if ((maxindex+1) >= priden) {
			shift = 0;
		}
		else {
			shift = somachild[maxindex+1];  /* check */
		}
		if (slice[i].lchild > slice[i].rchild ) {
			topsub = slice[i].lchild - shift;
			botsub = slice[i].rchild - shift;
		}
		else {
			topsub = slice[i].rchild - shift;
			botsub = slice[i].lchild - shift;
		}
		diffsub = topsub - botsub;

		if ((minwidth < botsub) || (minwidth < diffsub)) {
			looper1 = 1;
		}

		while (looper1 == 1) {
			looper1 = 0;

			if (botsub >= diffsub) {
				backup = botsub + shift;
				if (diffsub > minwidth) {
					minwidth = diffsub;
				}
				if (botsub > maxwidth) {
					maxwidth = botsub;
				}
				if (botsub > currentwidth) {
					currentwidth = botsub;
				}
			}
			else {	
				looper2 = 1;
				backup = topsub + shift;
                		if (botsub > minwidth) {
                        		minwidth = botsub;
                		}
                		if (diffsub > maxwidth) {
                        		maxwidth = diffsub;
                		}
                		if (diffsub > currentwidth) {
                        		currentwidth = diffsub;
                		}

        			k = 1;
        			while (renum[k] != st1) {
                			k++;
        			}
				st1 = k;

				j = botsub + 1;
				for (k = 1; k < diffsub+1; k++) {
					renumtmp[k] = renum[j];
					j++;
				}
				j = 1;
				for (k = diffsub+1; k < st1; k++) {
					renumtmp[k] = renum[j];
                			j++;
        			}
        			for (k = st1; k <= track[2]; k++) {
                			renumtmp[k] = renum[k];
        			}
				for (k = 1; k <= track[2]; k++) {
					renum[k] = renumtmp[k];
				}
				shift = shift + botsub;
			}

		i = 0;
       		while ((slice[i].compno != backup) && (i < track[0])) {
               		i++;
       		}
		
       		while (((slice[i].lchild == 0) || (slice[i].rchild == 0))
               		&& (i < track[0])) {
               		i++;
                	if ((slice[i+1].parent == track[2] - 1) ||
                        	(i+1 >= track[0])){
				i = track[0];
			}
       		}
       		st1 = slice[i].compno;
	
		if (i < (track[0] - 2)) {
                	if (slice[i].lchild > slice[i].rchild ) {
                        	topsub = slice[i].lchild - shift;
                        	botsub = slice[i].rchild - shift;
                	}
                	else {
                        	topsub = slice[i].rchild - shift;
                        	botsub = slice[i].lchild - shift;
                	}
                	diffsub = topsub - botsub;
                	if ((minwidth < diffsub) || (minwidth < botsub)) {
                        	looper1 = 1;
                	}
        	}
	}
	}
}

free(renumtmp);
free(bandchild);

if (looper2 == 1) {

ew = (struct slicepts *) malloc((sizeof(struct slicepts)*(track[0]+1)));
        if (ew == NULL){
                puts("Memory allocation error 20");
                exit(1);
        }

        for (i = 0; i < priden; i++) {
		flg1 = 0;
		for (j = 0; j < track[2]; j++) {	
			if ((somachild[i] == renum[j]) && 
				(flg1 == 0)) {
				somachild[i] = j;
				flg1 = 1;
			}
		}
        }

	printf("\t\tRenumbering will take a bit...\n");

	ewnum = 0;
	for (j = 1; j < track[2]; j++) {
		for (i = 0; i<track[0]; i++) {
			if (slice[i].compno == renum[j]) {
				ew[ewnum].compno = slice[i].compno; 
				ew[ewnum].slicerad = slice[i].slicerad; 
				ew[ewnum].slicelength = slice[i].slicelength; 
				ew[ewnum].slicepos[0] = slice[i].slicepos[0]; 
				ew[ewnum].slicepos[1] = slice[i].slicepos[1]; 
				ew[ewnum].slicepos[2] = slice[i].slicepos[2]; 
				ew[ewnum].parent = slice[i].parent; 
                                ew[ewnum].sib = slice[i].sib; 
                                ew[ewnum].lchild = slice[i].lchild; 
                                ew[ewnum].rchild = slice[i].rchild; 
				ewnum++;
			}
		}
	}

	for (i = 0; i<track[0]; i++) {
		flg1 = 0;
		flg2 = 0;
		flg3 = 0;
		flg4 = 0;
		flg5 = 0;
		for (j = 0; j < track[2]; j++) {
			if ((ew[i].compno == renum[j]) &&
				(flg1 == 0)) { 
				ew[i].compno = j;
				flg1 = 1;
			}
			if ((ew[i].parent == renum[j]) &&
				(flg2 == 0)) {
                                ew[i].parent = j;
				flg2 = 1;
                        }
			if ((ew[i].sib == renum[j]) &&
				(flg3 == 0)) { 
                                ew[i].sib = j;
				flg3 = 1;
                        }
			if ((ew[i].lchild == renum[j]) &&
				(flg4 == 0)) { 
                                ew[i].lchild = j;
				flg4 = 1;
                        }
			if ((ew[i].rchild == renum[j]) &&
				(flg5 == 0)) {
                                ew[i].rchild = j;
				flg5 = 1;
                        }
		}
	}

	j = track[0]-1;
	for (i = 0; i<track[0]; i++) {
		slice[i].compno = ew[j].compno;
		slice[i].parent = ew[j].parent;
		slice[i].sib = ew[j].sib;
		slice[i].lchild = ew[j].lchild;
		slice[i].rchild = ew[j].rchild;
		slice[i].slicerad = ew[j].slicerad;
		slice[i].slicelength = ew[j].slicelength;
		slice[i].slicepos[0] = ew[j].slicepos[0];
		slice[i].slicepos[1] = ew[j].slicepos[1];
		slice[i].slicepos[2] = ew[j].slicepos[2];
		j--;
	}

	free(ew);

/* Brag */
Rb = maxwidth - minwidth;
low = minwidth;
high = maxwidth;
summer = 0;
for (k = low; k<= high+1; k++){
	summer = summer - k;
}
savings = 2 * (((track[2]-1)*Rb) + summer);

if ((savings < 0) || (minwidth <= 1)) {
printf("\nSuspicious results in the optimization routine.");
printf("\nPlease rerun your file with a slightly different set");
printf("\nof parameters and contact the author for a fix.");
printf("\nThank you.\n\n");
exit(1);
}

/* Comment out all the brag statements. */

/*
printf("\nThe potential savings gained by implementing the Eichler West and Wilcox\n");
printf("renumbering algorithm on this dataset is\n");
printf("\t%d calculations per matrix operation.\n\n",savings);
*/

Rb = currentwidth - minwidth;
low = minwidth;
high = currentwidth;
summer = 0;
for (k = low; k<= high+1; k++){
	summer = summer - k;
}
savings = 2 * (((track[2]-1)*Rb) + summer);

/*
printf("The actual savings gained by implementing the Eichler West and Wilcox\n");
printf("renumbering algorithm on this dataset is\n");
printf("\t%d calculations per matrix operation.\n\n",savings);
printf("Any difference between these numbers is due to a somewhat\n");
printf("optimal ordering of points during the original data collection.\n");
printf("\t(Things could have been worse!)\n");
*/

summer = track[2]-1;
for (k = track[2]-2; k>= minwidth; k--){
        summer = summer + (2*k);
}

timesave = savings/60000.0;
summer = summer+savings;
persave = savings*100.0/summer;

/*
if (timesave > 0.5) {
printf("\nFor one second of simulated neuron time using a fixed step of\n");
printf("100 microseconds running on a machine capable of 10 Mflops,\n");
printf("the actual savings represents cutting approx. %0.1f minutes\n",timesave);
printf("per simulation, or approximatley %0.1f percent of the time\n", persave);
printf("spent in matrix-vector calculations.\n");
}
*/
 
}

printf("\nThe banding of this matrix is %d.\n", minwidth);

free(renum);

}
}

/* ****************************************************************** */
void output(void)
{
FILE *CLfile, *LFSfile, *DFSfile;
FILE *grafile, *A, *NEURON;
int i,j,Cno;
int *parent, *genparent;
double *length, *radius, *LFS, *DFS;
double L, tempdist;
double Ltotal = 0.0;
double Lmin = 10.0;
double Lmax = 0.0;
double a, srad, dist;
int tempi1, tempi2;

printf("\nBegin writing output files...\n");

length = (double *) malloc(sizeof(double) * (track[2]+1));
radius = (double *) malloc(sizeof(double) * (track[2]+1));
parent = (int *) malloc(sizeof(int) * (track[2]+1));
if ((length == NULL) ||
	(radius == NULL) ||
	(parent == NULL)) {
        puts("Memory allocation error 21");
        exit(1);
}

	dist = 0.0;
	tempi1 = track[0]-1;
for (i = track[0]-1; i > 0; i--) {
        while (slice[i].compno == slice[i-1].compno) {
                dist += distance(slice[i].slicepos[0],slice[i-1].slicepos[0],
                             slice[i].slicepos[1],slice[i-1].slicepos[1],
                             slice[i].slicepos[2],slice[i-1].slicepos[2]);
        i--;
        }
	tempi2 = i;
	radius[slice[i].compno] = slice[i].slicerad;
	for (j = tempi1; j <= tempi2; j++) {
		slice[j].slicerad = radius[slice[i].compno];
	}
	tempi1 = i;

	length[slice[i].compno] = slice[i].slicelength;	
	if (CSA != 1) { length[slice[i].compno] = dist; }
	parent[slice[i].compno] = slice[i].parent;

        if (dist == 0.0) {
        printf("Warning! Length of zero!\n");
        printf("This has not yet been fixed.\n");
        printf("Run your data again with slightly different parameters\n");
        printf("to correct this problem. Please contact the author\n");
        printf("for more information.\n");
        }
}

parent[track[2]-1] = slice[0].parent;
radius[track[2]-1] = slice[0].slicerad;
length[track[2]-1] = slice[0].slicerad;
 

if (outtype[0] == 1) { custom_out(radius, length); } 	/* custom code */
else if 
   (outtype[0] == 2) { genesis_out(radius, length); }			/* genesis */	


/* Format for Neuron */

else if (outtype[0] == 3) {
	NEURON = fopen("neuron-script.dat","w");

	fprintf(NEURON,"// Script file for the NEURON simulator\n");
	fprintf(NEURON,"// produced by the Mesh Generator for Neurons\n");
	fprintf(NEURON,"// copyright @ 1995-7 Rogene M. Eichler West\n");
	fprintf(NEURON,"// \n");
	fprintf(NEURON,"// Original data by: <........>\n");
	fprintf(NEURON,"// Description of data: <........>\n");
	fprintf(NEURON,"\n");

	fprintf(NEURON,"// The compartments were created with the following:\n");
	fprintf(NEURON,"forall Ra = %f\n",Ri);
	fprintf(NEURON,"forall g_pas = %f\n",1.0/Rm);
	fprintf(NEURON,"\n");
 
	j = track[2]-1;
	fprintf(NEURON,"ndend = %d\n", track[2] - 2);
	fprintf(NEURON,"\n");
	fprintf(NEURON,"create soma, dend[ndend]\n");
	fprintf(NEURON,"\n");
	fprintf(NEURON,"soma { \n");
	fprintf(NEURON,"nseg = 1\n");
	fprintf(NEURON,"L = %lf\n", 1000.0*length[j]);
	fprintf(NEURON,"diam = %lf\n",2000.0*radius[j]);
	fprintf(NEURON,"insert pas\n");
	fprintf(NEURON,"}\n");
	fprintf(NEURON,"\n");

	j = track[2]-2;
	for (i = 1; i < track[0] - 1; i++) {
		while (slice[i].compno == slice[i+1].compno) {
                        i++;
                }
		fprintf(NEURON,"dend[%d] {\n", j-1);
		fprintf(NEURON,"nseg = 1\n");
        	fprintf(NEURON,"diam = %lf\n",2000.0*radius[j]);
		fprintf(NEURON,"L = %lf\n", 1000.0*length[j]);
        	fprintf(NEURON,"insert pas\n");
        	fprintf(NEURON,"}\n");
		fprintf(NEURON,"\n");
		j--;
	}
	
	j = track[2]-2;	
	for (i = 1; i < track[0] - 1; i++) {
		while (slice[i].compno == slice[i+1].compno) {
                        i++;
                }
		if (slice[i].parent != track[2]-1) {
			fprintf(NEURON,"connect dend[%d](0), dend[%d](1)\n",j-1, 
				slice[i].parent - 1);
		}
		else {
			fprintf(NEURON,"connect dend[%d](0), soma(1)\n",j-1);
		}
		j--;
	}

	fclose(NEURON);	
}

/* Graphical coordinates for GL routines */

if (outtype[1] == 1) {
grafile = fopen("graph.dat","w");
if (grafile == NULL) {
	printf("Data file could not be written to.\n");
        exit(1);
}
/* first line needs to be the total number of cylinder objects */ 

fprintf(grafile, "%d\n", track[0]-1);
fprintf(grafile, "%d\n", track[2]-2);

srad = (rangesoma[0]+rangesoma[1]+rangesoma[2])/3.0;
fprintf(grafile, "%f %f %f %f\n", srad,midpt[0],midpt[1],midpt[2]); 

for (i = track[0]-1; i > 0; i--) {
fprintf(grafile, "%d   %f\t", slice[i].compno, slice[i].slicepos[0]);
fprintf(grafile, "%f   %f\t", slice[i].slicepos[1], slice[i].slicepos[2]);
fprintf(grafile, "%f\n", slice[i].slicerad);
}
fclose(grafile);
}

/* Misc files for electrotonic length */ 
 
if ((outtype[2] == 1) ||
	(outtype[3] == 1)) {
	if (outtype[2] == 1) {
        	CLfile = fopen("compL.dat","w");
        	if (CLfile == NULL) {
                	printf("Data file could not be written to.\n");
                	exit(1);
       	 	}
	}
	if (outtype[3] == 1) {
        	LFSfile = fopen("LfS.dat","w");
        	if (LFSfile == NULL) {
                	printf("Data file could not be written to.\n");
                	exit(1);
        	}
		LFS = (double *) malloc(sizeof(double) * (track[2]+1));
        	if (LFS == NULL) { 
                	puts("Memory allocation error 22");
                	exit(1);
		}
	}
	for (i = 1; i<track[2]; i++) {
	        L = length[i]/(sqrt((Rm*radius[i])/(2.0*Ri)));
       		Ltotal += L;
	if (i != (track[2]-1)) {
        	if (L > Lmax) {
               		 Lmax = L;
			 Cno = i;
        	}
        	if (L < Lmin) {
                	Lmin = L;
        	}
	}
		if (outtype[2] == 1) {
			fprintf(CLfile,"%lf\n",L);
		}
		if (outtype[3] == 1) {
		LFS[i] = L;
		}
	}
	if (outtype[2] == 1) {
		fclose(CLfile);
	}
printf("\nAverage electrotonic length %f\t", Ltotal/(track[2]-1));
printf("(includes the soma compartment)\n");
printf("Max electrotonic length %f\n", Lmax);
if (Cno == (track[2]-1)) {
printf("\t\t(at soma compartment)\n");
}
printf("Min electrotonic length %f\n", Lmin);
	if (outtype[3] == 1) {
		LFS[0] = 0.0;
        	Lmax = 0.0;
		for (i = track[2]-1; i>0; i--) {
		LFS[i] = LFS[i] + LFS[parent[i]];
		if (LFS[i] > Lmax) {
                        Lmax = LFS[i];
                        Cno = i;
                }
		}
	        for (i = 1; i<track[2]; i++) {
			fprintf(LFSfile,"%lf\n",LFS[i]);
		}		
		printf("\nMax L from Soma %lf at comp %d \n", Lmax, Cno);
		fclose(LFSfile);
	}
}

/* Physical distance from the soma in cm */

if (outtype[4] == 1) {
        DFSfile = fopen("DfS.dat","w");
        if (DFSfile == NULL) {
                printf("Data file could not be written to.\n");
                exit(1);
        }
	DFS = (double *) malloc(sizeof(double) * (track[2]+1));
                if (DFS == NULL) {
                        puts("Memory allocation error 23");
                        exit(1);
                }

	DFS[0] = 0.0;
	Lmax = 0.0;
	for (i = track[2]-1; i>0; i--) {
                DFS[i] = length[i] + DFS[parent[i]];
		if (DFS[i] > Lmax) {
			Lmax = DFS[i];
                	Cno = i;
		}
         }
         for (i = 1; i<track[2]; i++) {
                        fprintf(DFSfile,"%lf\n",DFS[i]);
         }

	printf("\nMax distance from Soma %lf cm at comp %d\n", Lmax, Cno);

	fclose(DFSfile);
}

/* Topography file. Uncomment for debugging or for graph structures */

#ifdef notdef

A = fopen("topo.dat","w");
        if (A == NULL) {
                printf("Data file could not be written to.\n");
                exit(1);
        }
for (i = track[0]-1; i > 0; i--) {
fprintf(A,"c %d p %d s %d ", slice[i].compno, slice[i].parent,slice[i].sib);
fprintf(A,"lc %d rc %d\n",slice[i].lchild, slice[i].rchild);
}
fclose(A);

#endif

}

void custom_out(double *radius, double *length)
{
FILE *radfile, *lengthfile, *connectfile;
int i, j, prtflg;

        radfile = fopen("radii.dat","w");
        if (radfile == NULL) {
                printf("Data file could not be written to.\n");
                exit(1);
        }
        lengthfile = fopen("length.dat","w");
        if (lengthfile == NULL) {
                printf("Data file could not be written to.\n");
                exit(1);
        }
        connectfile = fopen("connect.dat","w");
        if (connectfile == NULL) {
                printf("Data file could not be written to.\n");
                exit(1);
        }
 
        /* Write  radius, length */
 
        for (i = 1; i<track[2]; i++) {
                fprintf(lengthfile,"%20.19lf\n",length[i]);
                fprintf(radfile,"%11.10lf\n",radius[i]);
        }
        fclose(radfile);
        fclose(lengthfile);
 
 
        /* write connection */
 
        for (i = track[0]-1; i > 0; i--) {
                for (j = 1; j < track[2]; j++) {
                        if ((j == slice[i].parent) ||
                                (j ==  slice[i].rchild) ||
                                (j ==  slice[i].lchild)) {
                                fprintf(connectfile,"1.0\n");
                        }
                        else {
                                fprintf(connectfile,"0.0\n");
                        }
                }
                while (slice[i].compno == slice[i-1].compno) {
                        i--;
                }
        }
 
        /* write somachild array to connection matrix */
 
        for (j = 1; j < track[2]; j++) {
                prtflg = 0;
                for (i = 0; i < priden; i++) {
                        if ( j == somachild[i]) {
                                prtflg = 1;
                        }
                }
                if (prtflg == 1) {
                        fprintf(connectfile,"1.0\n");
                }
                else {
                        fprintf(connectfile,"0.0\n");
                }
        }
        fclose(connectfile);
}

void genesis_out(double *radius, double *length)
{
FILE *GENESIS;
double *gv0, *gv1, *gv2, v0, v1, v2, norm;
int i,j, *genparent;

        GENESIS = fopen("genesis-script.p","w");
 
        fprintf(GENESIS,"// Script file for the GENESIS simulator\n");
        fprintf(GENESIS,"// produced by the Mesh Generator for Neurons\n");
        fprintf(GENESIS,"// copyright @ 1995-7 Rogene M. Eichler West\n");
        fprintf(GENESIS,"// \n");
        fprintf(GENESIS,"// Original data by: <........>\n");
        fprintf(GENESIS,"// Description of data: <........>\n");
        fprintf(GENESIS,"\n");

        gv0 = (double *) malloc(sizeof(double) * (track[2]+1));
        gv1 = (double *) malloc(sizeof(double) * (track[2]+1));
        gv2 = (double *) malloc(sizeof(double) * (track[2]+1));
        genparent = (int *) malloc(sizeof(int) * (track[2]+1));
 
        j = track[2]-2;
        for (i = 1; i < track[0] - 1; i++) {
		v0 = slice[i].slicepos[0];
		v1 = slice[i].slicepos[1];
		v2 = slice[i].slicepos[2];

                while (slice[i].compno == slice[i+1].compno) {
                        i++;
                }
		gv0[j] = v0 - slice[i].slicepos[0];
		gv1[j] = v1 - slice[i].slicepos[1];
		gv2[j] = v2 - slice[i].slicepos[2];
		norm = distance(gv0[j], 0.0, gv1[j], 0.0, gv2[j], 0.0);
		gv0[j] = 10000.0*(gv0[j]/norm)*slice[i].slicelength;
		gv1[j] = 10000.0*(gv1[j]/norm)*slice[i].slicelength;
		gv2[j] = 10000.0*(gv2[j]/norm)*slice[i].slicelength;
                genparent[j] = slice[i].parent;
                j--;
        }

/* Genesis absolute option */
 
if (GENopt == 1) {
 
        fprintf(GENESIS,"*absolute\n");
        fprintf(GENESIS,"\n");
        fprintf(GENESIS,"// These values were used to create the compartments\n");
        fprintf(GENESIS,"*set_global RM %f\n",Rm/10000.0);
        fprintf(GENESIS,"*set_global RA %f\n",Ri/100.0);
        fprintf(GENESIS,"\n");
 
        j = track[2]-1;
        fprintf(GENESIS,"soma   none %0.3f %0.3f %0.3f %0.2f\n",
                10000.0*midpt[0], 10000.0*midpt[1], 10000.0*midpt[2], 20000.0*radius[j]);
 
        fprintf(GENESIS,"\n");
 
        for (j = track[2]-2; j > 0; j--) {
                if (genparent[j] != track[2]-1) {
			gv0[j] += gv0[genparent[j]];
			gv1[j] += gv1[genparent[j]];
			gv2[j] += gv2[genparent[j]];
                        fprintf(GENESIS,"dend_%05d dend_%05d %0.3f %0.3f %0.3f %0.2f\n",
                        j, genparent[j], gv0[j], gv1[j], gv2[j], 20000.0*radius[j]);
                }
                else {
			gv0[j] += 10000.0*midpt[0];
			gv1[j] += 10000.0*midpt[1];
			gv2[j] += 10000.0*midpt[2];
                        fprintf(GENESIS,"dend_%05d soma %0.3f %0.3f %0.3f %0.2f\n",
                        j, gv0[j], gv1[j], gv2[j], 20000.0*radius[j]);
                }
        }
}
 
/* Genesis relative option */
 
else if (GENopt == 2) {
        fprintf(GENESIS,"*relative\n");
        fprintf(GENESIS,"\n");


	fprintf(GENESIS,"*set_compt_param RM     {RMs}\n*set_compt_param RA     {RA}\n*set_compt_param CM     {CM}\n*set_compt_param ELEAK {ELEAK}\n*set_compt_param EREST_ACT {EREST_ACT}\n*compt /library/Purk_soma\n");

/*
        fprintf(GENESIS,"// These values were used to create the compartments\n");
        fprintf(GENESIS,"*set_global RM %f\n",Rm/10000.0);
        fprintf(GENESIS,"*set_global RA %f\n",Ri/100.0);
*/
        fprintf(GENESIS,"\n");
 
        j = track[2]-1;
 
        fprintf(GENESIS,"soma   none %0.3f %0.3f %0.3f %0.2f\n", 0.0, 0.0, 0.0,
                20000.0*radius[j]);
        j--;
        fprintf(GENESIS,"\n");

	fprintf(GENESIS,"*set_compt_param RM     {RMd}\n*compt /library/Purk_maind\n");
 
        for (j = track[2]-2; j > 0; j--) {
                if (genparent[j] != track[2]-1) {
                        fprintf(GENESIS,"dend_%05d dend_%05d %0.3f %0.3f %0.3f %0.2f\n",
                        j, genparent[j], gv0[j], gv1[j], gv2[j], 20000.0*radius[j]);
                }
                else {
                        fprintf(GENESIS,"dend_%05d soma %0.3f %0.3f %0.3f %0.2f\n",
                        j, gv0[j], gv1[j], gv2[j], 20000.0*radius[j]);
                }
        }
    }
    fclose(GENESIS);
}