diff options
| author | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
|---|---|---|
| committer | Joseph Hunkeler <jhunkeler@gmail.com> | 2015-07-08 20:46:52 -0400 |
| commit | fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4 (patch) | |
| tree | bdda434976bc09c864f2e4fa6f16ba1952b1e555 /math/bevington | |
| download | iraf-linux-fa080de7afc95aa1c19a6e6fc0e0708ced2eadc4.tar.gz | |
Initial commit
Diffstat (limited to 'math/bevington')
54 files changed, 2775 insertions, 0 deletions
diff --git a/math/bevington/README b/math/bevington/README new file mode 100644 index 00000000..d0b3f778 --- /dev/null +++ b/math/bevington/README @@ -0,0 +1,13 @@ +This directory cotnains the original Fortran source for the Bevington +routines, as copied from his book. + +Jun85 Removed the nonstandard & in col 1 UNIX style continuation and + replaced it with standard Fortran (col 6) continuation. + +Oct85 Added an extran continuation statement to chifit.f to avoid ambiguous + transfer of control to statement 70 + +Nov85 It was possible to get a floating point divide by zero in regres.f + when rmul=1 in the calculation in statement 135 of ftest. The value + of rmul (the correlation) is now tested before the division and a + value of ftest = -99999. is returned when rmul=1. diff --git a/math/bevington/agauss.f b/math/bevington/agauss.f new file mode 100644 index 00000000..6aefc503 --- /dev/null +++ b/math/bevington/agauss.f @@ -0,0 +1,40 @@ +c function agauss.f +c +c source +c Bevington, page 48. +c +c purpose +c evaluate integral of gaussian probability function +c +c usage +c result = agauss (x, averag, sigma) +c +c description of parameters +c x - limit for integral +c averag - mean of distribution +c sigma - standard deviation of distribution +c integration range is averag +/- z*sigma +c where z = abs(x-averag)/sigma +c +c subroutines and function subprograms required +c none +c + function agauss (x,averag,sigma) + double precision z,y2,term,sum,denom +11 z=abs(x-averag)/sigma + agauss=0. + if (z) 42,42,21 +21 term=0.7071067812*z +22 sum=term + y2=(z**2)/2. + denom=1. +c +c accumulate sums of terms +c +31 denom=denom+2. +32 term=term*(y2*2./denom) +33 sum=sum+term + if (term/sum-1.e-10) 41,41,31 +41 agauss=1.128379167*sum*dexp(-y2) +42 return + end diff --git a/math/bevington/area.f b/math/bevington/area.f new file mode 100644 index 00000000..ae89d35a --- /dev/null +++ b/math/bevington/area.f @@ -0,0 +1,79 @@ +c function area.f +c +c source +c Bevington, pages 272-273. +c +c purpose +c integrate the area beneath a set of data points +c +c usage +c result = area (x, y, npts, nterms) +c +c description of parameters +c x - array of data points for independent variable +c y - array of data points for dependent variable +c npts - number of pairs of data points +c nterms - number of terms in fitting polynomial +c +c subroutines and function subprograms required +c integ (x, y, nterms, i1, x1, x2 sum) +c fits a polynomial with nterms starting at i1 +c and integrates area from x1 to x2 +c + function area (x,y,npts,nterms) + double precision sum + dimension x(1),y(1) +11 sum=0. + if (npts-nterms) 21,21,13 +13 neven=2*(nterms/2) + idelta=nterms/2-1 + if (nterms-neven) 31,31,51 +c +c fit all points with one curve +c +21 x1=x(1) + x2=x(npts) +23 call integ (x,y,npts,1,x1,x2,sum) + goto 71 +c +c even number of terms +c +31 x1=x(1) + j=nterms-idelta + x2=x(j) + call integ (x,y,nterms,1,x1,x2,sum) + i1=npts-nterms+1 + j=i1+idelta + x1=x(j) + x2=x(npts) +39 call integ (x,y,nterms,i1,x1,x2,sum) + if (i1-2) 71,71,41 +41 imax=i1-1 + do 46 i=2,imax + j=i+idelta + x1=x(j) + x2=x(j+1) +46 call integ (x,y,nterms,i,x1,x2,sum) + goto 71 +c +c odd number of terms +c +51 x1=x(1) + j=nterms-idelta + x2=(x(j)+x(j-1))/2. + call integ (x,y,nterms,1,x1,x2,sum) + i1=npts-nterms+1 + j=i1+idelta + x1=(x(j)+x(j+1))/2. + x2=x(npts) +59 call integ (x,y,nterms,i1,x1,x2,sum) + if (i1-2) 71,71,61 +61 imax=i1-1 + do 66 i=2,imax + j=i+idelta + x1=(x(j+1)+x(j))/2. + x2=(x(j+2)+x(j+1))/2. +66 call integ (x,y,nterms,i,x1,x2,sum) +71 area=sum + return + end diff --git a/math/bevington/chifit.f b/math/bevington/chifit.f new file mode 100644 index 00000000..ff4b0669 --- /dev/null +++ b/math/bevington/chifit.f @@ -0,0 +1,169 @@ +C SUBROUTINE CHIFIT.F +C +C SOURCE +C BEVINGTON, PAGES 228-231. +C +C PURPOSE +C MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION +C WITH A PARABOLIC EXPANSION OF CHI SQUARE +C +C USAGE +C CALL CHIFIT (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA, +C SIGMAA, YFIT, CHISQR) +C +C DESCRIPTION OF PARAMETERS +C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE +C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE +C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS +C NPTS - NUMBER OF PAIRS OF DATA POINTS +C NTERMS - NUMBER OF PARAMETERS +C MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT +C +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2 +C 0 (NO WEIGHTING) WEIGHT(I) = 1. +C -1 (STATISTICAL) WEIGHT(I) = 1./Y(I) +C A - ARRAY OF PARAMETERS +C DELTAA - ARRAY OF INCREMENTS FOR PARAMETERS A +C SIGMAA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A +C YFIT - ARRAY OF CALCULATED VALUES OF Y +C CHISQR - REDUCED CHI SQUARE FOR FIT +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C FUNCTN (X, I, A) +C EVALUATES THE FITTING FUNCTION FOR THE ITH TERM +C FCHISQ (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) +C EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA +C MATINV (ARRAY, NTERMS, DET) +C INVERTS A SYMMETRIC TWO-DIMENSIONAL MATRIX OF DEGREE NTERMS +C AND CALCULATES ITS DETERMINANT +C +C COMMENTS +C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10 +C + SUBROUTINE CHIFIT (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA, + *SIGMAA,YFIT,CHISQR) + DOUBLE PRECISION ALPHA + DIMENSION X(1),Y(1),SIGMAY(1),A(1),DELTAA(1),SIGMAA(1), + *YFIT(1) + DIMENSION ALPHA(10,10),BETA(10),DA(10) + REAL FUNCTN + EXTERNAL FUNCTN +C +11 NFREE=NPTS-NTERMS + FREE=NFREE + IF (NFREE) 14,14,16 +14 CHISQR=0. + GOTO 120 +16 DO 17 I=1,NPTS +17 YFIT(I)=FUNCTN (X,I,A) + CHISQ1=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +C +C EVALUATE ALPHA AND BETA MATRICES +C +20 DO 60 J=1,NTERMS +C +C A(J) + DELTAA(J) +C +21 AJ=A(J) + A(J)=AJ+DELTAA(J) + DO 24 I=1,NPTS +24 YFIT(I)=FUNCTN (X,I,A) + CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) + ALPHA(J,J)=CHISQ2-2.*CHISQ1 + BETA(J)=-CHISQ2 +31 DO 50 K=1,NTERMS + IF (K-J) 33,50,36 +33 ALPHA(K,J)=(ALPHA(K,J)-CHISQ2)/2. + ALPHA(J,K)=ALPHA(K,J) + GOTO 50 +36 ALPHA(J,K)=CHISQ1-CHISQ2 +C +C A(J) + DELTAA(J) AND A(K) + DELTAA(K) +C +41 AK=A(K) + A(K)=AK+DELTAA(K) + DO 44 I=1,NPTS +44 YFIT(I)=FUNCTN (X,I,A) + CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) + ALPHA(J,K)=ALPHA(J,K)+CHISQ3 + A(K)=AK +50 CONTINUE +C +C A(J) - DELTAA(J) +C +51 A(J)=AJ-DELTAA(J) + DO 53 I=1,NPTS +53 YFIT(I)=FUNCTN (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. +60 CONTINUE +C +C ELIMINATE NEGATIVE CURVATURE +C +61 DO 70 J=1,NTERMS + IF (ALPHA(J,J)) 63,65,70 +63 ALPHA(J,J)=-ALPHA(J,J) + GOTO 66 +65 ALPHA(J,J)=0.01 +C Changed from DO 70 +66 DO 700 K=1,NTERMS +C Changed from 68, 70, 68 + IF (K-J) 68,700,68 +68 ALPHA(J,K)=0. + ALPHA(K,J)=0. +C New continuation statement +700 CONTINUE +70 CONTINUE +C +C INVERT MATRIX AND EVALUATE PARAMETER INCREMENTS +C +71 CALL MATINV (ALPHA,NTERMS,DET) + DO 76 J=1,NTERMS + DA(J)=0. +74 DO 75 K=1,NTERMS +75 DA(J)=DA(J)+BETA(K)*ALPHA(J,K) +76 DA(J)=0.2*DA(J)*DELTAA(J) +C +C MAKE SURE CHI SQUARE DECREASES +C +81 DO 82 J=1,NTERMS +82 A(J)=A(J)+DA(J) +83 DO 84 I=1,NPTS +84 YFIT(I)=FUNCTN (X,I,A) + CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) + IF (CHISQ1-CHISQ2) 87,91,91 +87 DO 89 J=1,NTERMS + DA(J)=DA(J)/2. +89 A(J)=A(J)-DA(J) + GOTO 83 +C +C INCREMENT PARAMETERS UNTIL CHI SQUARE STARTS TO INCREASE +C +91 DO 92 J=1,NTERMS +92 A(J)=A(J)+DA(J) + DO 94 I=1,NPTS +94 YFIT(I)=FUNCTN (X,I,A) + CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) + IF (CHISQ3-CHISQ2) 97,101,101 +97 CHISQ1=CHISQ2 + CHISQ2=CHISQ3 +99 GOTO 91 +C +C FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS +C +101 DELTA=1./(1.+(CHISQ1-CHISQ2)/(CHISQ3-CHISQ2))+0.5 + DO 104 J=1,NTERMS + A(J)=A(J)-DELTA*DA(J) +104 SIGMAA(J)=DELTAA(J)*SQRT(FREE*ALPHA(J,J)) + DO 106 I=1,NPTS +106 YFIT(I)=FUNCTN (X,I,A) + CHISQR=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +111 IF (CHISQ2-CHISQR) 112,120,120 +112 DO 113 J=1,NTERMS +113 A(J)=A(J)+(DELTA-1.)*DA(J) + DO 115 I=1,NPTS +115 YFIT(I)=FUNCTN (X,I,A) + CHISQR=CHISQ2 +120 RETURN + END diff --git a/math/bevington/curfit.f b/math/bevington/curfit.f new file mode 100644 index 00000000..523a4c17 --- /dev/null +++ b/math/bevington/curfit.f @@ -0,0 +1,128 @@ +C SUBROUTINE CURFIT.F +C +C SOURCE +C BEVINGTON, PAGES 237-239. +C +C PURPOSE +C MAKE A LEAST-SQUARES FIT TO A NON-LINEAR FUNCTION +C WITH A LINEARIZATION OF THE FITTING FUNCTION +C +C USAGE +C CALL CURFIT (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA, +C SIGMAA, FLAMDA, YFIT, CHISQR) +C +C DESCRIPTION OF PARAMETERS +C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE +C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE +C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS +C NPTS - NUMBER OF PAIRS OF DATA POINTS +C NTERMS - NUMBER OF PARAMETERS +C MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT +C +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2 +C 0 (NO WEIGHTING) WEIGHT(I) = 1. +C -1 (STATISTICAL) WEIGHT(I) = 1./Y(I) +C A - ARRAY OF PARAMETERS +C DELTAA - ARRAY OF INCREMENTS FOR PARAMETERS A +C SIGMAA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A +C FLAMDA - PROPORTION OF GRADIENT SEARCH INCLUDED +C YFIT - ARRAY OF CALCULATED VALUES OF Y +C CHISQR - REDUCED CHI SQUARE FOR FIT +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C FUNCTN (X, I, A) +C EVALUATES THE FITTING FUNCTION FOR THE ITH TERM +C FCHISQ (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) +C EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA +C FDERIV (X, I, A, DELTAA, NTERMS, DERIV) +C EVALUATES THE DERIVATIVES OF THE FITTING FUNCTION +C FOR THE ITH TERM WITH RESPECT TO EACH PARAMETER +C MATINV (ARRAY, NTERMS, DET) +C INVERTS A SYMMETRIC TWO-DIMENSIONAL MATRIX OF DEGREE NTERMS +C AND CALCULATES ITS DETERMINANT +C +C COMMENTS +C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10 +C SET FLAMDA = 0.001 AT BEGINNING OF SEARCH +C + SUBROUTINE CURFIT (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA, + *SIGMAA,FLAMDA,YFIT,CHISQR) + DOUBLE PRECISION ARRAY + DIMENSION X(1),Y(1),SIGMAY(1),A(1),DELTAA(1),SIGMAA(1), + *YFIT(1) + DIMENSION WEIGHT(100),ALPHA(10,10),BETA(10),DERIV(10), + *ARRAY(10,10),B(10) + REAL FUNCTN + EXTERNAL FUNCTN + +C +11 NFREE=NPTS-NTERMS + IF (NFREE) 13,13,20 +13 CHISQR=0. + GOTO 110 +C +C EVALUATE WEIGHTS +C +20 DO 30 I=1,NPTS +21 IF (MODE) 22,27,29 +22 IF (Y(I)) 25,27,23 +23 WEIGHT(I)=1./Y(I) + GOTO 30 +25 WEIGHT(I)=1./(-Y(I)) + GOTO 30 +27 WEIGHT(I)=1. + GOTO 30 +29 WEIGHT(I)=1./SIGMAY(I)**2 +30 CONTINUE +C +C EVALUATE ALPHA AND BETA MATRICES +C +31 DO 34 J=1,NTERMS + BETA(J)=0. + DO 34 K=1,J +34 ALPHA(J,K)=0. +41 DO 50 I=1,NPTS + CALL FDERIV (X,I,A,DELTAA,NTERMS,DERIV) + DO 46 J=1,NTERMS + BETA(J)=BETA(J)+WEIGHT(I)*(Y(I)-FUNCTN(X,I,A))*DERIV(J) + DO 46 K=1,J +46 ALPHA(J,K)=ALPHA(J,K)+WEIGHT(I)*DERIV(J)*DERIV(K) +50 CONTINUE +51 DO 53 J=1,NTERMS + DO 53 K=1,J +53 ALPHA(K,J)=ALPHA(J,K) +C +C EVALUATE CHI SQUARE AT STARTING POINT +C +61 DO 62 I=1,NPTS +62 YFIT(I)=FUNCTN (X,I,A) +63 CHISQ1=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +C +C INVERT MODIFIED CURVATURE MATRIX TO FIND NEW PARAMETERS +C +71 DO 74 J=1,NTERMS + DO 73 K=1,NTERMS +73 ARRAY(J,K)=ALPHA(J,K)/SQRT(ALPHA(J,J)*ALPHA(K,K)) +74 ARRAY(J,J)=1.+FLAMDA +80 CALL MATINV (ARRAY,NTERMS,DET) +81 DO 84 J=1,NTERMS + B(J)=A(J) + DO 84 K=1,NTERMS +84 B(J)=B(J)+BETA(K)*ARRAY(J,K)/SQRT(ALPHA(J,J)*ALPHA(K,K)) +C +C IF CHI SQUARE INCREASED, INCREASE FLAMDA AND TRY AGAIN +C +91 DO 92 I=1,NPTS +92 YFIT(I)=FUNCTN (X,I,B) +93 CHISQR=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) + IF (CHISQ1-CHISQR) 95,101,101 +95 FLAMDA=10.*FLAMDA + GOTO 71 +C +C EVALUATE PARAMETERS AND UNCERTAINTIES +C +101 DO 103 J=1,NTERMS + A(J)=B(J) +103 SIGMAA(J)=SQRT(ARRAY(J,J)/ALPHA(J,J)) + FLAMDA=FLAMDA/10. +110 RETURN + END diff --git a/math/bevington/determ.f b/math/bevington/determ.f new file mode 100644 index 00000000..4bdc2778 --- /dev/null +++ b/math/bevington/determ.f @@ -0,0 +1,54 @@ +C FUNCTION DETERM.F +C +C SOURCE +C BEVINGTON, PAGE 294. +C +C PURPOSE +C CALCULATE THE DETERMINANT OF A SQUARE MATRIX +C +C USAGE +C DET = DETERM (ARRAY, NORDER) +C +C DESCRIPTION OF PARAMETERS +C ARRAY - MATRIX +C NORDER - ORDER OF DETERMINANT (DEGREE OF MATRIX) +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C NONE +C +C COMMENTS +C THIS SUBPROGRAM DESTROYS THE INPUT MATRIX ARRAY +C DIMENSION STATEMENT VALID FOR NORDER UP TO 10 +C + FUNCTION DETERM (ARRAY,NORDER) + DOUBLE PRECISION ARRAY,SAVE + DIMENSION ARRAY(10,10) +C +10 DETERM=1. +11 DO 50 K=1,NORDER +C +C INTERCHANGE COLUMNS IF DIAGONAL ELEMENT IS ZERO +C + IF (ARRAY(K,K)) 41,21,41 +21 DO 23 J=K,NORDER + IF (ARRAY(K,J)) 31,23,31 +23 CONTINUE + DETERM=0. + GOTO 60 +31 DO 34 I=K,NORDER + SAVE=ARRAY(I,J) + ARRAY(I,J)=ARRAY(I,K) +34 ARRAY(I,K)=SAVE + DETERM=-DETERM +C +C SUBTRACT ROW K FROM LOWER ROWS TO GET DIAGONAL MATRIX +C +41 DETERM=DETERM*ARRAY(K,K) + IF (K-NORDER) 43,50,50 +43 K1=K+1 + DO 46 I=K1,NORDER + DO 46 J=K1,NORDER +46 ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K) +50 CONTINUE +60 RETURN + END diff --git a/math/bevington/factor.f b/math/bevington/factor.f new file mode 100644 index 00000000..23b23307 --- /dev/null +++ b/math/bevington/factor.f @@ -0,0 +1,39 @@ +c function factor.f +c +c source +c Bevington, page 32. +c +c purpose +c calculates factorial function for integers +c +c usage +c result = factor (n) +c +c description of parameters +c n - integer argument +c +c subroutines and function subprograms required +c none +c + function factor (n) + double precision fi,sum +11 factor=1. + if (n-1) 40,40,13 +13 if (n-10) 21,21,31 +c +c n less than 11 +c +21 do 23 i=2,n + fi=i +23 factor=factor*fi + goto 40 +c +c n greater than 10 +c +31 sum=0. + do 34 i=11,n + fi=i +34 sum=sum+dlog(fi) +35 factor=3628800.*dexp(sum) +40 return + end diff --git a/math/bevington/fchisq.f b/math/bevington/fchisq.f new file mode 100644 index 00000000..1f9a8dc2 --- /dev/null +++ b/math/bevington/fchisq.f @@ -0,0 +1,54 @@ +C FUNCTION FCHISQ.F +C +C SOURCE +C BEVINGTON, PAGE 194. +C +C PURPOSE +C EVALUATE REDUCED CHI SQUARE FOR FIT OF DATA +C FCHISQ = SUM ((Y-YFIT)**2 / SIGMA**2) / NFREE +C +C USAGE +C RESULT = FCHISQ (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) +C +C DESCRIPTION OF PARAMETERS +C Y - ARRAY OF DATA POINTS +C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR DATA POINTS +C NPTS - NUMBER OF DATA POINTS +C NFREE - NUMBER OF DEGREES OF FREEDOM +C MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT +C +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2 +C 0 (NO WEIGHTING) WEIGHT(I) = 1. +C -1 (STATISTICAL) WEIGHT(I) = 1./Y(I) +C YFIT - ARRAY OF CALCULATED VALUES OF Y +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C NONE +C + FUNCTION FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) + DOUBLE PRECISION CHISQ,WEIGHT + DIMENSION Y(1),SIGMAY(1),YFIT(1) +11 CHISQ=0. +12 IF (NFREE) 13,13,20 +13 FCHISQ=0. + GOTO 40 +C +C ACCUMULATE CHI SQUARE +C +20 DO 30 I=1,NPTS +21 IF (MODE) 22,27,29 +22 IF (Y(I)) 25,27,23 +23 WEIGHT=1./Y(I) + GOTO 30 +25 WEIGHT=1./(-Y(I)) + GOTO 30 +27 WEIGHT=1. + GOTO 30 +29 WEIGHT=1./SIGMAY(I)**2 +30 CHISQ=CHISQ+WEIGHT*(Y(I)-YFIT(I))**2 +C +C DIVIDE BY NUMBER OF DEGREES OF FREEDOM +C +31 FREE=NFREE +32 FCHISQ=CHISQ/FREE +40 RETURN + END diff --git a/math/bevington/fderiv.f b/math/bevington/fderiv.f new file mode 100644 index 00000000..720953cb --- /dev/null +++ b/math/bevington/fderiv.f @@ -0,0 +1,39 @@ +c subroutine fderiv.f (non-analytical) +c +c source +c Bevington, page 242. +c +c purpose +c evaluate derivatives of function for least-squares search +c for arbitrary function given by functn +c +c usage +c call fderiv (x, i, a, deltaa, nterms, deriv) +c +c description of parameters +c x - array of data points for independent variable +c i - index of data points +c a - array of parameters +c deltaa - array of parameter increments +c nterms - number of parameters +c deriv - derivatives of function +c +c subroutines and function subprograms required +c functn (x, i, a) +c evaluates the fitting function for the ith term +c + subroutine fderiv (x,i,a,deltaa,nterms,deriv) + dimension x(1),a(1),deltaa(1),deriv(1) + real FUNCTN + external FUNCTN + +11 do 18 j=1,nterms + aj=a(j) + delta=deltaa(j) + a(j)=aj+delta + yfit=functn(x,i,a) + a(j)=aj-delta + deriv(j)=(yfit-functn(x,i,a))/(2.*delta) +18 a(j)=aj + return + end diff --git a/math/bevington/gamma.f b/math/bevington/gamma.f new file mode 100644 index 00000000..9c393649 --- /dev/null +++ b/math/bevington/gamma.f @@ -0,0 +1,49 @@ +c function gamma.f +c +c source +c Bevington, page 126. +c +c purpose +c calculate the gamma function for integers and half-integers +c +c usage +c result = gamma (x) +c +c description of parameters +c x - integer or half-integer +c +c subroutines or function subprograms required +c factor (n) +c calculates n factorial for integers +c + function gamma (x) + double precision prod,sum,fi +c +c integerize argument +c +11 n=x-0.25 + xn=n +13 if (x-xn-0.75) 31,31,21 +c +c argument is integer +c +21 gamma=factor(n) + goto 60 +c +c argument is half-integer +c +31 prod=1.77245385 + if (n) 44,44,33 +33 if (n-10) 41,41,51 +41 do 43 i=1,n + fi=i +43 prod=prod*(fi-0.5) +44 gamma=prod + goto 60 +51 sum=0. + do 54 i=11,n + fi=i +54 sum=sum+dlog(fi-0.5) +55 gamma=prod*639383.8623*dexp(sum) +60 return + end diff --git a/math/bevington/gradls.f b/math/bevington/gradls.f new file mode 100644 index 00000000..f6b88d68 --- /dev/null +++ b/math/bevington/gradls.f @@ -0,0 +1,113 @@ +C SUBROUTINE GRADLS.F +C +C SOURCE +C BEVINGTON, PAGES 220-222. +C +C PURPOSE +C MAKE A GRADIENT-SEARCH LEAST-SQUARES FIT TO DATA WITH A +C SPECIFIED FUNCTION WHICH IS NOT LINEAR IN COEFFICIENTS +C +C USAGE +C CALL GRADLS (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA, +C YFIT, CHISQR) +C +C DESCRIPTION OF PARAMETERS +C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE +C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE +C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS +C NPTS - NUMBER OF PAIRS OF DATA POINTS +C NTERMS - NUMBER OF PARAMETERS +C MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT +C +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2 +C 0 (NO WEIGHTING) WEIGHT(I) = 1. +C -1 (STATISTICAL) WEIGHT(I) = 1./Y(I) +C A - ARRAY OF PARAMETERS +C DELTAA - ARRAY OF INCREMENTS FOR PARAMETERS A +C YFIT - ARRAY OF CALCULATED VALUES OF Y +C CHISQR - REDUCED CHI SQUARE FOR FIT +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C FUNCTN (X, I, A) +C EVALUATES THE FITTING FUNCTION FOR THE ITH TERM +C FCHISQ (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) +C EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA +C +C COMMENTS +C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10 +C + SUBROUTINE GRADLS (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA, + *YFIT,CHISQR) + DIMENSION X(1),Y(1),SIGMAY(1),A(1),DELTAA(1),YFIT(1) + DIMENSION GRAD(10) + REAL FUNCTN + EXTERNAL FUNCTN + +C +C EVALUATE CHI SQUARE AT BEGINNING +C +11 NFREE=NPTS-NTERMS + IF (NFREE) 13,13,21 +13 CHISQR=0. + GOTO 110 +21 DO 22 I=1,NPTS +22 YFIT(I)=FUNCTN (X,I,A) + CHISQ1=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +C +C EVALUATE GRADIENT OF CHI SQUARE +C +31 SUM=0. +32 DO 39 J=1,NTERMS + DELTA=0.1*DELTAA(J) + A(J)=A(J)+DELTA + DO 36 I=1,NPTS +36 YFIT(I)=FUNCTN (X,I,A) + A(J)=A(J)-DELTA + GRAD(J)=CHISQ1-FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +39 SUM=SUM+GRAD(J)**2 +41 DO 42 J=1,NTERMS +42 GRAD(J)=DELTAA(J)*GRAD(J)/SQRT(SUM) +C +C EVALUATE CHI SQUARE AT NEW POINT +C +51 DO 52 J=1,NTERMS +52 A(J)=A(J)+GRAD(J) +53 DO 54 I=1,NPTS +54 YFIT(I)=FUNCTN (X,I,A) + CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +C +C MAKE SURE CHI SQUARE DECREASES +C +61 IF (CHISQ1-CHISQ2) 62,62,71 +62 DO 64 J=1,NTERMS + A(J)=A(J)-GRAD(J) +64 GRAD(J)=GRAD(J)/2. + GOTO 51 +C +C INCREMENT PARAMETERS UNTIL CHI SQUARE STARTS TO INCREASE +C +71 DO 72 J=1,NTERMS +72 A(J)=A(J)+GRAD(J) + DO 74 I=1,NPTS +74 YFIT(I)=FUNCTN (X,I,A) +75 CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +76 IF (CHISQ3-CHISQ2) 81,91,91 +81 CHISQ1=CHISQ2 +82 CHISQ2=CHISQ3 + GOTO 71 +C +C FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS +C +91 DELTA=1./(1.+(CHISQ1-CHISQ2)/(CHISQ3-CHISQ2))+0.5 + DO 93 J=1,NTERMS +93 A(J)=A(J)-DELTA*GRAD(J) + DO 95 I=1,NPTS +95 YFIT(I)=FUNCTN (X,I,A) + CHISQR=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +101 IF (CHISQ2-CHISQR) 102,110,110 +102 DO 103 J=1,NTERMS +103 A(J)=A(J)+(DELTA-1.)*GRAD(J) +104 DO 105 I=1,NPTS +105 YFIT(I)=FUNCTN (X,I,A) +106 CHISQR=CHISQ2 +110 RETURN + END diff --git a/math/bevington/gridls.f b/math/bevington/gridls.f new file mode 100644 index 00000000..e62d3219 --- /dev/null +++ b/math/bevington/gridls.f @@ -0,0 +1,102 @@ +C SUBROUTINE GRIDLS.F +C +C SOURCE +C BEVINGTON, PAGES 212-213. +C +C PURPOSE +C MAKE A GRID-SEARCH LEAST-SQUARES FIT TO DATA WITH A SPECIFIED +C FUNCTION WHICH IS NOT LINEAR IN COEFFICIENTS +C +C USAGE +C CALL GRIDLS (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, DELTAA, +C SIGMAA, YFIT, CHISQR) +C +C DESCRIPTION OF PARAMETERS +C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE +C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE +C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS +C NPTS - NUMBER OF PAIRS OF DATA POINTS +C NTERMS - NUMBER OF PARAMETERS +C MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT +C +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2 +C 0 (NO WEIGHTING) WEIGHT(I) = 1. +C -1 (STATISTICAL) WEIGHT(I) = 1./Y(I) +C A - ARRAY OF PARAMETERS +C DELTAA - ARRAY OF INCREMENTS FOR PARAMETERS A +C SIGMAA - ARRAY OF STANDARD DEVIATIONS FOR PARAMETERS A +C YFIT - ARRAY OF CALCULATED VALUES OF Y +C CHISQR - REDUCED CHI SQUARE FOR FIT +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C FUNCTN (X, I, A) +C EVALUATES THE FITTING FUNCTION FOR THE ITH TERM +C FCHISQ (Y, SIGMAY, NPTS, NFREE, MODE, YFIT) +C EVALUATES REDUCED CHI SQUARED FOR FIT TO DATA +C +C COMMENTS +C DELTAA VALUES ARE MODIFIED BY THE PROGRAM +C + SUBROUTINE GRIDLS (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA, + *SIGMAA,YFIT,CHISQR) + DIMENSION X(1),Y(1),SIGMAY(1),A(1),DELTAA(1),SIGMAA(1), + *YFIT(1) + REAL FUNCTN + EXTERNAL FUNCTN + +C +11 NFREE=NPTS-NTERMS + FREE=NFREE + CHISQR=0. + IF (NFREE) 100,100,20 +20 DO 90 J=1,NTERMS +C +C EVALUATE CHI SQUARE AT FIRST TWO SEARCH POINTS +C +21 DO 22 I=1,NPTS +22 YFIT(I)=FUNCTN (X,I,A) +23 CHISQ1=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) + FN=0. + DELTA=DELTAA(J) +41 A(J)=A(J)+DELTA + DO 43 I=1,NPTS +43 YFIT(I)=FUNCTN (X,I,A) +44 CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +45 IF (CHISQ1-CHISQ2) 51,41,61 +C +C REVERSE DIRECTION OF SEARCH IF CHI SQUARE IS INCREASING +C +51 DELTA=-DELTA + A(J)=A(J)+DELTA + DO 54 I=1,NPTS +54 YFIT(I)=FUNCTN (X,I,A) + SAVE=CHISQ1 + CHISQ1=CHISQ2 +57 CHISQ2=SAVE +C +C INCREMENT A(J) UNTIL CHI SQUARE INCREASES +C +61 FN=FN+1. + A(J)=A(J)+DELTA + DO 64 I=1,NPTS +64 YFIT(I)=FUNCTN (X,I,A) + CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +66 IF (CHISQ3-CHISQ2) 71,81,81 +71 CHISQ1=CHISQ2 + CHISQ2=CHISQ3 + GOTO 61 +C +C FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS +C +81 DELTA=DELTA*(1./(1.+(CHISQ1-CHISQ2)/(CHISQ3-CHISQ2))+0.5) +82 A(J)=A(J)-DELTA +83 SIGMAA(J)=DELTAA(J)*SQRT(2./(FREE*(CHISQ3-2.*CHISQ2+CHISQ1))) +84 DELTAA(J)=DELTAA(J)*FN/3. +90 CONTINUE +C +C EVALUATE FIT AND CHI SQUARE FOR FINAL PARAMETERS +C +91 DO 92 I=1,NPTS +92 YFIT(I)=FUNCTN (X,I,A) +93 CHISQR=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT) +100 RETURN + END diff --git a/math/bevington/integ.f b/math/bevington/integ.f new file mode 100644 index 00000000..153bafe1 --- /dev/null +++ b/math/bevington/integ.f @@ -0,0 +1,58 @@ +c subroutine integ.f +c +c source +c Bevington, page 274. +c +c purpose +c integrate the area beneath two data points +c +c usage +c call integ (x, y, nterms, i1, x1, x2, sum) +c +c description of parameters +c x - array of data points for independent variable +c y - array of data points for dependent variable +c nterms - number of terms in fitting polymonial +c i1 - first data point for fitting polynomial +c x1 - first value of x for integration +c x2 - final value of x for integration +c +c subroutines and function subprograms required +c none +c +c comments +c dimension statement valid for nterms up to 10 +c + subroutine integ (x,y,nterms,i1,x1,x2,sum) + double precision xjk,array,a,denom,deltax,sum + dimension x(1),y(1) + dimension array(10,10) +c +c construct square matrix and invert +c +11 do 17 j=1,nterms + i=j+i1-1 + deltax=x(i)-x(i1) + xjk=1. + do 17 k=1,nterms + array(j,k)=xjk +17 xjk=xjk*deltax +21 call matinv (array,nterms,det) + if (det) 31,23,31 +23 imid=i1+nterms/2 + sum=sum+y(imid)*(x2-x1) + goto 40 +c +c evaluate coefficients and integrate +c +31 dx1=x1-x(i1) + dx2=x2-x(i1) +33 do 39 j=1,nterms + i=j+i1-1 + a=0. + do 37 k=1,nterms +37 a=a+y(i)*array(j,k) + denom=j +39 sum=sum+(a/denom)*(dx2**j-dx1**j) +40 return + end diff --git a/math/bevington/interp.f b/math/bevington/interp.f new file mode 100644 index 00000000..75ca1b23 --- /dev/null +++ b/math/bevington/interp.f @@ -0,0 +1,85 @@ +C SUBROUTINE INTERP.F +C +C SOURCE +C BEVINGTON, PAGES 266-267. +C +C PURPOSE +C INTERPOLATE BETWEEN DATA POINTS TO EVALUATE A FUNCTION +C +C USAGE +C CALL INTERP (X, Y, NPTS, NTERMS, XIN, YOUT) +C +C DESCRIPTION OF PARAMETERS +C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE +C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE +C NPTS - NUMBER OF PAIRS OF DATA POINTS +C NTERMS - NUMBER OF TERMS IN FITTING POLYNOMIAL +C XIN - INPUT VALUE OF X +C YOUT - INTERPOLATED VALUE OF Y +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C NONE +C +C COMMENTS +C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10 +C VALUE OF NTERMS MAY BE MODIFIED BY THE PROGRAM +C + SUBROUTINE INTERP (X,Y,NPTS,NTERMS,XIN,YOUT) + DOUBLE PRECISION DELTAX,DELTA,A,PROD,SUM + DIMENSION X(1),Y(1) + DIMENSION DELTA(10),A(10) +C +C SEARCH FOR APPROPRIATE VALUE OF X(1) +C +11 DO 19 I=1,NPTS + IF (XIN-X(I)) 13,17,19 +13 I1=I-NTERMS/2 + IF (I1) 15,15,21 +15 I1=1 + GOTO 21 +17 YOUT=Y(I) +18 GOTO 61 +19 CONTINUE + I1=NPTS-NTERMS+1 +21 I2=I1+NTERMS-1 + IF (NPTS-I2) 23,31,31 +23 I2=NPTS + I1=I2-NTERMS+1 +25 IF (I1) 26,26,31 +26 I1=1 +27 NTERMS=I2-I1+1 +C +C EVALUATE DEVIATIONS DELTA +C +31 DENOM=X(I1+1)-X(I1) + DELTAX=(XIN-X(I1))/DENOM + DO 35 I=1,NTERMS + IX=I1+I-1 +35 DELTA(I)=(X(IX)-X(I1))/DENOM +C +C ACCUMULATE COEFFICIENTS A +C +40 A(1)=Y(I1) +41 DO 50 K=2,NTERMS + PROD=1. + SUM=0. + IMAX=K-1 + IXMAX=I1+IMAX + DO 49 I=1,IMAX + J=K-I + PROD=PROD*(DELTA(K)-DELTA(J)) +49 SUM=SUM-A(J)/PROD +50 A(K)=SUM+Y(IXMAX)/PROD +C +C ACCUMULATE SUM OF EXPANSION +C +51 SUM=A(1) + DO 57 J=2,NTERMS + PROD=1. + IMAX=J-1 + DO 56 I=1,IMAX +56 PROD=PROD*(DELTAX-DELTA(I)) +57 SUM=SUM+A(J)*PROD +60 YOUT=SUM +61 RETURN + END diff --git a/math/bevington/legfit.f b/math/bevington/legfit.f new file mode 100644 index 00000000..eb00b824 --- /dev/null +++ b/math/bevington/legfit.f @@ -0,0 +1,173 @@ +c subroutine legfit.f +c +c source +c Bevington, pages 155-157. +c +c purpose +c make a least-squares fit to data with a legendre polynomial +c y = a(1) + a(2)*x + a(3)*(3x**2-1)/2 + ... +c = a(1) * (1. + b(2)*x + b(3)*(3x**2-1)/2 + ... ) +c where x = cos(theta) +c +c usage +c call legfit (theta, y, sigmay, npts, norder, neven, mode, ftest, +c yfit, a, sigmaa, b, sigmab, chisqr) +c +c description of parameters +c theta - array of angles (in degrees) of the data points +c y - array of data points for dependent variable +c sigmay - array of standard deivations for y data points +c npts - number of pairs of data points +c norder - highest order of polynomial (number of terms - 1) +c neven - determines odd or even character of polynomial +c +1 fits only to even terms +c 0 fits to all terms +c -1 fits only to odd terms (plus constant term) +c mode - determines mode of weighting least-squares fit +c +1 (instrumental) weight(i) = 1./sigmay(i)**2 +c 0 (no weighting) weight(i) = 1. +c -1 (statistical) weight(i) = 1./y(i) +c ftest - array of values of f(l) for an f test +c yfit - array of calculated values of y +c a - array of coefficients of polynomial +c sigmaa - array of standard deviations for coefficients +c b - array of normalized relative coefficients +c sigmab - array of standard deviations for relative coefficients +c chisqr - reduced chi square for fit +c +c subroutines and function subprograms required +c matinv (array, nterms, det) +c inverts a symmetric two-dimensional matrix of degree nterms +c and calculates its determinant +c +c comments +c dimension statement valid for npts up to 100 and order up to 9 +c dcos changed to cos in statement 31 +c + subroutine legfit (theta,y,sigmay,npts,norder,neven,mode, + *ftest,yfit,a,sigmaa,b,sigmab,chisqr) + double precision cosine,p,beta,alpha,chisq + dimension theta(1),y(1),sigmay(1),ftest(1),yfit(1), + *a(1),sigmaa(1),b(1),sigmab(1) + dimension weight(100),p(100,10),beta(10),alpha(10,10) +c +c accumulate weights and legendre polynomials +c +11 nterms=1 + ncoeff=1 + jmax=norder+1 +20 do 40 i=1,npts +21 if (mode) 22,27,29 +22 if (y(i)) 25,27,23 +23 weight(i)=1./y(i) + goto 31 +25 weight(i)=1./(-y(i)) + goto 31 +27 weight(i)=1. + goto 31 +29 weight(i)=1./sigmay(i)**2 +31 cosine=cos(0.01745329252*theta(i)) + p(i,1)=1. + p(i,2)=cosine + do 36 l=2,norder + fl=l +36 p(i,l+1)=((2.*fl-1.)*cosine*p(i,l)-(fl-1.)*p(i,l-1))/fl +40 continue +c +c accumulate matrices alpha and beta +c +51 do 54 j=1,nterms + beta(j)=0. + do 54 k=1,nterms +54 alpha(j,k)=0. +61 do 66 i=1,npts + do 66 j=1,nterms + beta(j)=beta(j)+p(i,j)*y(i)*weight(i) + do 66 k=j,nterms + alpha(j,k)=alpha(j,k)+p(i,j)*p(i,k)*weight(i) +66 alpha(k,j)=alpha(j,k) +c +c delete fixed coefficients +c +70 if (neven) 71,91,81 +71 do 76 j=3,nterms,2 + beta(j)=0. + do 75 k=1,nterms + alpha(j,k)=0. +75 alpha(k,j)=0. +76 alpha(j,j)=1. + goto 91 +81 do 86 j=2,nterms,2 + beta(j)=0. + do 85 k=1,nterms + alpha(j,k)=0. +85 alpha(k,j)=0. +86 alpha(j,j)=1. +c +c invert curvature matrix alpha +c +91 do 95 j=1,jmax + a(j)=0. + sigmaa(j)=0. + b(j)=0. +95 sigmab(j)=0. + do 97 i=1,npts +97 yfit(i)=0. +101 call matinv (alpha,nterms,det) + if (det) 111,103,111 +103 chisqr=0. + goto 170 +c +c calculate coefficients, fit, and chi square +c +111 do 115 j=1,nterms + do 113 k=1,nterms +113 a(j)=a(j)+beta(k)*alpha(j,k) + do 115 i=1,npts +115 yfit(i)=yfit(i)+a(j)*p(i,j) +121 chisq=0. + do 123 i=1,npts +123 chisq=chisq+(y(i)-yfit(i))**2*weight(i) + free=npts-ncoeff + chisqr=chisq/free +c +c test for end of fit +c +131 if (nterms-jmax) 132,151,151 +132 if (ncoeff-2) 133,134,141 +133 if (neven) 137,137,135 +134 if (neven) 135,137,135 +135 nterms=nterms+2 + goto 138 +137 nterms=nterms+1 +138 ncoeff=ncoeff+1 + chisq1=chisq + goto 51 +141 fvalue=(chisq1-chisq)/chisqr + if (ftest(nterms)-fvalue) 134,143,143 +143 if (neven) 144,146,144 +144 nterms=nterms-2 + goto 147 +146 nterms=nterms-1 +147 ncoeff=ncoeff-1 + jmax=nterms +149 goto 51 +c +c calculate remainder of output +c +151 if (mode) 152,154,152 +152 varnce=1. + goto 155 +154 varnce=chisqr +155 do 156 j=1,nterms +156 sigmaa(j)=dsqrt(varnce*alpha(j,j)) +161 if (a(1)) 162,170,162 +162 do 166 j=2,nterms + if (a(j)) 164,166,164 +164 b(j)=a(j)/a(1) +165 sigmab(j)=b(j)*dsqrt((sigmaa(j)/a(j))**2+(sigmaa(1)/a(1))**2 + *-2.*varnce*alpha(j,1)/(a(j)*a(1))) +166 continue + b(1)=1. +170 return + end diff --git a/math/bevington/linfit.f b/math/bevington/linfit.f new file mode 100644 index 00000000..132c62b0 --- /dev/null +++ b/math/bevington/linfit.f @@ -0,0 +1,79 @@ +c subroutine linfit.f +c +c source +c Bevington, pages 104-105. +c +c purpose +c make a least-squares fit to a data set with a straight line +c +c usage +c call linfit (x, y, sigmay, npts, mode, a, sigmaa, b, sigmab, r) +c +c description of parameters +c x - array of data points for independent variable +c y - array of data points for dependent variable +c sigmay - array of standard deviations for y data points +c npts - number of pairs of data points +c mode - determines method of weighting least-squares fit +c +1 (instrumental) weight(i) = 1./sigmay(i)**2 +c 0 (no weighting) weight(i) = 1. +c -1 (statistical) weight(i) = 1./y(i) +c a - y intercept of fitted straight line +c sigmaa - standard deviation of a +c b - slope of fitted straight line +c sigmab - standard deviation of b +c r - linear correlation coefficient +c +c subroutines and function subprograms required +c none +c + subroutine linfit (x,y,sigmay,npts,mode,a,sigmaa,b,sigmab,r) + double precision sum,sumx,sumy,sumx2,sumxy,sumy2 + double precision xi,yi,weight,delta,varnce + dimension x(1),y(1),sigmay(1) +c +c accumulate weighted sums +c +11 sum=0. + sumx=0. + sumy=0. + sumx2=0. + sumxy=0. + sumy2=0. +21 do 50 i=1,npts + xi=x(i) + yi=y(i) + if (mode) 31,36,38 +31 if (yi) 34,36,32 +32 weight=1./yi + goto 41 +34 weight=1./(-yi) + goto 41 +36 weight=1. + goto 41 +38 weight=1./sigmay(i)**2 +41 sum=sum+weight + sumx=sumx+weight*xi + sumy=sumy+weight*yi + sumx2=sumx2+weight*xi*xi + sumxy=sumxy+weight*xi*yi + sumy2=sumy2+weight*yi*yi +50 continue +c +c calculate coefficients and standard deviations +c +51 delta=sum*sumx2-sumx*sumx + a=(sumx2*sumy-sumx*sumxy)/delta +53 b=(sumxy*sum-sumx*sumy)/delta +61 if (mode) 62,64,62 +62 varnce=1. + goto 67 +64 c=npts-2 + varnce=(sumy2+a*a*sum+b*b*sumx2 + *-2.*(a*sumy+b*sumxy-a*b*sumx))/c +67 sigmaa=dsqrt(varnce*sumx2/delta) +68 sigmab=dsqrt(varnce*sum/delta) +71 r=(sum*sumxy-sumx*sumy)/ + *dsqrt(delta*(sum*sumy2-sumy*sumy)) + return + end diff --git a/math/bevington/man/agauss.3m b/math/bevington/man/agauss.3m new file mode 100644 index 00000000..4f66e826 --- /dev/null +++ b/math/bevington/man/agauss.3m @@ -0,0 +1,24 @@ +.TH AGAUSS 3M +.SH NAME +agauss +.SH DESCRIPTION +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 diff --git a/math/bevington/man/area.3m b/math/bevington/man/area.3m new file mode 100644 index 00000000..1505ef4d --- /dev/null +++ b/math/bevington/man/area.3m @@ -0,0 +1,25 @@ +.TH AREA 3M +.SH NAME +area +.SH DESCRIPTION +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 diff --git a/math/bevington/man/chifit.3m b/math/bevington/man/chifit.3m new file mode 100644 index 00000000..ec8f7f2c --- /dev/null +++ b/math/bevington/man/chifit.3m @@ -0,0 +1,44 @@ +.TH CHIFIT 3M +.SH NAME +chifit +.SH DESCRIPTION +subroutine chifit.f + +source + Bevington, pages 228-231. + +purpose + make 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 + functn (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 symmetric two-dimension matrix of degree nterms + and calculates its determinant + +comments + valid for nterms up to 10 diff --git a/math/bevington/man/curfit.3m b/math/bevington/man/curfit.3m new file mode 100644 index 00000000..799f0545 --- /dev/null +++ b/math/bevington/man/curfit.3m @@ -0,0 +1,49 @@ +.TH CURFIT 3M +.SH NAME +curfit +.SH DESCRIPTION +subroutine curfit.f + +source + Bevington, pages 237-239. + +purpose + make 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 + functn (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 symmetric two-dimension matrix of degree nterms + and calculates its determinant + +comments + valid for nterms up to 10 + set flamda = 0.001 at beginning of search diff --git a/math/bevington/man/determ.3m b/math/bevington/man/determ.3m new file mode 100644 index 00000000..5ec86d06 --- /dev/null +++ b/math/bevington/man/determ.3m @@ -0,0 +1,25 @@ +.TH DETERM 3M +.SH NAME +determ +.SH DESCRIPTION +function determ.f + +source + Bevington, page 294. + +purpose + calculate the determinant of a square matrix + +usage + result = 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 + valid for norder up to 10 diff --git a/math/bevington/man/factor.3m b/math/bevington/man/factor.3m new file mode 100644 index 00000000..1127c5d0 --- /dev/null +++ b/math/bevington/man/factor.3m @@ -0,0 +1,20 @@ +.TH FACTOR 3M +.SH NAME +factor +.SH DESCRIPTION +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 diff --git a/math/bevington/man/fchisq.3m b/math/bevington/man/fchisq.3m new file mode 100644 index 00000000..0511409b --- /dev/null +++ b/math/bevington/man/fchisq.3m @@ -0,0 +1,29 @@ +.TH FCHISQ 3M +.SH NAME +fchisq +.SH DESCRIPTION +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 diff --git a/math/bevington/man/fderiv.3m b/math/bevington/man/fderiv.3m new file mode 100644 index 00000000..d9c88aeb --- /dev/null +++ b/math/bevington/man/fderiv.3m @@ -0,0 +1,27 @@ +.TH FDERIV 3M +.SH NAME +fderiv +.SH DESCRIPTION +subroutine fderiv.f (non-analytical) + +source + Bevington, page 242. + +purpose + evaluate derivatives of function for least-squares search + for arbitrary function given by functn + +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 + functn (x, i, a) + evaluates the fitting function for the ith term diff --git a/math/bevington/man/gamma.3m b/math/bevington/man/gamma.3m new file mode 100644 index 00000000..31939ef3 --- /dev/null +++ b/math/bevington/man/gamma.3m @@ -0,0 +1,21 @@ +.TH GAMMA 3M +.SH NAME +gamma +.SH DESCRIPTION +function gamma.f + +source + Bevington, page 126. + +purpose + calculate gamma function for integers and 1/2-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 diff --git a/math/bevington/man/gradls.3m b/math/bevington/man/gradls.3m new file mode 100644 index 00000000..681c93f4 --- /dev/null +++ b/math/bevington/man/gradls.3m @@ -0,0 +1,40 @@ +.TH GRADLS 3M +.SH NAME +gradls +.SH DESCRIPTION +subroutine gradls.f + +source + Bevington, pages 220-222. + +purpose + make 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 + functn (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 + valid for nterms up to 10 diff --git a/math/bevington/man/gridls.3m b/math/bevington/man/gridls.3m new file mode 100644 index 00000000..331a566c --- /dev/null +++ b/math/bevington/man/gridls.3m @@ -0,0 +1,41 @@ +.TH GRIDLS 3M +.SH NAME +gridls +.SH DESCRIPTION +subroutine gridls.f + +source + Bevington, pages 212-213. + +purpose + make grid-search least-squares fit to data with 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 + functn (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 diff --git a/math/bevington/man/integ.3m b/math/bevington/man/integ.3m new file mode 100644 index 00000000..e6200645 --- /dev/null +++ b/math/bevington/man/integ.3m @@ -0,0 +1,28 @@ +.TH INTEG 3M +.SH NAME +integ +.SH DESCRIPTION +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 + valid for nterms up to 10 diff --git a/math/bevington/man/interp.3m b/math/bevington/man/interp.3m new file mode 100644 index 00000000..0a25b2b0 --- /dev/null +++ b/math/bevington/man/interp.3m @@ -0,0 +1,29 @@ +.TH INTERP 3M +.SH NAME +interp +.SH DESCRIPTION +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 + valid for nterms up to 10 + value of nterms may be modified by the program diff --git a/math/bevington/man/legfit.3m b/math/bevington/man/legfit.3m new file mode 100644 index 00000000..efe61c07 --- /dev/null +++ b/math/bevington/man/legfit.3m @@ -0,0 +1,49 @@ +.TH LEGFIT 3M +.SH NAME +legfit +.SH DESCRIPTION +subroutine legfit.f + +source + Bevington, pages 155-157. + +purpose + make 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 std. deviations of relative coefficients + chisqr - reduced chi square for fit + +subroutines and function subprograms required + matinv (array, nterms, det) + inverts symmetric two-dimension matrix of degree nterms + and calculates its determinant + +comments + valid for npts up to 100 and order up to 9 + one source line modified - cos substituted for dcos diff --git a/math/bevington/man/linfit.3m b/math/bevington/man/linfit.3m new file mode 100644 index 00000000..e6e63f75 --- /dev/null +++ b/math/bevington/man/linfit.3m @@ -0,0 +1,33 @@ +.TH LINFIT 3M +.SH NAME +linfit +.SH DESCRIPTION +subroutine linfit.f + +source + Bevington, pages 104-105. + +purpose + make 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 diff --git a/math/bevington/man/matinv.3m b/math/bevington/man/matinv.3m new file mode 100644 index 00000000..6b2d1c44 --- /dev/null +++ b/math/bevington/man/matinv.3m @@ -0,0 +1,25 @@ +.TH MATINV 3M +.SH NAME +matinv +.SH DESCRIPTION +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 + +comments + valid for norder up to 10 diff --git a/math/bevington/man/pbinom.3m b/math/bevington/man/pbinom.3m new file mode 100644 index 00000000..344a62b3 --- /dev/null +++ b/math/bevington/man/pbinom.3m @@ -0,0 +1,23 @@ +.TH PBINOM 3M +.SH NAME +pbinom +.SH DESCRIPTION +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 diff --git a/math/bevington/man/pchisq.3m b/math/bevington/man/pchisq.3m new file mode 100644 index 00000000..1f4c6742 --- /dev/null +++ b/math/bevington/man/pchisq.3m @@ -0,0 +1,26 @@ +.TH PCHISQ 3M +.SH NAME +pchisq +.SH DESCRIPTION +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 1/2-integers + +comments + calculation is approximate for nfree odd and + chi square greater than 50 diff --git a/math/bevington/man/pcorre.3m b/math/bevington/man/pcorre.3m new file mode 100644 index 00000000..832c6640 --- /dev/null +++ b/math/bevington/man/pcorre.3m @@ -0,0 +1,22 @@ +.TH PCORRE 3M +.SH NAME +pcorre +.SH DESCRIPTION +function pcorre.f + +source + Bevington, pages 124-125. + +purpose + evaluate probability of no correlation between 2 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 1/2-integers diff --git a/math/bevington/man/pgauss.3m b/math/bevington/man/pgauss.3m new file mode 100644 index 00000000..4f0c7c7d --- /dev/null +++ b/math/bevington/man/pgauss.3m @@ -0,0 +1,22 @@ +.TH PGAUSS 3M +.SH NAME +pgauss +.SH DESCRIPTION +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 diff --git a/math/bevington/man/ploren.3m b/math/bevington/man/ploren.3m new file mode 100644 index 00000000..dfd485f5 --- /dev/null +++ b/math/bevington/man/ploren.3m @@ -0,0 +1,22 @@ +.TH PLOREN 3M +.SH NAME +ploren +.SH DESCRIPTION +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 diff --git a/math/bevington/man/polfit.3m b/math/bevington/man/polfit.3m new file mode 100644 index 00000000..9553e619 --- /dev/null +++ b/math/bevington/man/polfit.3m @@ -0,0 +1,36 @@ +.TH POLFIT 3M +.SH NAME +polfit +.SH DESCRIPTION +subroutine polfit.f + +source + Bevington, pages 140-142. + +purpose + make 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-dimension matrix of order norder + +comments + valid for nterms up to 10 diff --git a/math/bevington/man/ppoiss.3m b/math/bevington/man/ppoiss.3m new file mode 100644 index 00000000..0197768f --- /dev/null +++ b/math/bevington/man/ppoiss.3m @@ -0,0 +1,22 @@ +.TH PPOISS 3M +.SH NAME +ppoiss +.SH DESCRIPTION +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 diff --git a/math/bevington/man/regres.3m b/math/bevington/man/regres.3m new file mode 100644 index 00000000..e8cf076d --- /dev/null +++ b/math/bevington/man/regres.3m @@ -0,0 +1,49 @@ +.TH REGRES 3M +.SH NAME +regres +.SH DESCRIPTION +subroutine regres.f + +source + Bevington, pages 172-175. + +purpose + make mulitple linear regression fit to data with 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 fctn + 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 + fctn (x, i, j, m) + evaluates function for the jth term and the ith data + point using array m to specify terms in function + matinv (array, nterms, det) + inverts symmetric two-dimension matrix of degree nterms + and calculates its determinant + +comments + valid for npts up to 100 and nterms up to 10 + one source line modified - sigmaa substituted for sigmaag diff --git a/math/bevington/man/smooth.3m b/math/bevington/man/smooth.3m new file mode 100644 index 00000000..6db69d22 --- /dev/null +++ b/math/bevington/man/smooth.3m @@ -0,0 +1,21 @@ +.TH SMOOTH 3M +.SH NAME +smooth +.SH DESCRIPTION +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 diff --git a/math/bevington/man/xfit.3m b/math/bevington/man/xfit.3m new file mode 100644 index 00000000..53ed8ae8 --- /dev/null +++ b/math/bevington/man/xfit.3m @@ -0,0 +1,29 @@ +.TH XFIT 3M +.SH NAME +xfit +.SH DESCRIPTION +subroutine xfit.f + +source + Bevington, page 76. + +purpose + calculate mean and estimated errors for 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 diff --git a/math/bevington/matinv.f b/math/bevington/matinv.f new file mode 100644 index 00000000..5b1815c2 --- /dev/null +++ b/math/bevington/matinv.f @@ -0,0 +1,96 @@ +C SUBROUTINE MATINV.F +C +C SOURCE +C BEVINGTON, PAGES 302-303. +C +C PURPOSE +C INVERT A SYMMETRIC MATRIX AND CALCULATE ITS DETERMINANT +C +C USAGE +C CALL MATINV (ARRAY, NORDER, DET) +C +C DESCRIPTION OF PARAMETERS +C ARRAY - INPUT MATRIX WHICH IS REPLACED BY ITS INVERSE +C NORDER - DEGREE OF MATRIX (ORDER OF DETERMINANT) +C DET - DETERMINANT OF INPUT MATRIX +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C NONE +C +C COMMENT +C DIMENSION STATEMENT VALID FOR NORDER UP TO 10 +C + SUBROUTINE MATINV (ARRAY,NORDER,DET) + DOUBLE PRECISION ARRAY,AMAX,SAVE + DIMENSION ARRAY(10,10),IK(10),JK(10) +C +10 DET=1. +11 DO 100 K=1,NORDER +C +C FIND LARGEST ELEMENT ARRAY(I,J) IN REST OF MATRIX +C + AMAX=0. +21 DO 30 I=K,NORDER + DO 30 J=K,NORDER +23 IF (DABS(AMAX)-DABS(ARRAY(I,J))) 24,24,30 +24 AMAX=ARRAY(I,J) + IK(K)=I + JK(K)=J +30 CONTINUE +C +C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARRAY(K,K) +C +31 IF (AMAX) 41,32,41 +32 DET=0. + GOTO 140 +41 I=IK(K) + IF (I-K) 21,51,43 +43 DO 50 J=1,NORDER + SAVE=ARRAY(K,J) + ARRAY(K,J)=ARRAY(I,J) +50 ARRAY(I,J)=-SAVE +51 J=JK(K) + IF (J-K) 21,61,53 +53 DO 60 I=1,NORDER + SAVE=ARRAY(I,K) + ARRAY(I,K)=ARRAY(I,J) +60 ARRAY(I,J)=-SAVE +C +C ACCUMULATE ELEMENTS OF INVERSE MATRIX +C +61 DO 70 I=1,NORDER + IF (I-K) 63,70,63 +63 ARRAY(I,K)=-ARRAY(I,K)/AMAX +70 CONTINUE +71 DO 80 I=1,NORDER + DO 80 J=1,NORDER + IF (I-K) 74,80,74 +74 IF (J-K) 75,80,75 +75 ARRAY(I,J)=ARRAY(I,J)+ARRAY(I,K)*ARRAY(K,J) +80 CONTINUE +81 DO 90 J=1,NORDER + IF (J-K) 83,90,83 +83 ARRAY(K,J)=ARRAY(K,J)/AMAX +90 CONTINUE + ARRAY(K,K)=1./AMAX +100 DET=DET*AMAX +C +C RESTORE ORDERING OF MATRIX +C +101 DO 130 L=1,NORDER + K=NORDER-L+1 + J=IK(K) + IF (J-K) 111,111,105 +105 DO 110 I=1,NORDER + SAVE=ARRAY(I,K) + ARRAY(I,K)=-ARRAY(I,J) +110 ARRAY(I,J)=SAVE +111 I=JK(K) + IF (I-K) 130,130,113 +113 DO 120 J=1,NORDER + SAVE=ARRAY(K,J) + ARRAY(K,J)=-ARRAY(I,J) +120 ARRAY(I,J)=SAVE +130 CONTINUE +140 RETURN + END diff --git a/math/bevington/mkpkg b/math/bevington/mkpkg new file mode 100644 index 00000000..5d9bef3f --- /dev/null +++ b/math/bevington/mkpkg @@ -0,0 +1,35 @@ +# Make the Bevington library. + +$checkout libbev.a lib$ +$update libbev.a +$checkin libbev.a lib$ +$exit + +libbev.a: + agauss.f + area.f + chifit.f + curfit.f + determ.f + factor.f + fchisq.f + fderiv.f + gamma.f + gradls.f + gridls.f + integ.f + interp.f + legfit.f + linfit.f + matinv.f + pbinom.f + pchisq.f + pcorre.f + pgauss.f + ploren.f + polfit.f + ppoiss.f + regres.f + smooth.f + xfit.f + ; diff --git a/math/bevington/pbinom.f b/math/bevington/pbinom.f new file mode 100644 index 00000000..325f6f70 --- /dev/null +++ b/math/bevington/pbinom.f @@ -0,0 +1,26 @@ +c function pbinom.f +c +c source +c Bevington, page 31. +c +c purpose +c evaluate binomial probability coefficient +c +c usage +c result = pbinom (nobs, ntotal, prob) +c +c description of parameters +c nobs - number of items observed +c ntotal - total number of items +c prob - probability of observing each item +c +c subroutines and function subprograms required +c factor (n) +c calculates n factorial for integers +c + function pbinom (nobs,ntotal,prob) +1 notobs=ntotal-nobs +2 pbinom=factor(ntotal)/(factor(nobs)*factor(notobs))* + * (prob**nobs)*(1.-prob)**notobs + return + end diff --git a/math/bevington/pchisq.f b/math/bevington/pchisq.f new file mode 100644 index 00000000..0dbe6d9b --- /dev/null +++ b/math/bevington/pchisq.f @@ -0,0 +1,62 @@ +c function pchisq.f +c +c source +c Bevington, pages 192-193. +c +c purpose +c evaluate probability for exceeding chi square +c +c usage +c result = pchisq (chisqr, nfree) +c +c description of parameters +c chisqr - comparison value of reduced chi square +c nfree - number of degrees of freedom +c +c subroutines and function subprograms required +c gamma (x) +c calculates gamma function for integers and half-integers +c +c comments +c calculation is approximate for nfree odd and +c chi square greater than 50 +c + function pchisq (chisqr,nfree) + double precision z,term,sum +11 if (nfree) 12,12,14 +12 pchisq=0. + goto 60 +14 free=nfree + z=chisqr*free/2. + neven=2*(nfree/2) + if (nfree-neven) 21,21,41 +c +c number of degrees of freedom is even +c +21 imax=nfree/2 + term=1. + sum=0. +31 do 34 i=1,imax + fi=i + sum=sum+term +34 term=term*z/fi +35 pchisq=sum*dexp(z) + goto 60 +c +c number of degrees of freedom is odd +c +41 if (z-25) 44,44,42 +42 z=chisqr*(free-1.)/2. + goto 21 +44 pwr=free/2. + term=1. + sum=term/pwr +51 do 56 i=1,1000 + fi=i + term=-term*z/fi + sum=sum+term/(pwr+fi) +55 if (dabs(term/sum)-0.00001) 57,57,56 +56 continue +57 pchisq=1.-(z**pwr)*sum/gamma(pwr) +60 return + end diff --git a/math/bevington/pcorre.f b/math/bevington/pcorre.f new file mode 100644 index 00000000..9c0b0ae3 --- /dev/null +++ b/math/bevington/pcorre.f @@ -0,0 +1,69 @@ +c function pcorre.f +c +c source +c Bevington, pages 124-125. +c +c purpose +c evaluate probability for no correlation between two variables +c +c usage +c result = pcorre (r, npts) +c +c description of parameters +c r - linear correlation coefficient +c npts - number of data points +c +c subroutines and function subprograms required +c gamma (x) +c calculates gamma function for integers and half-integers +c + function pcorre (r,npts) + double precision r2,term,sum,fi,fnum,denom +c +c evaluate number of degrees of freedom +c +11 nfree=npts-2 + if (nfree) 13,13,15 +13 pcorre=0. + goto 60 +15 r2=r**2 + if (1.-r2) 13,13,17 +17 neven=2*(nfree/2) + if (nfree-neven) 21,21,41 +c +c number of degrees of freedom is even +c +21 imax=(nfree-2)/2 + free=nfree +23 term=abs(r) + sum=term + if (imax) 60,26,31 +26 pcorre=1.-term + goto 60 +31 do 36 i=1,imax + fi=i + fnum=imax-i+1 + denom=2*i+1 + term=-term*r2*fnum/fi +36 sum=sum+term/denom + pcorre=1.128379167*(gamma((free+1.)/2.)/gamma(free/2.)) + pcorre=1.-pcorre*sum + goto 60 +c +c number of degrees of freedom is odd +c +41 imax=(nfree-3)/2 +42 term=abs(r)*dsqrt(1.-r2) +43 sum=datan(r2/term) + if (imax) 57,45,51 +45 sum=sum+term + goto 57 +51 sum=sum+term +52 do 56 i=1,imax + fnum=2*i + denom=2*i+1 + term=term*(1.-r2)*fnum/denom +56 sum=sum+term +57 pcorre=1.-0.6366197724*sum +60 return + end diff --git a/math/bevington/pgauss.f b/math/bevington/pgauss.f new file mode 100644 index 00000000..56cf5d4e --- /dev/null +++ b/math/bevington/pgauss.f @@ -0,0 +1,25 @@ +c function pgauss.f +c +c source +c Bevington, page 45. +c +c purpose +c evaluate gaussian probability function +c +c usage +c result = pgauss (x, averag, sigma) +c +c description of parameters +c x - value for which probability is to be evaluated +c averag - mean of distribution +c sigma - standard deviation of distribution +c +c subroutines and function subprograms required +c none +c + function pgauss (x,averag,sigma) + double precision z +1 z=(x-averag)/sigma +2 pgauss=0.3989422804/sigma*dexp(-(z**2)/2.) + return + end diff --git a/math/bevington/ploren.f b/math/bevington/ploren.f new file mode 100644 index 00000000..20fcb7c3 --- /dev/null +++ b/math/bevington/ploren.f @@ -0,0 +1,23 @@ +c function ploren.f +c +c source +c Bevington, page 51. +c +c purpose +c evaluate lorentzian probability function +c +c usage +c result = ploren (x, averag, width) +c +c description of parameters +c x - value for which probability is to be evaluated +c averag - mean of distribution +c width - full width at half maximum of distribution +c +c subroutines and function subprograms required +c none +c + function ploren (x,averag,width) +1 ploren=0.1591549431*width/((x-averag)**2+(width/2.)**2) + return + end diff --git a/math/bevington/polfit.f b/math/bevington/polfit.f new file mode 100644 index 00000000..cfbdb178 --- /dev/null +++ b/math/bevington/polfit.f @@ -0,0 +1,100 @@ +C SUBROUTINE POLFIT.F +C +C SOURCE +C BEVINGTON, PAGES 140-142. +C +C PURPOSE +C MAKE A LEAST-SQUARES FIT TO DATA WITH A POLYNOMIAL CURVE +C Y = A(1) + A(2)*X + A(3)*X**2 + A(3)*X**3 + . . . +C +C USAGE +C CALL POLFIT (X, Y, SIGMAY, NPTS, NTERMS, MODE, A, CHISQR) +C +C DESCRIPTION OF PARAMETERS +C X - ARRAY OF DATA POINTS FOR INDEPENDENT VARIABLE +C Y - ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE +C SIGMAY - ARRAY OF STANDARD DEVIATIONS FOR Y DATA POINTS +C NPTS - NUMBER OF PAIRS OF DATA POINTS +C NTERMS - NUMBER OF COEFFICIENTS (DEGREE OF POLYNOMIAL + 1) +C MODE - DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT +C +1 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2 +C 0 (NO WEIGHTING) WEIGHT(I) = 1. +C -1 (STATISTICAL) WEIGHT(I) = 1./Y(I) +C A - ARRAY OF COEFFICIENTS OF POLYNOMIAL +C CHISQR - REDUCED CHI SQUARE FOR FIT +C +C SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED +C DETERM (ARRAY, NORDER) +C EVALUATES THE DETERMINANT OF A SYMETRICAL +C TWO-DIMENSIONAL MATRIX OF ORDER NORDER +C +C COMMENTS +C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10 +C + SUBROUTINE POLFIT (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR) + DOUBLE PRECISION SUMX,SUMY,XTERM,YTERM,ARRAY,CHISQ + DIMENSION X(1),Y(1),SIGMAY(1),A(1) + DIMENSION SUMX(19),SUMY(10),ARRAY(10,10) +C +C ACCUMULATE WEIGHTED SUMS +C +11 NMAX=2*NTERMS-1 + DO 13 N=1,NMAX +13 SUMX(N)=0. + DO 15 J=1,NTERMS +15 SUMY(J)=0. + CHISQ=0. +21 DO 50 I=1,NPTS + XI=X(I) + YI=Y(I) +31 IF (MODE) 32,37,39 +32 IF (YI) 35,37,33 +33 WEIGHT=1./YI + GOTO 41 +35 WEIGHT=1./(-YI) + GOTO 41 +37 WEIGHT=1. + GOTO 41 +39 WEIGHT=1./SIGMAY(I)**2 +41 XTERM=WEIGHT + DO 44 N=1,NMAX + SUMX(N)=SUMX(N)+XTERM +44 XTERM=XTERM*XI +45 YTERM=WEIGHT*YI + DO 48 N=1,NTERMS + SUMY(N)=SUMY(N)+YTERM +48 YTERM=YTERM*XI +49 CHISQ=CHISQ+WEIGHT*YI**2 +50 CONTINUE +C +C CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS +C +51 DO 54 J=1,NTERMS + DO 54 K=1,NTERMS + N=J+K-1 +54 ARRAY(J,K)=SUMX(N) + DELTA=DETERM (ARRAY,NTERMS) + IF (DELTA) 61,57,61 +57 CHISQR=0. + DO 59 J=1,NTERMS +59 A(J)=0. + GOTO 80 +61 DO 70 L=1,NTERMS +62 DO 66 J=1,NTERMS + DO 65 K=1,NTERMS + N=J+K-1 +65 ARRAY(J,K)=SUMX(N) +66 ARRAY(J,L)=SUMY(J) +70 A(L)=DETERM (ARRAY,NTERMS) /DELTA +C +C CALCULATE CHI SQUARE +C +71 DO 75 J=1,NTERMS + CHISQ=CHISQ-2.*A(J)*SUMY(J) + DO 75 K=1,NTERMS + N=J+K-1 +75 CHISQ=CHISQ+A(J)*A(K)*SUMX(N) +76 FREE=NPTS-NTERMS +77 CHISQR=CHISQ/FREE +80 RETURN + END diff --git a/math/bevington/ppoiss.f b/math/bevington/ppoiss.f new file mode 100644 index 00000000..4782fde8 --- /dev/null +++ b/math/bevington/ppoiss.f @@ -0,0 +1,23 @@ +c function ppoiss.f +c +c source +c Bevington, page 39. +c +c purpose +c evaluate poisson probability function +c +c usage +c result = ppoiss (nobs, averag) +c +c description of parameters +c nobs - number of items observed +c averag - mean of distribution +c +c subroutines and function subprograms required +c factor (n) +c calculates n factorial for integers +c + function ppoiss (nobs,averag) +1 ppoiss=((averag**nobs)/factor(nobs))*exp(-averag) + return + end diff --git a/math/bevington/regres.f b/math/bevington/regres.f new file mode 100644 index 00000000..ef23a574 --- /dev/null +++ b/math/bevington/regres.f @@ -0,0 +1,173 @@ +c subroutine regres.f +c +c source +c Bevington, pages 172-175. +c +c purpose +c make a mulitple linear regression fit to data with a specified +c function which is linear in coefficients +c +c usage +c call regres (x, y, sigmay, npts, nterms, m, mode, yfit, +c a0, a, sigma0, sigmaa, r, rmul, chisqr, ftest) +c +c description of parameters +c x - array of points for independent variable +c y - array of points for dependent variable +c sigmay - array of standard deviations for y data points +c npts - number of pairs of data points +c nterms - number of coefficients +c m - array of inclusion/rejection criteria for fctn +c mode - determines method of weighting least-squares fit +c +1 (instrumental) weight(i) = 1./sigmay(i)**2 +c 0 (no weighting) weight(i) = 1. +c -1 (statistical) weight(i) = 1./y(i) +c yfit - array of calculated values of y +c a0 - constant term +c a - array of coefficients +c sigma0 - standard deviation of a0 +c sigmaa - array of standard deviations for coefficients +c r - array of linear correlation coefficients +c rmul - multiple linear correlation coefficient +c chisqr - reduced chi square for fit +c ftest - value of f for test of fit +c +c subroutines and function subprograms required +c fctn (x, i, j, m) +c evaluates the function for the jth term and the ith data point +c using the array m to specify terms in the function +c matinv (array, nterms, det) +c inverts a symmetric two-dimensional matrix of degree nterms +c and calculates its determinant +c +c comments +c (dim npts changed 100->1000 21-may-84 dct) +c dimension statement valid for npts up to 100 and nterms up to 10 +c sigmaag changed to sigmaa in statement following statement 132 +c + subroutine regres (x,y,sigmay,npts,nterms,m,mode,yfit, + *a0,a,sigma0,sigmaa,r,rmul,chisqr,ftest) + double precision array,sum,ymean,sigma,chisq,xmean,sigmax + dimension x(1),y(1),sigmay(1),m(1),yfit(1),a(1),sigmaa(1), + *r(1) + dimension weight(1000),xmean(10),sigmax(10),array(10,10) + REAL FCTN + EXTERNAL FCTN + +c +c initialize sums and arrays +c +11 sum=0. + ymean=0. + sigma=0. + chisq=0. + rmul=0. + do 17 i=1,npts +17 yfit(i)=0. +21 do 28 j=1,nterms + xmean(j)=0. + sigmax(j)=0. + r(j)=0. + a(j)=0. + sigmaa(j)=0. + do 28 k=1,nterms +28 array(j,k)=0. +c +c accumulate weighted sums +c +30 do 50 i=1,npts +31 if (mode) 32,37,39 +32 if (y(i)) 35,37,33 +33 weight(i)=1./y(i) + goto 41 +35 weight(i)=1./(-y(i)) + goto 41 +37 weight(i)=1. + goto 41 +39 weight(i)=1./sigmay(i)**2 +41 sum=sum+weight(i) + ymean=ymean+weight(i)*y(i) + do 44 j=1,nterms +44 xmean(j)=xmean(j)+weight(i)*fctn(x,i,j,m) +50 continue +51 ymean=ymean/sum + do 53 j=1,nterms +53 xmean(j)=xmean(j)/sum + fnpts=npts + wmean=sum/fnpts + do 57 i=1,npts +57 weight(i)=weight(i)/wmean +c +c accumulate matrices r and array +c +61 do 67 i=1,npts + sigma=sigma+weight(i)*(y(i)-ymean)**2 + do 67 j=1,nterms + sigmax(j)=sigmax(j)+weight(i)*(fctn(x,i,j,m)-xmean(j))**2 + r(j)=r(j)+weight(i)*(fctn(x,i,j,m)-xmean(j))*(y(i)-ymean) + do 67 k=1,j +67 array(j,k)=array(j,k)+weight(i)*(fctn(x,i,j,m)-xmean(j))* + *(fctn(x,i,k,m)-xmean(k)) +71 free1=npts-1 +72 sigma=dsqrt(sigma/free1) + do 78 j=1,nterms +74 sigmax(j)=dsqrt(sigmax(j)/free1) + r(j)=r(j)/(free1*sigmax(j)*sigma) + do 78 k=1,j + array(j,k)=array(j,k)/(free1*sigmax(j)*sigmax(k)) +78 array(k,j)=array(j,k) +c +c invert symmetric matrix +c +81 call matinv (array,nterms,det) + if (det) 101,91,101 +91 a0=0. + sigma0=0. + rmul=0. + chisqr=0. + ftest=0. + goto 150 +c +c calculate coefficients, fit, and chi square +c +101 a0=ymean +102 do 108 j=1,nterms + do 104 k=1,nterms +104 a(j)=a(j)+r(k)*array(j,k) +105 a(j)=a(j)*sigma/sigmax(j) +106 a0=a0-a(j)*xmean(j) +107 do 108 i=1,npts +108 yfit(i)=yfit(i)+a(j)*fctn(x,i,j,m) +111 do 113 i=1,npts + yfit(i)=yfit(i)+a0 +113 chisq=chisq+weight(i)*(y(i)-yfit(i))**2 + freen=npts-nterms-1 +115 chisqr=chisq*wmean/freen +c +c calculate uncertainties +c +121 if (mode) 122,124,122 +122 varnce=1./wmean + goto 131 +124 varnce=chisqr +131 do 133 j=1,nterms +132 sigmaa(j)=array(j,j)*varnce/(free1*sigmax(j)**2) + sigmaa(j)=sqrt(sigmaa(j)) +133 rmul=rmul+a(j)*r(j)*sigmax(j)/sigma + freej=nterms +c +noao: When rmul = 1, the following division (stmt 135) would blow up. +c It has been changed so ftest is set to -99999. in this case. + if (1. - rmul) 135, 935, 135 +135 ftest=(rmul/freej)/((1.-rmul)/freen) + goto 136 +935 ftest = -99999. +c -noao +136 rmul=sqrt(rmul) +141 sigma0=varnce/fnpts + do 145 j=1,nterms + do 145 k=1,nterms +145 sigma0=sigma0+varnce*xmean(j)*xmean(k)*array(j,k)/ + *(free1*sigmax(j)*sigmax(k)) +146 sigma0=sqrt(sigma0) +150 return + end diff --git a/math/bevington/smooth.f b/math/bevington/smooth.f new file mode 100644 index 00000000..19c2c457 --- /dev/null +++ b/math/bevington/smooth.f @@ -0,0 +1,29 @@ +c subroutine smooth.f +c +c source +c Bevington, page 260. +c +c purpose +c smooth a set of data points by averaging adjacent channels +c +c usage +c call smooth (y, npts) +c +c description of parameters +c y - array of data points +c npts - number of data points +c +c subroutines and function subprograms required +c none +c + subroutine smooth (y,npts) + dimension y(1) +11 imax=npts-1 + yi=y(1) +21 do 24 i=1,imax + ynew=(yi+2.*y(i)+y(i+1))/4. + yi=y(i) +24 y(i)=ynew +25 y(npts)=(yi+3.*y(npts))/4. + return + end diff --git a/math/bevington/xfit.f b/math/bevington/xfit.f new file mode 100644 index 00000000..a097d4af --- /dev/null +++ b/math/bevington/xfit.f @@ -0,0 +1,59 @@ +c subroutine xfit.f +c +c source +c Bevington, page 76. +c +c purpose +c calculate the mean and estimated errors for a set of data points +c +c usage +c call xfit (x, sigmax, npts, mode, xmean, sigmam, sigma) +c +c description of parameters +c x - array of data points +c sigmax - array of standard deviations for data points +c npts - number of data points +c mode - determines method of weighting +c +1 (instrumental) weight(i) = 1./sigmax(i)**2 +c 0 (no weighting) weight(i) = 1. +c -1 (statistical) weight(i) = 1. +c xmean - weighted mean +c sigmam - standard deviation of mean +c sigma - standard deviation of data +c +c subroutines and function subprograms required +c none +c + subroutine xfit (x,sigmax,npts,mode,xmean,sigmam,sigma) + double precision sum,sumx,weight,free + dimension x(1),sigmax(1) +c +c accumulate weighted sums +c +11 sum=0. + sumx=0. + sigma=0. + sigmam=0. +20 do 32 i=1,npts +21 if (mode) 22,22,24 +22 weight=1. + goto 31 +24 weight=1./sigmax(i)**2 +31 sum=sum+weight +32 sumx=sumx+weight*x(i) +c +c evaluate mean and standard deviations +c +41 xmean=sumx/sum +51 do 52 i=1,npts +52 sigma=sigma+(x(i)-xmean)**2 + free=npts-1 +54 sigma=dsqrt(sigma/free) +61 if (mode) 62,64,66 +62 sigmam=dsqrt(xmean/sum) + goto 70 +64 sigmam=sigma/dsqrt(sum) + goto 70 +66 sigmam=dsqrt(1./sum) +70 return + end |