/* This file is part of Cloudy and is copyright (C) 1978-2004 by Gary J. Ferland.
 * For conditions of distribution and use, see copyright notice in license.txt */
/* linfit do linear least squares fit to arrays */
#include "cddefines.h"
#include "bevington.h"

double agauss(double x, 
	  double averag, 
	  double sigma)
{
	/*
	 * function agauss.f
	 *
	 * source
	 *   Bevington, page 48.
	 *
	 * purpose
	 *   evaluate integral of gaussian probability function
	 *
	 * usage
	 *   result = agauss (x, averag, sigma)
	 *
	 * description of parameters
	 *   x      - limit for integral
	 *   averag - mean of distribution
	 *   sigma  - standard deviation of distribution
	 *   integration range is averag +/- z*sigma
	 *      where z = abs(x-averag)/sigma
	 *
	 * subroutines and function subprograms required
	 *   none
	 *
	 */
	double agauss_v;
	double denom, 
	  sum, 
	  term, 
	  y2, 
	  z;

#	ifdef DEBUG_FUN
	fputs( "<+>agauss()\n", debug_fp );
#	endif
	z = fabs(x-averag)/sigma;
	agauss_v = 0.;
	if( z > 0 )
	{
		term = 0.7071067812*z;
		sum = term;
		y2 = (z*z)/2.;
		denom = 1.;

		/* accumulate sums of terms
		 * */
		while( TRUE )
		{
			denom += 2.;
			term *= y2*2./denom;
			sum += term;
			if( term/sum - 1.e-10 <= 0 )
				break;
		}
		agauss_v = 1.128379167*sum*exp(-y2);
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->agauss()\n", debug_fp );
#	endif
	return( agauss_v );
}

double area(double x[], 
	  double y[], 
	  long int npts, 
	  long int nterms)
{
	/*
	 * function area.f
	 *
	 * source
	 *   Bevington, pages 272-273.
	 *
	 * purpose
	 *   integrate the area beneath a set of data points
	 *
	 * usage
	 *   result = area (x, y, npts, nterms)
	 *
	 * description of parameters
	 *   x	   - array of data points for independent variable
	 *   y	   - array of data points for dependent variable
	 *   npts   - number of pairs of data points
	 *   nterms - number of terms in fitting polynomial
	 *
	 * subroutines and function subprograms required
	 *   integ (x, y, nterms, i1, x1, x2 sum)
	 *      fits a polynomial with nterms starting at i1
	 *      and integrates area from x1 to x2
	 *
	 */
	long int i, 
	  i1, 
	  idelta, 
	  imax, 
	  j, 
	  neven;
	double area_v, 
	  x1, 
	  x2;
	double sum;

#	ifdef DEBUG_FUN
	fputs( "<+>area()\n", debug_fp );
#	endif
	sum = 0.;
	if( npts <= nterms )
	{

		/* fit all points with one curve
		 * */
		x1 = x[0];
		x2 = x[npts-1];
		integ(x,y,npts,1,x1,x2,&sum);
	}
	else
	{
		neven = 2*(nterms/2);
		idelta = nterms/2 - 1;
		if( nterms > neven )
		{

			/* odd number of terms
			 * */
			x1 = x[0];
			j = nterms - idelta;
			x2 = (x[j-1] + x[j-2])/2.;
			integ(x,y,nterms,1,x1,x2,&sum);
			i1 = npts - nterms + 1;
			j = i1 + idelta;
			x1 = (x[j-1] + x[j])/2.;
			x2 = x[npts-1];
			integ(x,y,nterms,i1,x1,x2,&sum);
			if( i1 > 2 )
			{
				imax = i1 - 1;
				for( i=2; i <= imax; i++ )
				{
					j = i + idelta;
					x1 = (x[j] + x[j-1])/2.;
					x2 = (x[j+1] + x[j])/2.;
					integ(x,y,nterms,i,x1,x2,&sum);
				}
			}
		}
		else
		{

			/* even number of terms
			 * */
			x1 = x[0];
			j = nterms - idelta;
			x2 = x[j-1];
			integ(x,y,nterms,1,x1,x2,&sum);
			i1 = npts - nterms + 1;
			j = i1 + idelta;
			x1 = x[j-1];
			x2 = x[npts-1];
			integ(x,y,nterms,i1,x1,x2,&sum);
			if( i1 > 2 )
			{
				imax = i1 - 1;
				for( i=2; i <= imax; i++ )
				{
					j = i + idelta;
					x1 = x[j-1];
					x2 = x[j];
					integ(x,y,nterms,i,x1,x2,&sum);
				}
			}
		}
	}
	area_v = sum;
	
#	ifdef DEBUG_FUN
	fputs( " <->area()\n", debug_fp );
#	endif
	return( area_v );
}

void chifit(double X[], 
	  double Y[], 
	  double SIGMAY[], 
	  long int NPTS, 
	  long int nterms, 
	  long int MODE, 
	  double A[], 
	  double deltaa[], 
	  double SIGMAA[], 
	  double YFIT[], 
	  double *CHISQR)
{
	/*
	 * SUBROUTINE CHIFIT.F 
	 *
	 * SOURCE
	 *   BEVINGTON, PAGES 228-231.
	 *
	 * PURPOSE
	 *   MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION
	 *      WITH A PARABOLIC EXPANSION OF CHI SQUARE 
	 *
	 * USAGE 
	 *   CALL CHIFIT (X, Y, SIGMAY, NPTS, nterms, MODE, A, deltaa,
	 *      SIGMAA, YFIT, CHISQR)
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   X      - array OF DATA POINTS FOR INDEPENDENT VARIABLE
	 *   Y      - array OF DATA POINTS FOR DEPENDENT VARIABLE
	 *   SIGMAY - array OF STANDARD DEVIATIONS FOR Y DATA POINTS
	 *   NPTS   - NUMBER OF PAIRS OF DATA POINTS
	 *   nterms - NUMBER OF PARAMETERS
	 *   MODE   - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
	 *            +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
	 *             0 (NO WEIGHTING) WEIGHT(I) = 1.
	 *            -1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
	 *   A      - array OF PARAMETERS
	 *   deltaa - array OF INCREMENTS FOR PARAMETERS A
	 *   SIGMAA - array OF STANDARD DEVIATIONS FOR PARAMETERS A
	 *   YFIT   - array OF CALCULATED VALUES OF Y
	 *   CHISQR - REDUCED CHI SQUARE FOR FIT 
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   BevFunctn (X, I, A)
	 *      EVALUATES THE FITTING FUNCTION FOR THE ITH TERM
	 *   fchisq (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 
	 *      EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA
	 *   matinv (array, nterms, det) 
	 *      INVERTS A SYMMETRIC TWO-DIMENSIONAL MATRIX OF DEGREE nterms
	 *      AND CALCULATES ITS DETERMINANT
	 *
	 * COMMENTS
	 *   DIMENSION STATEMENT VALID FOR nterms UP TO 10
	 *
	 */
	long int I, 
	  I_, 
	  j, 
	  J_, 
	  k, 
	  K_, 
	  NFREE;
	double AJ, 
	  AK, 
	  beta[10], 
	  chisq1, 
	  chisq2, 
	  chisq3, 
	  daArray[10], 
	  DELTA, 
	  det, 
	  FREE;
	double alpha[10][10];

#	ifdef DEBUG_FUN
	fputs( "<+>CHIFIT()\n", debug_fp );
#	endif

	NFREE = NPTS - nterms;
	FREE = (double)NFREE;
	if( NFREE > 0 )
	{
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			YFIT[I_] = BevFunctn(X,I,A);
		}
		chisq1 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);

		/* EVALUATE alpha AND beta MATRICES
		 * */
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;

			AJ = A[J_];
			A[J_] = AJ + deltaa[J_];
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			chisq2 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
			alpha[J_][J_] = chisq2 - 2.*chisq1;
			beta[J_] = -chisq2;
			for( k=1; k <= nterms; k++ )
			{
				K_ = k - 1;
				if( k != j )
				{
					if( k > j )
					{
						alpha[K_][J_] = chisq1 - chisq2;

						AK = A[K_];
						A[K_] = AK + deltaa[K_];
						for( I=1; I <= NPTS; I++ )
						{
							I_ = I - 1;
							YFIT[I_] = BevFunctn(X,I,A);
						}
						chisq3 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
						alpha[K_][J_] += chisq3;
						A[K_] = AK;
					}
					else
					{
						alpha[J_][K_] = (alpha[J_][K_] - chisq2)/2.;
						alpha[K_][J_] = alpha[J_][K_];
					}
				}
			}

			A[J_] = AJ - deltaa[J_];
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			chisq3 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
			A[J_] = AJ;
			alpha[J_][J_] = (alpha[J_][J_] + chisq3)/2.;
			beta[J_] = (beta[J_] + chisq3)/4.;
		}

		/* ELIMINATE NEGATIVE CURVATURE
		 * */
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			if( alpha[J_][J_] == 0 )
			{
				alpha[J_][J_] = 0.01;
			}
			else if( alpha[J_][J_] > 0 )
			{
				continue;
			}
			else
			{
				alpha[J_][J_] = -alpha[J_][J_];
			}
			/* Changed from DO 70 */
			for( k=1; k <= nterms; k++ )
			{
				K_ = k - 1;
				/* Changed from 68, 70, 68 */
				if( k != j )
				{
					alpha[K_][J_] = 0.;
					alpha[J_][K_] = 0.;
				}
				/* New continuation statement */
			}
		}

		/* INVERT MATRIX AND EVALUATE PARAMETER INCREMENTS
		 * */
		matinv(alpha,nterms,&det);
		for( j=0; j < nterms; j++ )
		{
			daArray[j] = 0.;
			for( k=0; k < nterms; k++ )
			{
				daArray[j] += beta[k]*alpha[k][j];
			}
			daArray[j] *= 0.2*deltaa[j];
		}

		/* MAKE SURE CHI SQUARE DECREASES
		 * */
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			A[J_] += daArray[J_];
		}
		while( TRUE )
		{
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			chisq2 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
			if( chisq1 >= chisq2 )
				break;
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				daArray[J_] /= 2.;
				A[J_] -= daArray[J_];
			}
		}

		/* INCREMENT PARAMETERS UNTIL CHI SQUARE STARTS TO INCREASE
		 * */
		chisq1 = -DBL_MAX;
		chisq2 = -DBL_MAX;
		chisq3 = -DBL_MAX;
		while( TRUE )
		{
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				A[J_] += daArray[J_];
			}
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			chisq3 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
			if( chisq3 >= chisq2 )
				break;
			chisq1 = chisq2;
			chisq2 = chisq3;
		}

		/* FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS 
		 * */
		DELTA = 1./(1. + (chisq1 - chisq2)/(chisq3 - chisq2)) + 0.5;
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			A[J_] += -DELTA*daArray[J_];
			SIGMAA[J_] = deltaa[J_]*sqrt(FREE*alpha[J_][J_]);
		}
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			YFIT[I_] = BevFunctn(X,I,A);
		}
		*CHISQR = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
		if( chisq2 < *CHISQR )
		{
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				A[J_] += (DELTA - 1.)*daArray[J_];
			}
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			*CHISQR = chisq2;
		}
	}
	else
	{
		*CHISQR = 0.;
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->CHIFIT()\n", debug_fp );
#	endif
	return;
}

void curfit(double X[], 
	  double Y[], 
	  double SIGMAY[], 
	  long int NPTS, 
	  long int nterms, 
	  long int MODE, 
	  double A[], 
	  double deltaa[], 
	  double SIGMAA[], 
	  double *FLAMDA, 
	  double YFIT[], 
	  double *CHISQR)
{
	/*
	 * SUBROUTINE CURFIT.F 
	 *
	 * SOURCE
	 *   BEVINGTON, PAGES 237-239.
	 *
	 * PURPOSE
	 *   MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION
	 *      WITH A LINEARIZATION OF THE FITTING FUNCTION
	 *
	 * USAGE 
	 *   CALL CURFIT (X, Y, SIGMAY, NPTS, nterms, MODE, A, deltaa,
	 *      SIGMAA, FLAMDA, YFIT, CHISQR)
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   X	   - array OF DATA POINTS FOR INDEPENDENT VARIABLE
	 *   Y	   - array OF DATA POINTS FOR DEPENDENT VARIABLE
	 *   SIGMAY - array OF STANDARD DEVIATIONS FOR Y DATA POINTS
	 *   NPTS   - NUMBER OF PAIRS OF DATA POINTS
	 *   nterms - NUMBER OF PARAMETERS
	 *   MODE   - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
	 *	     +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
	 *	      0 (NO WEIGHTING) WEIGHT(I) = 1.
	 *	     -1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
	 *   A	   - array OF PARAMETERS
	 *   deltaa - array OF INCREMENTS FOR PARAMETERS A
	 *   SIGMAA - array OF STANDARD DEVIATIONS FOR PARAMETERS A
	 *   FLAMDA - PROPORTION OF GRADIENT SEARCH INCLUDED
	 *   YFIT   - array OF CALCULATED VALUES OF Y
	 *   CHISQR - REDUCED CHI SQUARE FOR FIT 
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   BevFunctn (X, I, A)
	 *      EVALUATES THE FITTING FUNCTION FOR THE ITH TERM
	 *   fchisq (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 
	 *      EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA
	 *   fderiv (X, I, A, deltaa, nterms, DERIV)
	 *      EVALUATES THE DERIVATIVES OF THE FITTING FUNCTION
	 *      FOR THE ITH TERM WITH RESPECT TO EACH PARAMETER
	 *   matinv (array, nterms, det) 
	 *      INVERTS A SYMMETRIC TWO-DIMENSIONAL MATRIX OF DEGREE nterms
	 *      AND CALCULATES ITS DETERMINANT
	 *
	 * COMMENTS
	 *   DIMENSION STATEMENT VALID FOR nterms UP TO 10
	 *   SET FLAMDA = 0.001 AT BEGINNING OF SEARCH
	 *
	 */
	long int I, 
	  I_, 
	  j, 
	  J_, 
	  k, 
	  K_, 
	  NFREE;
	double alpha[10][10], 
	  B[10], 
	  beta[10], 
	  chisq1, 
	  DERIV[10], 
	  det, 
	  WEIGHT[100];
	double array[10][10];

#	ifdef DEBUG_FUN
	fputs( "<+>CURFIT()\n", debug_fp );
#	endif

	NFREE = NPTS - nterms;
	if( NFREE > 0 )
	{

		/* EVALUATE WEIGHTS
		 * */
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			if( MODE != 0 )
			{
				if( MODE > 0 )
				{
					WEIGHT[I_] = 1./(SIGMAY[I_] * SIGMAY[I_]);
					continue;
				}
				else if( Y[I_] < 0 )
				{
					WEIGHT[I_] = 1./(-Y[I_]);
					continue;
				}
				else if( Y[I_] != 0 )
				{
					WEIGHT[I_] = 1./Y[I_];
					continue;
				}
			}
			WEIGHT[I_] = 1.;
		}

		/* EVALUATE alpha AND beta MATRICES
		 * */
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			beta[J_] = 0.;
			for( k=1; k <= j; k++ )
			{
				K_ = k - 1;
				alpha[K_][J_] = 0.;
			}
		}
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			fderiv(X,I,A,deltaa,nterms,DERIV);
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				beta[J_] += WEIGHT[I_]*(Y[I_] - BevFunctn(X,I,A))*DERIV[J_];
				for( k=1; k <= j; k++ )
				{
					K_ = k - 1;
					alpha[K_][J_] += WEIGHT[I_]*DERIV[J_]*DERIV[K_];
				}
			}
		}
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			for( k=1; k <= j; k++ )
			{
				K_ = k - 1;
				alpha[J_][K_] = alpha[K_][J_];
			}
		}

		/* EVALUATE CHI SQUARE AT STARTING POINT 
		 * */
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			YFIT[I_] = BevFunctn(X,I,A);
		}
		chisq1 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);

		/* INVERT MODIFIED CURVATURE MATRIX TO FIND NEW PARAMETERS
		 * */
		while( TRUE )
		{
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				for( k=1; k <= nterms; k++ )
				{
					K_ = k - 1;
					array[K_][J_] = alpha[K_][J_]/sqrt(alpha[J_][J_]*
					  alpha[K_][K_]);
				}
				array[J_][J_] = 1. + *FLAMDA;
			}
			matinv(array,nterms,&det);
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				B[J_] = A[J_];
				for( k=1; k <= nterms; k++ )
				{
					K_ = k - 1;
					B[J_] += beta[K_]*array[K_][J_]/sqrt(alpha[J_][J_]*
					  alpha[K_][K_]);
				}
			}

			/* IF CHI SQUARE INCREASED, INCREASE FLAMDA AND TRY AGAIN
			 * */
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,B);
			}
			*CHISQR = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
			if( chisq1 >= *CHISQR )
				break;
			*FLAMDA *= 10.;
		}

		/* EVALUATE PARAMETERS AND UNCERTAINTIES 
		 * */
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			A[J_] = B[J_];
			SIGMAA[J_] = sqrt(array[J_][J_]/alpha[J_][J_]);
		}
		*FLAMDA /= 10.;
	}
	else
	{
		*CHISQR = 0.;
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->CURFIT()\n", debug_fp );
#	endif
	return;
}

double determ(double array[][10], 
	  long int NORDER)
{
	/*
	 * FUNCTION determ.F
	 *
	 * SOURCE
	 *   BEVINGTON, PAGE 294.
	 *
	 * PURPOSE
	 *   CALCULATE THE DETERMINANT OF A SQUARE MATRIX
	 *
	 * USAGE 
	 *   det = determ (array, NORDER)
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   array  - MATRIX
	 *   NORDER - ORDER OF DETERMINANT (DEGREE OF MATRIX)
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   NONE
	 *
	 * COMMENTS
	 *   THIS SUBPROGRAM DESTROYS THE INPUT MATRIX array
	 *   DIMENSION STATEMENT VALID FOR NORDER UP TO 10
	 *
	 */
	int lgHIT;
	long int I, 
	  I_, 
	  j, 
	  J_, 
	  k, 
	  K1, 
	  K_;
	double DETERM_v;
	double SAVE;

#	ifdef DEBUG_FUN
	fputs( "<+>determ()\n", debug_fp );
#	endif

	DETERM_v = 1.;
	for( k=1; k <= NORDER; k++ )
	{
		K_ = k - 1;

		/* INTERCHANGE COLUMNS IF DIAGONAL ELEMENT IS ZERO
		 * */
		if( array[K_][K_] == 0 )
		{
			lgHIT = FALSE;
			for( j=k; j <= NORDER; j++ )
			{
				J_ = j - 1;
				if( array[J_][K_] != 0 )
					lgHIT = TRUE;
			}

			if( !lgHIT )
			{
				DETERM_v = 0.;
				
#	ifdef DEBUG_FUN
				fputs( " <->determ()\n", debug_fp );
#	endif
				return( DETERM_v );
			}

			for( I=k; I <= NORDER; I++ )
			{
				I_ = I - 1;
				SAVE = array[j-1][I_];
				array[j-1][I_] = array[K_][I_];
				array[K_][I_] = SAVE;
			}
			DETERM_v = -DETERM_v;
		}

		/* SUBTRACT ROW k FROM LOWER ROWS TO GET DIAGONAL MATRIX 
		 * */
		DETERM_v *= array[K_][K_];
		if( k < NORDER )
		{
			K1 = k + 1;
			for( I=K1; I <= NORDER; I++ )
			{
				I_ = I - 1;
				for( j=K1; j <= NORDER; j++ )
				{
					J_ = j - 1;
					array[J_][I_] += -array[K_][I_]*array[J_][K_]/
					  array[K_][K_];
				}
			}
		}
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->determ()\n", debug_fp );
#	endif
	return( DETERM_v );

}

double factor(long int n)
{
	/*
	 * function factor.f
	 *
	 * source
	 *   Bevington, page 32.
	 *
	 * purpose
	 *   calculates factorial function for integers
	 *
	 * usage
	 *   result = factor (n)
	 *
	 * description of parameters
	 *   n      - integer argument
	 *
	 * subroutines and function subprograms required
	 *   none
	 *
	 */
	long int i;
	double factor_v;
	double fi, 
	  sum;

#	ifdef DEBUG_FUN
	fputs( "<+>factor()\n", debug_fp );
#	endif
	factor_v = 1.;
	if( n > 1 )
	{
		if( n > 10 )
		{

			/* n greater than 10
			 * */
			sum = 0.;
			for( i=11; i <= n; i++ )
			{
				fi = (double)i;
				sum += log(fi);
			}
			factor_v = 3628800.*exp(sum);
		}
		else
		{

			/* n less than 11
			 * */
			for( i=2; i <= n; i++ )
			{
				fi = (double)i;
				factor_v *= fi;
			}
		}
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->factor()\n", debug_fp );
#	endif
	return( factor_v );
}

double fchisq(double Y[], 
	  double SIGMAY[], 
	  long int NPTS, 
	  long int NFREE, 
	  long int MODE, 
	  double YFIT[])
{
	/*
	 * FUNCTION fchisq.F
	 *
	 * SOURCE
	 *   BEVINGTON, PAGE 194.
	 *
	 * PURPOSE
	 *   EVALUATE REDUCED CHI SQUARE FOR FIT OF DATA 
	 *     fchisq = SUM ((Y-YFIT)**2 / SIGMA**2) / NFREE
	 *
	 * USAGE 
	 *   RESULT = fchisq (Y, SIGMAY, NPTS, NFREE, MODE, YFIT)
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   Y	   - array OF DATA POINTS
	 *   SIGMAY - array OF STANDARD DEVIATIONS FOR DATA POINTS
	 *   NPTS   - NUMBER OF DATA POINTS
	 *   NFREE  - NUMBER OF DEGREES OF FREEDOM
	 *   MODE   - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
	 *	     +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
	 *	      0 (NO WEIGHTING) WEIGHT(I) = 1.
	 *	     -1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
	 *   YFIT   - array OF CALCULATED VALUES OF Y
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   NONE
	 *
	 */
	long int I, 
	  I_;
	double FCHISQ_v, 
	  FREE;
	double CHISQ, 
	  WEIGHT;

#	ifdef DEBUG_FUN
	fputs( "<+>fchisq()\n", debug_fp );
#	endif
	CHISQ = 0.;
	if( NFREE > 0 )
	{

		/* ACCUMULATE CHI SQUARE 
		 * */
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			WEIGHT = 1.;
			if( MODE != 0 )
			{
				if( MODE > 0 )
				{
					WEIGHT = 1./(SIGMAY[I_]*SIGMAY[I_]);
				}
				else if( Y[I_] < 0 )
				{
					WEIGHT = 1./(-Y[I_]);
				}
				else if( Y[I_] != 0 )
				{
					WEIGHT = 1./Y[I_];
				}
			}
			CHISQ += WEIGHT*(Y[I_] - YFIT[I_])*(Y[I_] - YFIT[I_]);
		}

		/*	DIVIDE BY NUMBER OF DEGREES OF FREEDOM
		 * */
		FREE = (double)NFREE;
		FCHISQ_v = CHISQ/FREE;
	}
	else
	{
		FCHISQ_v = 0.;
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->fchisq()\n", debug_fp );
#	endif
	return( FCHISQ_v );
}

void fderiv(double x[], 
	  long int i, 
	  double a[], 
	  double deltaa[], 
	  long int nterms, 
	  double deriv[])
{
	/*
	 * subroutine fderiv.f (non-analytical)
	 *
	 * source
	 *   Bevington, page 242.
	 *
	 * purpose
	 *   evaluate derivatives of function for least-squares search
	 *      for arbitrary function given by BevFunctn
	 *
	 * usage
	 *   call fderiv (x, i, a, deltaa, nterms, deriv)
	 *
	 * description of parameters
	 *   x	   - array of data points for independent variable
	 *   i	   - index of data points
	 *   a	   - array of parameters
	 *   deltaa - array of parameter increments
	 *   nterms - number of parameters
	 *   deriv  - derivatives of function
	 *
	 * subroutines and function subprograms required
	 *   BevFunctn (x, i, a)
	 *      evaluates the fitting function for the ith term
	 *
	 */
	long int j, 
	  j_;
	double aj, 
	  delta, 
	  yfit;

#	ifdef DEBUG_FUN
	fputs( "<+>fderiv()\n", debug_fp );
#	endif
	for( j=1; j <= nterms; j++ )
	{
		j_ = j - 1;
		aj = a[j_];
		delta = deltaa[j_];
		a[j_] = aj + delta;
		yfit = BevFunctn(x,i,a);
		a[j_] = aj - delta;
		deriv[j_] = (yfit - BevFunctn(x,i,a))/(2.*delta);
		a[j_] = aj;
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->fderiv()\n", debug_fp );
#	endif
	return;
}

double bev_gamma(double x)
{
	/*
	 *
	 * source
	 *   Bevington, page 126.
	 *
	 * purpose
	 *   calculate the gamma function for integers and half-integers
	 *
	 * usage
	 *   result = bev_gamma (x)
	 *
	 * description of parameters
	 *   x      - integer or half-integer
	 *
	 * subroutines or function subprograms required
	 *   factor (n)
	 *      calculates n factorial for integers
	 *
	 */
	long int i, 
	  n;
	double gamma_v, 
	  xn;
	double fi, 
	  prod, 
	  sum;

#	ifdef DEBUG_FUN
	fputs( "<+>bev_gamma()\n", debug_fp );
#	endif

	/* integerize argument
	 * */
	n = (long)(x - 0.25);
	xn = (double)n;
	if( x - xn - 0.75 <= 0 )
	{

		/* argument is half-integer
		 * */
		prod = 1.77245385;
		if( n > 0 )
		{
			if( n > 10 )
			{
				sum = 0.;
				for( i=11; i <= n; i++ )
				{
					fi = (double)i;
					sum += log(fi-0.5);
				}
				gamma_v = prod*639383.8623*exp(sum);
				
#	ifdef DEBUG_FUN
				fputs( " <->bev_gamma()\n", debug_fp );
#	endif
				return( gamma_v );
			}
			else
			{
				for( i=1; i <= n; i++ )
				{
					fi = (double)i;
					prod *= fi - 0.5;
				}
			}
		}
		gamma_v = prod;
	}
	else
	{

		/* argument is integer
		 * */
		gamma_v = factor(n);
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->bev_gamma()\n", debug_fp );
#	endif
	return( gamma_v );
}

void gradls(double X[], 
	  double Y[], 
	  double SIGMAY[], 
	  long int NPTS, 
	  long int nterms, 
	  long int MODE, 
	  double A[], 
	  double deltaa[], 
	  double YFIT[], 
	  double *CHISQR)
{
	/*
	 * SUBROUTINE GRADLS.F 
	 *
	 * SOURCE
	 *   BEVINGTON, PAGES 220-222.
	 *
	 * PURPOSE
	 *   MAKE A GRADIENT-SEARCH LEAST-SQUARES FIT TO DATA WITH A
	 *      SPECIFIED FUNCTION WHICH IS NOT LINEAR IN COEFFICIENTS
	 *
	 * USAGE 
	 *   CALL GRADLS (X, Y, SIGMAY, NPTS, nterms, MODE, A, deltaa,
	 *      YFIT, CHISQR)
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   X	   - array OF DATA POINTS FOR INDEPENDENT VARIABLE
	 *   Y	   - array OF DATA POINTS FOR DEPENDENT VARIABLE
	 *   SIGMAY - array OF STANDARD DEVIATIONS FOR Y DATA POINTS
	 *   NPTS   - NUMBER OF PAIRS OF DATA POINTS
	 *   nterms - NUMBER OF PARAMETERS
	 *   MODE   - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
	 *	     +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
	 *	      0 (NO WEIGHTING) WEIGHT(I) = 1.
	 *	     -1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
	 *   A	   - array OF PARAMETERS
	 *   deltaa - array OF INCREMENTS FOR PARAMETERS A
	 *   YFIT   - array OF CALCULATED VALUES OF Y
	 *   CHISQR - REDUCED CHI SQUARE FOR FIT 
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   BevFunctn (X, I, A)
	 *      EVALUATES THE FITTING FUNCTION FOR THE ITH TERM
	 *   fchisq (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 
	 *      EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA
	 *
	 * COMMENTS
	 *   DIMENSION STATEMENT VALID FOR nterms UP TO 10
	 *
	 */
	long int I, 
	  I_, 
	  j, 
	  J_, 
	  NFREE;
	double chisq1, 
	  chisq2, 
	  chisq3, 
	  DELTA, 
	  GRAD[10], 
	  SUM;

#	ifdef DEBUG_FUN
	fputs( "<+>GRADLS()\n", debug_fp );
#	endif

	/* EVALUATE CHI SQUARE AT BEGINNING
	 * */
	NFREE = NPTS - nterms;
	if( NFREE > 0 )
	{
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			YFIT[I_] = BevFunctn(X,I,A);
		}
		chisq1 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);

		/* EVALUATE GRADIENT OF CHI SQUARE
		 * */
		SUM = 0.;
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			DELTA = 0.1*deltaa[J_];
			A[J_] += DELTA;
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			A[J_] -= DELTA;
			GRAD[J_] = chisq1 - fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
			SUM += GRAD[J_]*GRAD[J_];
		}
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			GRAD[J_] = deltaa[J_]*GRAD[J_]/sqrt(SUM);
		}

		/* EVALUATE CHI SQUARE AT NEW POINT
		 * */
		while( TRUE )
		{
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				A[J_] += GRAD[J_];
			}
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			chisq2 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);

			/* MAKE SURE CHI SQUARE DECREASES
			 * */
			if( chisq1 > chisq2 )
				break;
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				A[J_] -= GRAD[J_];
				GRAD[J_] /= 2.;
			}
		}

		/* INCREMENT PARAMETERS UNTIL CHI SQUARE STARTS TO INCREASE
		 * */
		chisq1 = -DBL_MAX;
		chisq2 = -DBL_MAX;
		chisq3 = -DBL_MAX;
		while( TRUE )
		{
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				A[J_] += GRAD[J_];
			}
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			chisq3 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
			if( chisq3 >= chisq2 )
				break;
			chisq1 = chisq2;
			chisq2 = chisq3;
		}

		/* FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS 
		 * */
		DELTA = 1./(1. + (chisq1 - chisq2)/(chisq3 - chisq2)) + 0.5;
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			A[J_] += -DELTA*GRAD[J_];
		}
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			YFIT[I_] = BevFunctn(X,I,A);
		}
		*CHISQR = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
		if( chisq2 < *CHISQR )
		{
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				A[J_] += (DELTA - 1.)*GRAD[J_];
			}
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			*CHISQR = chisq2;
		}
	}
	else
	{
		*CHISQR = 0.;
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->GRADLS()\n", debug_fp );
#	endif
	return;
}

void gridls(double X[], 
	  double Y[], 
	  double SIGMAY[], 
	  long int NPTS, 
	  long int nterms, 
	  long int MODE, 
	  double A[], 
	  double deltaa[], 
	  double SIGMAA[], 
	  double YFIT[], 
	  double *CHISQR)
{
	/*
	 * SUBROUTINE GRIDLS.F 
	 *
	 * SOURCE
	 *   BEVINGTON, PAGES 212-213.
	 *
	 * PURPOSE
	 *   MAKE A GRID-SEARCH LEAST-SQUARES FIT TO DATA WITH A SPECIFIED
	 *   FUNCTION WHICH IS NOT LINEAR IN COEFFICIENTS
	 *
	 * USAGE 
	 *   CALL GRIDLS (X, Y, SIGMAY, NPTS, nterms, MODE, A, deltaa,
	 *      SIGMAA, YFIT, CHISQR)
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   X	   - array OF DATA POINTS FOR INDEPENDENT VARIABLE
	 *   Y	   - array OF DATA POINTS FOR DEPENDENT VARIABLE
	 *   SIGMAY - array OF STANDARD DEVIATIONS FOR Y DATA POINTS
	 *   NPTS   - NUMBER OF PAIRS OF DATA POINTS
	 *   nterms - NUMBER OF PARAMETERS
	 *   MODE   - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
	 *	     +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
	 *	      0 (NO WEIGHTING) WEIGHT(I) = 1.
	 *	     -1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
	 *   A	   - array OF PARAMETERS
	 *   deltaa - array OF INCREMENTS FOR PARAMETERS A
	 *   SIGMAA - array OF STANDARD DEVIATIONS FOR PARAMETERS A
	 *   YFIT   - array OF CALCULATED VALUES OF Y
	 *   CHISQR - REDUCED CHI SQUARE FOR FIT 
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   BevFunctn (X, I, A)
	 *      EVALUATES THE FITTING FUNCTION FOR THE ITH TERM
	 *   fchisq (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) 
	 *      EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA
	 *
	 * COMMENTS
	 *   deltaa VALUES ARE MODIFIED BY THE PROGRAM
	 *
	 */
	long int I, 
	  I_, 
	  j, 
	  J_, 
	  NFREE;
	double chisq1, 
	  chisq2, 
	  chisq3, 
	  DELTA, 
	  FN, 
	  FREE, 
	  SAVE;

#	ifdef DEBUG_FUN
	fputs( "<+>GRIDLS()\n", debug_fp );
#	endif

	NFREE = NPTS - nterms;
	FREE = (double)NFREE;
	*CHISQR = 0.;
	if( NFREE > 0 )
	{
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;

			/* EVALUATE CHI SQUARE AT FIRST TWO SEARCH POINTS
			 * */
			for( I=1; I <= NPTS; I++ )
			{
				I_ = I - 1;
				YFIT[I_] = BevFunctn(X,I,A);
			}
			chisq1 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
			chisq2 = -DBL_MAX;
			FN = 0.;
			DELTA = deltaa[J_];
			while( TRUE )
			{
				A[J_] += DELTA;
				for( I=1; I <= NPTS; I++ )
				{
					I_ = I - 1;
					YFIT[I_] = BevFunctn(X,I,A);
				}
				chisq2 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
				if( chisq1 != chisq2 )
					break;
			}
			if( chisq1 <= chisq2 )
			{

				/* REVERSE DIRECTION OF SEARCH IF CHI SQUARE IS INCREASING
				 * */
				DELTA = -DELTA;
				A[J_] += DELTA;
				for( I=1; I <= NPTS; I++ )
				{
					I_ = I - 1;
					YFIT[I_] = BevFunctn(X,I,A);
				}
				SAVE = chisq1;
				chisq1 = chisq2;
				chisq2 = SAVE;
			}

			/* INCREMENT A(j) UNTIL CHI SQUARE INCREASES
			 * */
			chisq3 = -DBL_MAX;
			while( TRUE )
			{
				FN += 1.;
				A[J_] += DELTA;
				for( I=1; I <= NPTS; I++ )
				{
					I_ = I - 1;
					YFIT[I_] = BevFunctn(X,I,A);
				}
				chisq3 = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
				if( chisq3 >= chisq2 )
					break;
				chisq1 = chisq2;
				chisq2 = chisq3;
			}

			/* FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS 
			 * */
			DELTA *= 1./(1. + (chisq1 - chisq2)/(chisq3 - chisq2)) + 
			  0.5;
			A[J_] -= DELTA;
			SIGMAA[J_] = deltaa[J_]*sqrt(2./(FREE*(chisq3-2.*chisq2+
			  chisq1)));
			deltaa[J_] = deltaa[J_]*FN/3.;
		}

		/* EVALUATE FIT AND CHI SQUARE FOR FINAL PARAMETERS
		 * */
		for( I=1; I <= NPTS; I++ )
		{
			I_ = I - 1;
			YFIT[I_] = BevFunctn(X,I,A);
		}
		*CHISQR = fchisq(Y,SIGMAY,NPTS,NFREE,MODE,YFIT);
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->GRIDLS()\n", debug_fp );
#	endif
	return;
}

void integ(double x[], 
	  double y[], 
	  long int nterms, 
	  long int i1, 
	  double x1, 
	  double x2, 
	  double *sum)
{
	/*
	 * subroutine integ.f
	 *
	 * source
	 *   Bevington, page 274.
	 *
	 * purpose
	 *   integrate the area beneath two data points
	 *
	 * usage
	 *   call integ (x, y, nterms, i1, x1, x2, sum)
	 *
	 * description of parameters
	 *   x	   - array of data points for independent variable
	 *   y	   - array of data points for dependent variable
	 *   nterms - number of terms in fitting polymonial
	 *   i1	   - first data point for fitting polynomial
	 *   x1	   - first value of x for integration
	 *   x2	   - final value of x for integration
	 *
	 * subroutines and function subprograms required
	 *   none
	 *
	 * comments
	 *   dimension statement valid for nterms up to 10
	 *
	 */
	long int i, 
	  imid, 
	  j, 
	  j_, 
	  k, 
	  k_;
	double det, 
	  dx1, 
	  dx2;
	double a, 
	  array[10][10], 
	  deltax, 
	  denom, 
	  xjk;

#	ifdef DEBUG_FUN
	fputs( "<+>integ()\n", debug_fp );
#	endif

	/* construct square matrix and invert
	 * */
	for( j=1; j <= nterms; j++ )
	{
		j_ = j - 1;
		i = j + i1 - 1;
		deltax = x[i-1] - x[i1-1];
		xjk = 1.;
		for( k=1; k <= nterms; k++ )
		{
			k_ = k - 1;
			array[k_][j_] = xjk;
			xjk *= deltax;
		}
	}
	matinv(array,nterms,&det);
	if( det == 0 )
	{
		imid = i1 + nterms/2;
		*sum += y[imid-1]*(x2 - x1);
	}
	else
	{

		/* evaluate coefficients and integrate
		 * */
		dx1 = x1 - x[i1-1];
		dx2 = x2 - x[i1-1];
		for( j=1; j <= nterms; j++ )
		{
			j_ = j - 1;
			i = j + i1 - 1;
			a = 0.;
			for( k=1; k <= nterms; k++ )
			{
				k_ = k - 1;
				a += y[i-1]*array[k_][j_];
			}
			denom = (double)j;
			*sum += (a/denom)*(powi(dx2,j) - powi(dx1,j));
		}
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->integ()\n", debug_fp );
#	endif
	return;
}

void interp(double X[], 
	  double Y[], 
	  long int NPTS, 
	  long int *nterms, 
	  double XIN, 
	  double *YOUT)
{
	/*
	 * SUBROUTINE INTERP.F 
	 *
	 * SOURCE
	 *   BEVINGTON, PAGES 266-267.
	 *
	 * PURPOSE
	 *   INTERPOLATE BETWEEN DATA POINTS TO EVALUATE A FUNCTION
	 *
	 * USAGE 
	 *   CALL INTERP (X, Y, NPTS, nterms, XIN, YOUT) 
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   X	   - array OF DATA POINTS FOR INDEPENDENT VARIABLE
	 *   Y	   - array OF DATA POINTS FOR DEPENDENT VARIABLE
	 *   NPTS   - NUMBER OF PAIRS OF DATA POINTS
	 *   nterms - NUMBER OF TERMS IN FITTING POLYNOMIAL
	 *   XIN    - INPUT VALUE OF X
	 *   YOUT   - INTERPOLATED VALUE OF Y
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   NONE
	 *
	 * COMMENTS
	 *   DIMENSION STATEMENT VALID FOR nterms UP TO 10
	 *   VALUE OF nterms MAY BE MODIFIED BY THE PROGRAM
	 *
	 */
	long int I, 
	  I1, 
	  I2, 
	  IMAX, 
	  IX, 
	  IXMAX, 
	  I_, 
	  j, 
	  J_, 
	  k, 
	  K_;
	double DENOM;
	double A[10], 
	  DELTA[10], 
	  DELTAX, 
	  PROD, 
	  SUM;

#	ifdef DEBUG_FUN
	fputs( "<+>INTERP()\n", debug_fp );
#	endif

	/* SEARCH FOR APPROPRIATE VALUE OF X(1)
	 * */
	for( I=1; I <= NPTS; I++ )
	{
		I_ = I - 1;
		if( XIN == X[I_] )
		{
			*YOUT = Y[I_];
			
#	ifdef DEBUG_FUN
			fputs( " <->INTERP()\n", debug_fp );
#	endif
			return;
		}
		if( XIN <= X[I_] )
			goto L_13;
	}
	I1 = NPTS - *nterms + 1;
	goto L_21;
L_13:
	I1 = I - *nterms/2;
	if( I1 <= 0 )
		I1 = 1;
L_21:
	I2 = I1 + *nterms - 1;
	if( NPTS < I2 )
	{
		I2 = NPTS;
		I1 = I2 - *nterms + 1;
		if( I1 <= 0 )
		{
			I1 = 1;
			*nterms = I2 - I1 + 1;
		}
	}

	/* EVALUATE DEVIATIONS DELTA
	 * */
	DENOM = X[I1] - X[I1-1];
	DELTAX = (XIN - X[I1-1])/DENOM;
	for( I=1; I <= *nterms; I++ )
	{
		I_ = I - 1;
		IX = I1 + I - 1;
		DELTA[I_] = (X[IX-1] - X[I1-1])/DENOM;
	}

	/* ACCUMULATE COEFFICIENTS A
	 * */
	A[0] = Y[I1-1];
	for( k=2; k <= *nterms; k++ )
	{
		K_ = k - 1;
		PROD = 1.;
		SUM = 0.;
		IMAX = k - 1;
		IXMAX = I1 + IMAX;
		for( I=1; I <= IMAX; I++ )
		{
			I_ = I - 1;
			j = k - I;
			PROD *= DELTA[K_] - DELTA[j-1];
			SUM += -A[j-1]/PROD;
		}
		A[K_] = SUM + Y[IXMAX-1]/PROD;
	}

	/* ACCUMULATE SUM OF EXPANSION
	 * */
	SUM = A[0];
	for( j=2; j <= *nterms; j++ )
	{
		J_ = j - 1;
		PROD = 1.;
		IMAX = j - 1;
		for( I=1; I <= IMAX; I++ )
		{
			I_ = I - 1;
			PROD *= DELTAX - DELTA[I_];
		}
		SUM += A[J_]*PROD;
	}
	*YOUT = SUM;
	
#	ifdef DEBUG_FUN
	fputs( " <->INTERP()\n", debug_fp );
#	endif
	return;
}

void legfit(double theta[], 
	  double y[], 
	  double sigmay[], 
	  long int npts, 
	  long int norder, 
	  long int neven, 
	  long int mode, 
	  double ftest[], 
	  double yfit[], 
	  double a[], 
	  double sigmaa[], 
	  double b[], 
	  double sigmab[], 
	  double *chisqr)
{
	/*
	 * subroutine legfit.f
	 *
	 * source
	 *   Bevington, pages 155-157.
	 *
	 * purpose
	 *   make a least-squares fit to data with a legendre polynomial
	 *      y = a(1) + a(2)*x + a(3)*(3x**2-1)/2 + ...
	 *        = a(1) * (1. + b(2)*x + b(3)*(3x**2-1)/2 + ... )
	 *      where x = cos(theta)
	 *
	 * usage
	 *   call legfit (theta, y, sigmay, npts, norder, neven, mode, ftest,
	 *      yfit, a, sigmaa, b, sigmab, chisqr)
	 *
	 * description of parameters
	 *   theta  - array of angles (in degrees) of the data points
	 *   y      - array of data points for dependent variable
	 *   sigmay - array of standard deivations for y data points
	 *   npts   - number of pairs of data points
	 *   norder - highest order of polynomial (number of terms - 1)
	 *   neven  - determines odd or even character of polynomial
	 *            +1 fits only to even terms
	 *             0 fits to all terms
	 *            -1 fits only to odd terms (plus constant term)
	 *   mode   - determines mode of weighting least-squares fit
	 *            +1 (instrumental) weight(i) = 1./sigmay(i)**2
	 *             0 (no weighting) weight(i) = 1.
	 *            -1 (statistical)  weight(i) = 1./y(i)
	 *   ftest  - array of values of f(l) for an f test
	 *   yfit   - array of calculated values of y
	 *   a      - array of coefficients of polynomial
	 *   sigmaa - array of standard deviations for coefficients
	 *   b      - array of normalized relative coefficients
	 *   sigmab - array of standard deviations for relative coefficients
	 *   chisqr - reduced chi square for fit
	 *
	 * subroutines and function subprograms required
	 *   matinv (array, nterms, det)
	 *      inverts a symmetric two-dimensional matrix of degree nterms
	 *      and calculates its determinant
	 *
	 * comments
	 *   dimension statement valid for npts up to 100 and order up to 9
	 *   dcos changed to cos in statement 31
	 *
	 */
	long int i, 
	  i_, 
	  j, 
	  j_, 
	  jmax, 
	  k, 
	  k_, 
	  l, 
	  l_, 
	  ncoeff, 
	  nterms;
	double chisq1, 
	  det, 
	  fl, 
	  freedom, 
	  fvalue, 
	  varnce, 
	  weight[100];
	double alpha[10][10], 
	  beta[10], 
	  chisq, 
	  cosine; 
	  
	/* >>chng 01 sep 11, replace array allocation on stack with
	 * MALLOC to avoid bug in gcc 3.0 on Sparc platforms, PvH */
/*  	double p[10][100]; */
	double **p;

#	ifdef DEBUG_FUN
	fputs( "<+>legfit()\n", debug_fp );
#	endif

	if( ( p = (double**)MALLOC((size_t)(10*sizeof(double*))) ) == NULL )
		BadMalloc();
	for( i=0; i < 10; ++i )
	{
		if( ( p[i] = (double*)MALLOC((size_t)(100*sizeof(double)) ) ) == NULL )
			BadMalloc();
	}

	/* accumulate weights and legendre polynomials
	 * */
	nterms = 1;
	ncoeff = 1;
	jmax = norder + 1;
	for( i=1; i <= npts; i++ )
	{
		i_ = i - 1;
		weight[i_] = 1.;
		if( mode != 0 )
		{
			if( mode > 0 )
			{
				weight[i_] = 1./(sigmay[i_]*sigmay[i_]);
			}
			else if( y[i_] < 0 )
			{
				weight[i_] = 1./(-y[i_]);
			}
			else if( y[i_] != 0 )
			{
				weight[i_] = 1./y[i_];
			}
		}
		cosine = cos(0.01745329252*theta[i_]);
		p[0][i_] = 1.;
		p[1][i_] = cosine;
		for( l=2; l <= norder; l++ )
		{
			l_ = l - 1;
			fl = (double)l;
			p[l_+1][i_] = ((2.*fl - 1.)*cosine*p[l_][i_] - (fl - 1.)*
			  p[l_-1][i_])/fl;
		}
	}

	/* accumulate matrices alpha and beta
	 * */
	while( TRUE )
	{
		for( j=1; j <= nterms; j++ )
		{
			j_ = j - 1;
			beta[j_] = 0.;
			for( k=1; k <= nterms; k++ )
			{
				k_ = k - 1;
				alpha[k_][j_] = 0.;
			}
		}
		for( i=1; i <= npts; i++ )
		{
			i_ = i - 1;
			for( j=1; j <= nterms; j++ )
			{
				j_ = j - 1;
				beta[j_] += p[j_][i_]*y[i_]*weight[i_];
				for( k=j; k <= nterms; k++ )
				{
					k_ = k - 1;
					alpha[k_][j_] += p[j_][i_]*p[k_][i_]*weight[i_];
					alpha[j_][k_] = alpha[k_][j_];
				}
			}
		}

		/* delete fixed coefficients
		 * */
		if( neven != 0 )
		{
			if( neven > 0 )
			{
				for( j=2; j <= nterms; j += 2 )
				{
					j_ = j - 1;
					beta[j_] = 0.;
					for( k=1; k <= nterms; k++ )
					{
						k_ = k - 1;
						alpha[k_][j_] = 0.;
						alpha[j_][k_] = 0.;
					}
					alpha[j_][j_] = 1.;
				}
			}
			else
			{
				for( j=3; j <= nterms; j += 2 )
				{
					j_ = j - 1;
					beta[j_] = 0.;
					for( k=1; k <= nterms; k++ )
					{
						k_ = k - 1;
						alpha[k_][j_] = 0.;
						alpha[j_][k_] = 0.;
					}
					alpha[j_][j_] = 1.;
				}
			}
		}

		/* invert curvature matrix alpha
		 * */
		for( j=1; j <= jmax; j++ )
		{
			j_ = j - 1;
			a[j_] = 0.;
			sigmaa[j_] = 0.;
			b[j_] = 0.;
			sigmab[j_] = 0.;
		}
		for( i=1; i <= npts; i++ )
		{
			i_ = i - 1;
			yfit[i_] = 0.;
		}
		matinv(alpha,nterms,&det);
		if( det == 0. )
		{
			*chisqr = 0.;
			goto L_103;
		}

		/* calculate coefficients, fit, and chi square
		 * */
		for( j=1; j <= nterms; j++ )
		{
			j_ = j - 1;
			for( k=1; k <= nterms; k++ )
			{
				k_ = k - 1;
				a[j_] += beta[k_]*alpha[k_][j_];
			}
			for( i=1; i <= npts; i++ )
			{
				i_ = i - 1;
				yfit[i_] += a[j_]*p[j_][i_];
			}
		}
		chisq = 0.;
		for( i=1; i <= npts; i++ )
		{
			i_ = i - 1;
			chisq += (y[i_] - yfit[i_])*(y[i_] - yfit[i_])*weight[i_];
		}
		freedom = (double)(npts - ncoeff);
		*chisqr = chisq/freedom;

		/* test for end of fit
		 * */
		if( nterms >= jmax )
			break;
		/* gary add following to stop flag used before defined */
		chisq1 = -DBL_MAX;

		if( ncoeff != 2 )
		{
			if( ncoeff > 2 )
			{
				fvalue = (chisq1 - chisq)/ *chisqr;
				if( ftest[nterms-1] >= fvalue )
				{
					if( neven == 0 )
					{
						nterms -= 1;
					}
					else
					{
						nterms -= 2;
					}
					ncoeff -= 1;
					jmax = nterms;
					continue;
				}
			}
			else if( neven > 0 )
			{
				goto L_135;
			}
			else
			{
				goto L_137;
			}
		}
		if( neven != 0 )
			goto L_135;
L_137:
		nterms += 1;
		goto L_138;
L_135:
		nterms += 2;
L_138:
		ncoeff += 1;
		chisq1 = chisq;
	}

	/* calculate remainder of output
	 * */
	if( mode == 0 )
	{
		varnce = *chisqr;
	}
	else
	{
		varnce = 1.;
	}
	for( j=1; j <= nterms; j++ )
	{
		j_ = j - 1;
		sigmaa[j_] = sqrt(varnce*alpha[j_][j_]);
	}
	if( a[0] != 0 )
	{
		for( j=2; j <= nterms; j++ )
		{
			j_ = j - 1;
			if( a[j_] != 0 )
			{
				b[j_] = a[j_]/a[0];
				sigmab[j_] = 
					b[j_]*sqrt(
					(sigmaa[j_]/a[j_]*sigmaa[j_]/a[j_])+
					(sigmaa[0]/a[0]*sigmaa[0]/a[0])-
					2.*varnce*alpha[0][j_]/(a[j_]*a[0]));
			}
		}
		b[0] = 1.;
	}
	
L_103:
	for( i=0; i < 10; ++i )
	{
		free( p[i] );
	}
	free( p );

#	ifdef DEBUG_FUN
	fputs( " <->legfit()\n", debug_fp );
#	endif
	return;
}

/* linfit do linear least squares fit to arrays, ret 0 if ok, 1 for disaster */
int linfit(double x[], 
	  double y[], 
	  double sigmay[], 
	  long int npts, 
	  long int mode, 
	  double *yintcp, 
	  double *sigmaa, 
	  double *slope, 
	  double *sigmab, 
	  double *r)
{
	/* return error condition */
	int ierr;
	/*
	 * subroutine linfit.f
	 *
	 * source
	 *   Bevington, pages 104-105.
	 *
	 * purpose
	 *   make a least-squares fit to a data set with a straight line
	 *
	 * usage
	 *   call linfit (x, y, sigmay, npts, mode, a, sigmaa, b, sigmab, r)
	 *
	 * description of parameters
	 *   x      - array of data points for independent variable
	 *   y      - array of data points for dependent variable
	 *   sigmay - array of standard deviations for y data points
	 *   npts   - number of pairs of data points
	 *   mode   - determines method of weighting least-squares fit
	 *            +1 (instrumental) weight(i) = 1./sigmay(i)**2
	 *             0 (no weighting) weight(i) = 1.
	 *            -1 (statistical)  weight(i) = 1./y(i)
	 *   yintcp  - y intercept of fitted straight line
	 *   sigmaa - standard deviation of a
	 *   slope  - slope of fitted straight line
	 *   sigmab - standard deviation of b
	 *   r      - linear correlation coefficient
	 *
	 * subroutines and function subprograms required
	 *   none
	 *
	 */
	long int i;
	double c;
	double delta, 
	  sum, 
	  sumx, 
	  sumx2, 
	  sumxy, 
	  sumy, 
	  sumy2, 
	  varnce, 
	  weight, 
	  xi, 
	  yi;

#	ifdef DEBUG_FUN
	fputs( "<+>linfit()\n", debug_fp );
#	endif

	/* accumulate weighted sums
	 * */
	sum = 0.;
	sumx = 0.;
	sumy = 0.;
	sumx2 = 0.;
	sumxy = 0.;
	sumy2 = 0.;
	for( i=0; i < npts; ++i )
	{
		xi = x[i];
		yi = y[i];
		weight = 1.;
		if( mode != 0 )
		{
			if( mode > 0 )
			{
				weight = 1./(sigmay[i]*sigmay[i]);
			}
			else if( yi < 0 )
			{
				weight = 1./(-yi);
			}
			else if( yi != 0 )
			{
				weight = 1./yi;
			}
		}
		sum += weight;
		sumx += weight*xi;
		sumy += weight*yi;
		sumx2 += weight*xi*xi;
		sumxy += weight*xi*yi;
		sumy2 += weight*yi*yi;
	}

	/* calculate coefficients and standard deviations
	 * */
	delta = sum*sumx2 - sumx*sumx;
	if( fabs( delta )>SMALLFLOAT )
	{
		/* the y intercept */
		*yintcp = (sumx2*sumy - sumx*sumxy)/delta;
		/* the slope */
		*slope = (sumxy*sum - sumx*sumy)/delta;
		if( mode == 0 )
		{
			c = (double)(npts - 2);
			varnce = (sumy2 + *yintcp**yintcp*sum + *slope**slope*sumx2 - 2.*(*yintcp*sumy + 
			*slope*sumxy - *yintcp**slope*sumx))/c;
		}
		else
		{
			varnce = 1.;
		}
		/* variance can be slightly negative when points very close together */
		*sigmaa = sqrt(MAX2(0.,varnce*sumx2/delta));
		*sigmab = sqrt(MAX2(0.,varnce*sum/delta));
		delta = MAX2(SMALLFLOAT,delta*(sum*sumy2-sumy*sumy));
		*r = (sum*sumxy - sumx*sumy)/sqrt(delta);
		ierr = 0;
	}
	else
		ierr = 1;
	
#	ifdef DEBUG_FUN
	fputs( " <->linfit()\n", debug_fp );
#	endif
	return ierr;
}

void matinv(double array[][10], 
	  long int NORDER, 
	  double *det)
{
	/*
	 * SUBROUTINE matinv.F 
	 *
	 * SOURCE								       
	 *   BEVINGTON, PAGES 302-303.
	 *
	 * PURPOSE
	 *   INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT
	 *
	 * USAGE 
	 *   CALL matinv (array, NORDER, det)
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   array  - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE
	 *   NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT)
	 *   det    - DETERMINANT OF INPUT MATRIX
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   NONE
	 *
	 * COMMENT
	 *   DIMENSION STATEMENT VALID FOR NORDER UP TO 10
	 *
	 */
	long int I, 
	  IK[10], 
	  I_, 
	  j, 
	  JK[10], 
	  J_, 
	  k, 
	  K_, 
	  L;
	double AMAX, 
	  SAVE;

#	ifdef DEBUG_FUN
	fputs( "<+>matinv()\n", debug_fp );
#	endif

	*det = 1.;

	/* gary add following to stop use before defined */
	j = LONG_MAX;

	for( k=1; k <= NORDER; k++ )
	{
		K_ = k - 1;
		/* set following to impossible values to make sure set eventually within if */
		IK[K_] =  LONG_MAX;
		JK[K_] =  LONG_MAX;

		/* FIND LARGEST ELEMENT array(I,j) IN REST OF MATRIX
		 * */
		AMAX = 0.;
		while( TRUE )
		{
			for( I=k; I <= NORDER; I++ )
			{
				I_ = I - 1;
				for( j=k; j <= NORDER; j++ )
				{
					J_ = j - 1;
					if( fabs(AMAX) <= fabs(array[J_][I_]) )
					{
						AMAX = array[J_][I_];
						IK[K_] = I;
						JK[K_] = j;
					}
				}
			}

			/* INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN array(k,k)
			 * */
			if( AMAX == 0 )
				goto L_32;
			I = IK[K_];
			if( I >= k )
			{
				if( I != k )
				{
					for( j=1; j <= NORDER; j++ )
					{
						J_ = j - 1;
						SAVE = array[J_][K_];
						array[J_][K_] = array[J_][I-1];
						array[J_][I-1] = -SAVE;
					}
				}
				j = JK[K_];
				if( j >= k )
					break;
			}
		}

		if( j != k )
		{
			for( I=1; I <= NORDER; I++ )
			{
				I_ = I - 1;
				SAVE = array[K_][I_];
				array[K_][I_] = array[j-1][I_];
				array[j-1][I_] = -SAVE;
			}
		}

		/* ACCUMULATE ELEMENTS OF INVERSE MATRIX 
		 * */
		for( I=1; I <= NORDER; I++ )
		{
			I_ = I - 1;
			if( I != k )
				array[K_][I_] = -array[K_][I_]/AMAX;
		}
		for( I=1; I <= NORDER; I++ )
		{
			I_ = I - 1;
			for( j=1; j <= NORDER; j++ )
			{
				J_ = j - 1;
				if( I != k )
				{
					if( j != k )
						array[J_][I_] += array[K_][I_]*array[J_][K_];
				}
			}
		}
		for( j=1; j <= NORDER; j++ )
		{
			J_ = j - 1;
			if( j != k )
				array[J_][K_] /= AMAX;
		}
		array[K_][K_] = 1./AMAX;
		*det *= AMAX;
	}

	/* RESTORE ORDERING OF MATRIX
	 * */
	for( L=1; L <= NORDER; L++ )
	{
		k = NORDER - L + 1;
		j = IK[k-1];
		if( j > k )
		{
			for( I=1; I <= NORDER; I++ )
			{
				I_ = I - 1;
				SAVE = array[k-1][I_];
				array[k-1][I_] = -array[j-1][I_];
				array[j-1][I_] = SAVE;
			}
		}
		I = JK[k-1];
		if( I > k )
		{
			for( j=1; j <= NORDER; j++ )
			{
				J_ = j - 1;
				SAVE = array[J_][k-1];
				array[J_][k-1] = -array[J_][I-1];
				array[J_][I-1] = SAVE;
			}
		}
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->matinv()\n", debug_fp );
#	endif
	return;
L_32:
	*det = 0.;
	
#	ifdef DEBUG_FUN
	fputs( " <->matinv()\n", debug_fp );
#	endif
	return;
}

double pbinom(long int nobs, 
	  long int ntotal, 
	  double prob)
{
	/*
	 * function pbinom.f
	 *
	 * source
	 *   Bevington, page 31.
	 *
	 * purpose
	 *   evaluate binomial probability coefficient
	 *
	 * usage
	 *   result = pbinom (nobs, ntotal, prob)
	 *
	 * description of parameters
	 *   nobs   - number of items observed
	 *   ntotal - total number of items
	 *   prob   - probability of observing each item
	 *
	 * subroutines and function subprograms required
	 *   factor (n)
	 *      calculates n factorial for integers
	 *
	 */
	long int notobs;
	double pbinom_v;

#	ifdef DEBUG_FUN
	fputs( "<+>pbinom()\n", debug_fp );
#	endif
	notobs = ntotal - nobs;
	pbinom_v = 
		factor(ntotal)/
		(factor(nobs)*factor(notobs))*
		(powi(prob,nobs))*
		powi(1. - prob,notobs);
	
#	ifdef DEBUG_FUN
	fputs( " <->pbinom()\n", debug_fp );
#	endif
	return( pbinom_v );
}

double pchisq(double chisqr, 
	  long int nfree)
{
	/*
	 * function pchisq.f
	 *
	 * source
	 *   Bevington, pages 192-193.
	 *
	 * purpose
	 *   evaluate probability for exceeding chi square
	 *
	 * usage
	 *   result = pchisq (chisqr, nfree)
	 *
	 * description of parameters
	 *   chisqr - comparison value of reduced chi square
	 *   nfree  - number of degrees of freedom
	 *
	 * subroutines and function subprograms required
	 *   bev_gamma (x)
	 *      calculates gamma function for integers and half-integers
	 *
	 * comments
	 *   calculation is approximate for nfree odd and
	 *      chi square greater than 50
	 *
	 */
	long int i, 
	  imax, 
	  neven;
	double fi, 
	  freedom, 
	  pchisq_v, 
	  pwr;
	double sum, 
	  term, 
	  z;

#	ifdef DEBUG_FUN
	fputs( "<+>pchisq()\n", debug_fp );
#	endif
	if( nfree > 0 )
	{
		freedom = (double)nfree;
		z = chisqr*freedom/2.;
		neven = 2*(nfree/2);
		if( nfree > neven )
		{

			/* number of degrees of freedom is odd
			 * */
			if( z <= 25 )
			{
				pwr = freedom/2.;
				term = 1.;
				sum = term/pwr;
				for( i=1; i <= 1000; i++ )
				{
					fi = (double)i;
					term = -term*z/fi;
					sum += term/(pwr + fi);
					if( fabs(term/sum) <= 0.00001 )
						break;
				}
				pchisq_v = 1. - (pow(z,pwr))*sum/bev_gamma(pwr);
				
#	ifdef DEBUG_FUN
				fputs( " <->pchisq()\n", debug_fp );
#	endif
				return( pchisq_v );
			}
			else
			{
				z = chisqr*(freedom - 1.)/2.;
			}
		}

		/* number of degrees of freedom is even
		 * */
		imax = nfree/2;
		term = 1.;
		sum = 0.;
		for( i=1; i <= imax; i++ )
		{
			fi = (double)i;
			sum += term;
			term = term*z/fi;
		}
		pchisq_v = sum*exp(z);
	}
	else
	{
		pchisq_v = 0.;
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->pchisq()\n", debug_fp );
#	endif
	return( pchisq_v );
}

double pcorre(double r, 
	  long int npts)
{
	/*
	 * function pcorre.f
	 *
	 * source
	 *   Bevington, pages 124-125.
	 *
	 * purpose
	 *   evaluate probability for no correlation between two variables
	 *
	 * usage
	 *   result = pcorre (r, npts)
	 *
	 * description of parameters
	 *   r      - linear correlation coefficient
	 *   npts   - number of data points
	 *
	 * subroutines and function subprograms required
	 *   bev_gamma (x)
	 *      calculates gamma function for integers and half-integers
	 *
	 */
	long int i, 
	  imax, 
	  neven, 
	  nfree;
	double freedom, 
	  pcorre_v;
	double denom, 
	  fi, 
	  fnum, 
	  r2, 
	  sum, 
	  term;

#	ifdef DEBUG_FUN
	fputs( "<+>pcorre()\n", debug_fp );
#	endif

	/* gary add to make sure always defined */
	pcorre_v = -DBL_MAX;

	/* evaluate number of degrees of freedom
	 * */
	nfree = npts - 2;
	if( nfree > 0 )
	{
		r2 = r*r;
		if( 1. > r2 )
		{
			neven = 2*(nfree/2);
			if( nfree > neven )
			{

				/* number of degrees of freedom is odd
				 * */
				imax = (nfree - 3)/2;
				term = fabs(r)*sqrt(1.-r2);
				sum = atan(r2/term);
				if( imax >= 0 )
				{
					if( imax > 0 )
					{
						sum += term;
						for( i=1; i <= imax; i++ )
						{
							fnum = 2*i;
							denom = 2*i + 1;
							term = term*(1. - r2)*fnum/denom;
							sum += term;
						}
					}
					else
					{
						sum += term;
					}
				}
				pcorre_v = 1. - 0.6366197724*sum;
			}
			else
			{

				/* number of degrees of freedom is even
				 * */
				imax = (nfree - 2)/2;
				freedom = (double)nfree;
				term = fabs(r);
				sum = term;
				if( imax >= 0 )
				{
					if( imax > 0 )
					{
						for( i=1; i <= imax; i++ )
						{
							fi = (double)i;
							fnum = (double)(imax - i + 1);
							denom = 2*i + 1;
							term = -term*r2*fnum/fi;
							sum += term/denom;
						}
						pcorre_v = 1.128379167*(bev_gamma((freedom+1.)/2.)/
						  bev_gamma(freedom/2.));
						pcorre_v = 1. - pcorre_v*sum;
					}
					else
					{
						pcorre_v = 1. - term;
					}
				}
			}
			
#	ifdef DEBUG_FUN
			fputs( " <->pcorre()\n", debug_fp );
#	endif
			return( pcorre_v );
		}
	}
	pcorre_v = 0.;
	
#	ifdef DEBUG_FUN
	fputs( " <->pcorre()\n", debug_fp );
#	endif
	return( pcorre_v );
}

double pgauss(double x, 
	  double averag, 
	  double sigma)
{
	/*
	 * function pgauss.f
	 *
	 * source
	 *   Bevington, page 45.
	 *
	 * purpose
	 *   evaluate gaussian probability function
	 *
	 * usage
	 *   result = pgauss (x, averag, sigma)
	 *
	 * description of parameters
	 *   x      - value for which probability is to be evaluated
	 *   averag - mean of distribution
	 *   sigma  - standard deviation of distribution
	 *
	 * subroutines and function subprograms required
	 *   none
	 *
	 */
	double pgauss_v;
	double z;

#	ifdef DEBUG_FUN
	fputs( "<+>pgauss()\n", debug_fp );
#	endif
	z = (x - averag)/sigma;
	pgauss_v = 0.3989422804/sigma*exp(-z*z/2.);
	
#	ifdef DEBUG_FUN
	fputs( " <->pgauss()\n", debug_fp );
#	endif
	return( pgauss_v );
}

double ploren(double x, 
	  double averag, 
	  double width)
{
	/*
	 * function ploren.f
	 *
	 * source
	 *   Bevington, page 51.
	 *
	 * purpose
	 *   evaluate lorentzian probability function
	 *
	 * usage
	 *   result = ploren (x, averag, width)
	 *
	 * description of parameters
	 *   x      - value for which probability is to be evaluated
	 *   averag - mean of distribution
	 *   width  - full width at half maximum of distribution
	 *
	 * subroutines and function subprograms required
	 *   none
	 *
	 */
	double ploren_v;

#	ifdef DEBUG_FUN
	fputs( "<+>ploren()\n", debug_fp );
#	endif
	ploren_v = 0.1591549431*width/((x - averag)*(x - averag) + 
		width/2.*width/2.);
	
#	ifdef DEBUG_FUN
	fputs( " <->ploren()\n", debug_fp );
#	endif
	return( ploren_v );
}

void polfit(double X[], 
	  double Y[], 
	  double SIGMAY[], 
	  long int NPTS, 
	  long int nterms, 
	  long int MODE, 
	  double A[], 
	  double *CHISQR)
{
	/*
	 * SUBROUTINE POLFIT.F 
	 *
	 *
	 * SOURCE
	 *   BEVINGTON, PAGES 140-142.
	 *
	 *
	 * PURPOSE
	 *   MAKE A LEAST-SQUARES FIT TO DATA WITH A POLYNOMIAL CURVE
	 *      Y = A(1) + A(2)*X + A(3)*X**2 + A(3)*X**3 + . . .
	 *
	 * USAGE 
	 *   CALL POLFIT (X, Y, SIGMAY, NPTS, nterms, MODE, A, CHISQR)
	 *
	 * DESCRIPTION OF PARAMETERS
	 *   X	   - array OF DATA POINTS FOR INDEPENDENT VARIABLE
	 *   Y	   - array OF DATA POINTS FOR DEPENDENT VARIABLE
	 *   SIGMAY - array OF STANDARD DEVIATIONS FOR Y DATA POINTS
	 *   NPTS   - NUMBER OF PAIRS OF DATA POINTS
	 *   nterms - NUMBER OF COEFFICIENTS (DEGREE OF POLYNOMIAL + 1)
	 *   MODE   - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT
	 *	     +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
	 *	      0 (NO WEIGHTING) WEIGHT(I) = 1.
	 *	     -1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
	 *   A	   - array OF COEFFICIENTS OF POLYNOMIAL
	 *   CHISQR - REDUCED CHI SQUARE FOR FIT 
	 *
	 * SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED 
	 *   determ (array, NORDER)
	 *      EVALUATES THE DETERMINANT OF A SYMETRICAL
	 *      TWO-DIMENSIONAL MATRIX OF ORDER NORDER
	 *
	 * COMMENTS
	 *   DIMENSION STATEMENT VALID FOR nterms UP TO 10
	 *
	 */
	long int I, 
	  I_, 
	  j, 
	  J_, 
	  k, 
	  K_, 
	  L, 
	  L_, 
	  N, 
	  NMAX, 
	  N_;
	double DELTA, 
	  FREE, 
	  WEIGHT, 
	  XI, 
	  YI;
	double array[10][10], 
	  CHISQ, 
	  SUMX[19], 
	  SUMY[10], 
	  XTERM, 
	  YTERM;

#	ifdef DEBUG_FUN
	fputs( "<+>POLFIT()\n", debug_fp );
#	endif

	/* ACCUMULATE WEIGHTED SUMS
	 * */
	NMAX = 2*nterms - 1;
	for( N=1; N <= NMAX; N++ )
	{
		N_ = N - 1;
		SUMX[N_] = 0.;
	}
	for( j=1; j <= nterms; j++ )
	{
		J_ = j - 1;
		SUMY[J_] = 0.;
	}
	CHISQ = 0.;
	for( I=1; I <= NPTS; I++ )
	{
		I_ = I - 1;
		XI = X[I_];
		YI = Y[I_];
		WEIGHT = 1.;
		if( MODE != 0 )
		{
			if( MODE > 0 )
			{
				WEIGHT = 1./(SIGMAY[I_]*SIGMAY[I_]);
			}
			else if( YI < 0 )
			{
				WEIGHT = 1./(-YI);
			}
			else if( YI != 0 )
			{
				WEIGHT = 1./YI;
			}
		}
		XTERM = WEIGHT;
		for( N=1; N <= NMAX; N++ )
		{
			N_ = N - 1;
			SUMX[N_] += XTERM;
			XTERM *= XI;
		}
		YTERM = WEIGHT*YI;
		for( N=1; N <= nterms; N++ )
		{
			N_ = N - 1;
			SUMY[N_] += YTERM;
			YTERM *= XI;
		}
		CHISQ += WEIGHT*(YI*YI);
	}

	/* CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS 
	 * */
	for( j=1; j <= nterms; j++ )
	{
		J_ = j - 1;
		for( k=1; k <= nterms; k++ )
		{
			K_ = k - 1;
			N = j + k - 1;
			array[K_][J_] = SUMX[N-1];
		}
	}
	DELTA = determ(array,nterms);
	if( DELTA == 0 )
	{
		*CHISQR = 0.;
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			A[J_] = 0.;
		}
	}
	else
	{
		for( L=1; L <= nterms; L++ )
		{
			L_ = L - 1;
			for( j=1; j <= nterms; j++ )
			{
				J_ = j - 1;
				for( k=1; k <= nterms; k++ )
				{
					K_ = k - 1;
					N = j + k - 1;
					array[K_][J_] = SUMX[N-1];
				}
				array[L_][J_] = SUMY[J_];
			}
			A[L_] = determ(array,nterms)/DELTA;
		}

		/* CALCULATE CHI SQUARE
		 * */
		for( j=1; j <= nterms; j++ )
		{
			J_ = j - 1;
			CHISQ += -2.*A[J_]*SUMY[J_];
			for( k=1; k <= nterms; k++ )
			{
				K_ = k - 1;
				N = j + k - 1;
				CHISQ += A[J_]*A[K_]*SUMX[N-1];
			}
		}
		FREE = (double)(NPTS - nterms);
		*CHISQR = CHISQ/FREE;
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->POLFIT()\n", debug_fp );
#	endif
	return;
}

double ppoiss(long int nobs, 
	  double averag)
{
	/*
	 * function ppoiss.f
	 *
	 * source
	 *   Bevington, page 39.
	 *
	 * purpose
	 *   evaluate poisson probability function
	 *
	 * usage
	 *   result = ppoiss (nobs, averag)
	 *
	 * description of parameters
	 *   nobs   - number of items observed
	 *   averag - mean of distribution
	 *
	 * subroutines and function subprograms required
	 *   factor (n)
	 *      calculates n factorial for integers
	 *
	 */
	double ppoiss_v;

#	ifdef DEBUG_FUN
	fputs( "<+>ppoiss()\n", debug_fp );
#	endif
	ppoiss_v = ((powi(averag,nobs))/factor(nobs))*exp(-averag);
	
#	ifdef DEBUG_FUN
	fputs( " <->ppoiss()\n", debug_fp );
#	endif
	return( ppoiss_v );
}

void regres(double x[], 
	  double y[], 
	  double sigmay[], 
	  long int npts, 
	  long int nterms, 
	  long int m[], 
	  long int mode, 
	  double yfit[], 
	  double *a0, 
	  double a[], 
	  double *sigma0, 
	  double sigmaa[], 
	  double r[], 
	  double *rmul, 
	  double *chisqr, 
	  double *ftest)
{
	/*
	 * subroutine regres.f
	 *
	 * source
	 *   Bevington, pages 172-175.
	 *
	 * purpose
	 *   make a mulitple linear regression fit to data with a specified
	 *      function which is linear in coefficients
	 *
	 * usage
	 *   call regres (x, y, sigmay, npts, nterms, m, mode, yfit,
	 *      a0, a, sigma0, sigmaa, r, rmul, chisqr, ftest)
	 *
	 * description of parameters
	 *   x	   - array of points for independent variable
	 *   y	   - array of points for dependent variable
	 *   sigmay - array of standard deviations for y data points
	 *   npts   - number of pairs of data points
	 *   nterms - number of coefficients
	 *   m	   - array of inclusion/rejection criteria for BevFctn
	 *   mode   - determines method of weighting least-squares fit
	 *	     +1 (instrumental) weight(i) = 1./sigmay(i)**2
	 *	      0 (no weighting) weight(i) = 1.
	 *	     -1 (statistical)  weight(i) = 1./y(i)
	 *   yfit   - array of calculated values of y
	 *   a0	   - constant term
	 *   a	   - array of coefficients
	 *   sigma0 - standard deviation of a0
	 *   sigmaa - array of standard deviations for coefficients
	 *   r	   - array of linear correlation coefficients
	 *   rmul   - multiple linear correlation coefficient
	 *   chisqr - reduced chi square for fit
	 *   ftest  - value of f for test of fit
	 *
	 * subroutines and function subprograms required
	 *   BevFctn (x, i, j, m)
	 *      evaluates the function for the jth term and the ith data point
	 *      using the array m to specify terms in the function
	 *   matinv (array, nterms, det)
	 *      inverts a symmetric two-dimensional matrix of degree nterms
	 *      and calculates its determinant
	 *
	 * comments
	 * (dim npts changed 100->1000 21-may-84 dct)
	 *   dimension statement valid for npts up to 100 and nterms up to 10
	 *   sigmaag changed to sigmaa in statement following statement 132
	 *
	 */
	long int i, 
	  i_, 
	  j, 
	  j_, 
	  k, 
	  k_;
	double det, 
	  fnpts, 
	  free1, 
	  freej, 
	  freen, 
	  varnce, 
	  weight[1000], 
	  wmean;
	double array[10][10], 
	  chisq, 
	  sigma, 
	  sigmax[10], 
	  sum, 
	  xmean[10], 
	  ymean;

#	ifdef DEBUG_FUN
	fputs( "<+>regres()\n", debug_fp );
#	endif

	/* initialize sums and arrays
	 * */
	sum = 0.;
	ymean = 0.;
	sigma = 0.;
	chisq = 0.;
	*rmul = 0.;
	for( i=1; i <= npts; i++ )
	{
		i_ = i - 1;
		yfit[i_] = 0.;
	}
	for( j=1; j <= nterms; j++ )
	{
		j_ = j - 1;
		xmean[j_] = 0.;
		sigmax[j_] = 0.;
		r[j_] = 0.;
		a[j_] = 0.;
		sigmaa[j_] = 0.;
		for( k=1; k <= nterms; k++ )
		{
			k_ = k - 1;
			array[k_][j_] = 0.;
		}
	}

	/* accumulate weighted sums
	 * */
	for( i=1; i <= npts; i++ )
	{
		i_ = i - 1;
		weight[i_] = 1.;
		if( mode != 0 )
		{
			if( mode > 0 )
			{
				weight[i_] = 1./(sigmay[i_]*sigmay[i_]);
			}
			else if( y[i_] < 0 )
			{
				weight[i_] = 1./(-y[i_]);
			}
			else if( y[i_] != 0 )
			{
				weight[i_] = 1./y[i_];
			}
		}
		sum += weight[i_];
		ymean += weight[i_]*y[i_];
		for( j=1; j <= nterms; j++ )
		{
			j_ = j - 1;
			xmean[j_] += weight[i_]*BevFctn(x,i,j,m);
		}
	}
	ymean /= sum;
	for( j=1; j <= nterms; j++ )
	{
		j_ = j - 1;
		xmean[j_] /= sum;
	}
	fnpts = (double)npts;
	wmean = sum/fnpts;
	for( i=1; i <= npts; i++ )
	{
		i_ = i - 1;
		weight[i_] /= wmean;
	}

	/* accumulate matrices r and array
	 * */
	for( i=1; i <= npts; i++ )
	{
		i_ = i - 1;
		sigma += weight[i_]*(y[i_] - ymean)*(y[i_] - ymean);
		for( j=1; j <= nterms; j++ )
		{
			j_ = j - 1;
			sigmax[j_] += weight[i_]*powi(BevFctn(x,i,j,m) - xmean[j_],2);
			r[j_] += weight[i_]*(BevFctn(x,i,j,m) - xmean[j_])*(y[i_] - 
			  ymean);
			for( k=1; k <= j; k++ )
			{
				k_ = k - 1;
				array[k_][j_] += weight[i_]*(BevFctn(x,i,j,m) - xmean[j_])*
				  (BevFctn(x,i,k,m) - xmean[k_]);
			}
		}
	}
	free1 = (double)(npts - 1);
	sigma = sqrt(sigma/free1);
	for( j=1; j <= nterms; j++ )
	{
		j_ = j - 1;
		sigmax[j_] = sqrt(sigmax[j_]/free1);
		r[j_] /= free1*sigmax[j_]*sigma;
		for( k=1; k <= j; k++ )
		{
			k_ = k - 1;
			array[k_][j_] /= free1*sigmax[j_]*sigmax[k_];
			array[j_][k_] = array[k_][j_];
		}
	}

	/* invert symmetric matrix
	 * */
	matinv(array,nterms,&det);
	if( det == 0 )
	{
		*a0 = 0.;
		*sigma0 = 0.;
		*rmul = 0.;
		*chisqr = 0.;
		*ftest = 0.;
	}
	else
	{

		/* calculate coefficients, fit, and chi square
		 * */
		*a0 = ymean;
		for( j=1; j <= nterms; j++ )
		{
			j_ = j - 1;
			for( k=1; k <= nterms; k++ )
			{
				k_ = k - 1;
				a[j_] += r[k_]*array[k_][j_];
			}
			a[j_] = a[j_]*sigma/sigmax[j_];
			*a0 += -a[j_]*xmean[j_];
			for( i=1; i <= npts; i++ )
			{
				i_ = i - 1;
				yfit[i_] += a[j_]*BevFctn(x,i,j,m);
			}
		}
		for( i=1; i <= npts; i++ )
		{
			i_ = i - 1;
			yfit[i_] += *a0;
			chisq += weight[i_]*(y[i_] - yfit[i_])*(y[i_] - yfit[i_]);
		}
		freen = (double)(npts - nterms - 1);
		*chisqr = chisq*wmean/freen;

		/* calculate uncertainties
		 * */
		if( mode == 0 )
		{
			varnce = *chisqr;
		}
		else
		{
			varnce = 1./wmean;
		}
		for( j=1; j <= nterms; j++ )
		{
			j_ = j - 1;
			sigmaa[j_] = array[j_][j_]*varnce/
				(free1*sigmax[j_]*sigmax[j_]);
			sigmaa[j_] = sqrt(sigmaa[j_]);
			*rmul += a[j_]*r[j_]*sigmax[j_]/sigma;
		}
		freej = (double)nterms;
		/* +noao: When rmul = 1, the following division (stmt 135) would blow up.
		 *        It has been changed so ftest is set to -99999. in this case. */
		if( 1. == *rmul )
		{
			*ftest = -99999.;
		}
		else
		{
			*ftest = (*rmul/freej)/((1. - *rmul)/freen);
		}
		/* -noao */
		*rmul = sqrt(*rmul);
		*sigma0 = varnce/fnpts;
		for( j=1; j <= nterms; j++ )
		{
			j_ = j - 1;
			for( k=1; k <= nterms; k++ )
			{
				k_ = k - 1;
				*sigma0 += varnce*xmean[j_]*xmean[k_]*array[k_][j_]/
				  (free1*sigmax[j_]*sigmax[k_]);
			}
		}
		*sigma0 = sqrt(*sigma0);
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->regres()\n", debug_fp );
#	endif
	return;
}

void smooth(double y[], 
	  long int npts)
{
	/*
	 * subroutine smooth.f
	 *
	 * source
	 *   Bevington, page 260.
	 *
	 * purpose
	 *   smooth a set of data points by averaging adjacent channels
	 *
	 * usage
	 *   call smooth (y, npts)
	 *
	 * description of parameters
	 *   y	   - array of data points
	 *   npts   - number of data points
	 *
	 * subroutines and function subprograms required
	 *   none
	 *
	 */
	long int i, 
	  i_, 
	  imax;
	double yi, 
	  ynew;

#	ifdef DEBUG_FUN
	fputs( "<+>smooth()\n", debug_fp );
#	endif
	imax = npts - 1;
	yi = y[0];
	for( i=1; i <= imax; i++ )
	{
		i_ = i - 1;
		ynew = (yi + 2.*y[i_] + y[i_+1])/4.;
		yi = y[i_];
		y[i_] = ynew;
	}
	y[npts-1] = (yi + 3.*y[npts-1])/4.;
	
#	ifdef DEBUG_FUN
	fputs( " <->smooth()\n", debug_fp );
#	endif
	return;
}

void xfit(double x[], 
	  double sigmax[], 
	  long int npts, 
	  long int mode, 
	  double *xmean, 
	  double *sigmam, 
	  double *sigma)
{
	/*
	 * subroutine xfit.f
	 *
	 * source
	 *   Bevington, page 76.
	 *
	 * purpose
	 *   calculate the mean and estimated errors for a set of data points
	 *
	 * usage
	 *   call xfit (x, sigmax, npts, mode, xmean, sigmam, sigma)
	 *
	 * description of parameters
	 *   x      - array of data points
	 *   sigmax - array of standard deviations for data points
	 *   npts   - number of data points
	 *   mode   - determines method of weighting
	 *            +1 (instrumental) weight(i) = 1./sigmax(i)**2
	 *             0 (no weighting) weight(i) = 1.
	 *            -1 (statistical)  weight(i) = 1.
	 *   xmean  - weighted mean
	 *   sigmam - standard deviation of mean
	 *   sigma  - standard deviation of data
	 *
	 * subroutines and function subprograms required
	 *   none
	 *
	 */
	long int i, 
	  i_;
	double freedom, 
	  sum, 
	  sumx, 
	  weight;

#	ifdef DEBUG_FUN
	fputs( "<+>xfit()\n", debug_fp );
#	endif

	/* accumulate weighted sums
	 * */
	sum = 0.;
	sumx = 0.;
	*sigma = 0.;
	*sigmam = 0.;
	for( i=1; i <= npts; i++ )
	{
		i_ = i - 1;
		if( mode > 0 )
		{
			weight = 1./(sigmax[i_]*sigmax[i_]);
		}
		else
		{
			weight = 1.;
		}
		sum += weight;
		sumx += weight*x[i_];
	}

	/* evaluate mean and standard deviations
	 * */
	*xmean = sumx/sum;
	for( i=1; i <= npts; i++ )
	{
		i_ = i - 1;
		*sigma += (x[i_] - *xmean)*(x[i_] - *xmean);
	}
	freedom = (double)(npts - 1);
	*sigma = sqrt(*sigma/freedom);
	if( mode == 0 )
	{
		*sigmam = *sigma/sqrt(sum);
	}
	else if( mode > 0 )
	{
		*sigmam = sqrt(1./sum);
	}
	else
	{
		*sigmam = sqrt(*xmean/sum);
	}
	
#	ifdef DEBUG_FUN
	fputs( " <->xfit()\n", debug_fp );
#	endif
	return;
}

double BevFunctn( 
	  double x[], 
	  long i, 
	  double a[] )
{
	return( x[0] * i * a[0] );
}

double BevFctn(
	  double x[], 
	  long i, 
	  long j, 
	  long m[] )
{
	return ( x[0] * i * j * m[0] );
}
