Initial commit

1 files changed, 113 insertions, 0 deletions

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