Initial commit

1 files changed, 128 insertions, 0 deletions

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