From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001
From: Joe Hunkeler <jhunkeler@gmail.com>
Date: Tue, 11 Aug 2015 16:51:37 -0400
Subject: Repatch (from linux) of OSX IRAF

---
 math/bevington/gridls.f | 102 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 102 insertions(+)
 create mode 100644 math/bevington/gridls.f

(limited to 'math/bevington/gridls.f')

diff --git a/math/bevington/gridls.f b/math/bevington/gridls.f
new file mode 100644
index 00000000..e62d3219
--- /dev/null
+++ b/math/bevington/gridls.f
@@ -0,0 +1,102 @@
+C SUBROUTINE GRIDLS.F 
+C
+C SOURCE
+C   BEVINGTON, PAGES 212-213.
+C
+C PURPOSE
+C   MAKE A GRID-SEARCH LEAST-SQUARES FIT TO DATA WITH A SPECIFIED
+C   FUNCTION WHICH IS NOT LINEAR IN COEFFICIENTS
+C
+C USAGE 
+C   CALL GRIDLS (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
+C COMMENTS
+C   DELTAA VALUES ARE MODIFIED BY THE PROGRAM
+C
+      SUBROUTINE GRIDLS (X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA,
+     *SIGMAA,YFIT,CHISQR)
+      DIMENSION X(1),Y(1),SIGMAY(1),A(1),DELTAA(1),SIGMAA(1), 
+     *YFIT(1)
+      REAL FUNCTN
+      EXTERNAL FUNCTN
+
+C
+11    NFREE=NPTS-NTERMS
+      FREE=NFREE
+      CHISQR=0.
+      IF (NFREE) 100,100,20
+20    DO 90 J=1,NTERMS
+C
+C EVALUATE CHI SQUARE AT FIRST TWO SEARCH POINTS
+C
+21    DO 22 I=1,NPTS
+22    YFIT(I)=FUNCTN (X,I,A)
+23    CHISQ1=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
+      FN=0.
+      DELTA=DELTAA(J) 
+41    A(J)=A(J)+DELTA 
+      DO 43 I=1,NPTS
+43    YFIT(I)=FUNCTN (X,I,A)
+44    CHISQ2=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
+45    IF (CHISQ1-CHISQ2) 51,41,61
+C
+C REVERSE DIRECTION OF SEARCH IF CHI SQUARE IS INCREASING
+C
+51    DELTA=-DELTA
+      A(J)=A(J)+DELTA 
+      DO 54 I=1,NPTS
+54    YFIT(I)=FUNCTN (X,I,A)
+      SAVE=CHISQ1
+      CHISQ1=CHISQ2
+57    CHISQ2=SAVE
+C
+C INCREMENT A(J) UNTIL CHI SQUARE INCREASES
+C
+61    FN=FN+1.
+      A(J)=A(J)+DELTA 
+      DO 64 I=1,NPTS
+64    YFIT(I)=FUNCTN (X,I,A)
+      CHISQ3=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
+66    IF (CHISQ3-CHISQ2) 71,81,81
+71    CHISQ1=CHISQ2
+      CHISQ2=CHISQ3
+      GOTO 61 
+C
+C FIND MINIMUM OF PARABOLA DEFINED BY LAST THREE POINTS 
+C
+81    DELTA=DELTA*(1./(1.+(CHISQ1-CHISQ2)/(CHISQ3-CHISQ2))+0.5)
+82    A(J)=A(J)-DELTA 
+83    SIGMAA(J)=DELTAA(J)*SQRT(2./(FREE*(CHISQ3-2.*CHISQ2+CHISQ1)))
+84    DELTAA(J)=DELTAA(J)*FN/3.
+90    CONTINUE
+C
+C EVALUATE FIT AND CHI SQUARE FOR FINAL PARAMETERS
+C
+91    DO 92 I=1,NPTS
+92    YFIT(I)=FUNCTN (X,I,A)
+93    CHISQR=FCHISQ (Y,SIGMAY,NPTS,NFREE,MODE,YFIT)
+100   RETURN
+      END
-- 
cgit