diff options
author | Joe Hunkeler <jhunkeler@gmail.com> | 2015-08-11 16:51:37 -0400 |
---|---|---|
committer | Joe Hunkeler <jhunkeler@gmail.com> | 2015-08-11 16:51:37 -0400 |
commit | 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 (patch) | |
tree | 4464880c571602d54f6ae114729bf62a89518057 /math/bevington/polfit.f | |
download | iraf-osx-40e5a5811c6ffce9b0974e93cdd927cbcf60c157.tar.gz |
Repatch (from linux) of OSX IRAF
Diffstat (limited to 'math/bevington/polfit.f')
-rw-r--r-- | math/bevington/polfit.f | 100 |
1 files changed, 100 insertions, 0 deletions
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 |