Initial commit

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