/*Translated by FOR_C, v4.0 (W), on 01/13/1999 at 18:30:42 */ /*FOR_C Options SET: cd ct=a do=r fx=gmufilmstuxz" io=p op=ai s=djsz1b - prototypes */ #include #include #include #include "cddefines.h" #include "powi.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 = 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 = i; sum += log(fi); } factor_v = 3628800.*exp(sum); } else { /* n less than 11 * */ for( i=2; i <= n; i++ ) { fi = 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 = 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 gamma(double x) { /* * * source * Bevington, page 126. * * purpose * calculate the gamma function for integers and half-integers * * usage * result = 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( "<+>gamma()\n", debug_fp ); # endif /* integerize argument * */ n = (long)(x - 0.25); xn = 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 = i; sum += log(fi-0.5); } gamma_v = prod*639383.8623*exp(sum); # ifdef DEBUG_FUN fputs( " <->gamma()\n", debug_fp ); # endif return( gamma_v ); } else { for( i=1; i <= n; i++ ) { fi = i; prod *= fi - 0.5; } } } gamma_v = prod; } else { /* argument is integer * */ gamma_v = factor(n); } # ifdef DEBUG_FUN fputs( " <->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 = 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 = 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, p[10][100]; # ifdef DEBUG_FUN fputs( "<+>legfit()\n", debug_fp ); # endif /* 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 = 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 ) 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 = 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.; } # ifdef DEBUG_FUN fputs( " <->legfit()\n", debug_fp ); # endif return; L_103: *chisqr = 0.; # ifdef DEBUG_FUN fputs( " <->legfit()\n", debug_fp ); # endif return; } void linfit(double x[], double y[], double sigmay[], long int npts, long int mode, double *a, double *sigmaa, double *b, double *sigmab, double *r) { /* * 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) * a - y intercept of fitted straight line * sigmaa - standard deviation of a * b - slope of fitted straight line * sigmab - standard deviation of b * r - linear correlation coefficient * * subroutines and function subprograms required * none * */ long int i, 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=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; } } 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; *a = (sumx2*sumy - sumx*sumxy)/delta; *b = (sumxy*sum - sumx*sumy)/delta; if( mode == 0 ) { c = npts - 2; varnce = (sumy2 + *a**a*sum + *b**b*sumx2 - 2.*(*a*sumy + *b*sumxy - *a**b*sumx))/c; } else { varnce = 1.; } *sigmaa = sqrt(varnce*sumx2/delta); *sigmab = sqrt(varnce*sum/delta); *r = (sum*sumxy - sumx*sumy)/sqrt(delta*(sum*sumy2-sumy*sumy)); # ifdef DEBUG_FUN fputs( " <->linfit()\n", debug_fp ); # endif return; } 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 * 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 = 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 = i; term = -term*z/fi; sum += term/(pwr + fi); if( fabs(term/sum) <= 0.00001 ) break; } pchisq_v = 1. - (pow(z,pwr))*sum/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 = 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 * 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 = nfree; term = fabs(r); sum = term; if( imax >= 0 ) { if( imax > 0 ) { for( i=1; i <= imax; i++ ) { fi = i; fnum = imax - i + 1; denom = 2*i + 1; term = -term*r2*fnum/fi; sum += term/denom; } pcorre_v = 1.128379167*(gamma((freedom+1.)/2.)/ 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 = 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 = 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 = 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 = 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 = 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 = 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] ); }