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/chifit.f | |
download | iraf-osx-40e5a5811c6ffce9b0974e93cdd927cbcf60c157.tar.gz |
Repatch (from linux) of OSX IRAF
Diffstat (limited to 'math/bevington/chifit.f')
-rw-r--r-- | math/bevington/chifit.f | 169 |
1 files changed, 169 insertions, 0 deletions
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 |