Initial commit

1 files changed, 205 insertions, 0 deletions

diff --git a/noao/onedspec/fortran/polft1.f b/noao/onedspec/fortran/polft1.f
new file mode 100644
index 00000000..625b69c7
--- /dev/null
+++ b/noao/onedspec/fortran/polft1.f
@@ -0,0 +1,205 @@
+C+
+C
+C  SUBROUTINE POLFT1
+C
+C  PURPOSE
+C	MAKE A LEAST SQUARES FIT TO DATA WITH A POLYNOMIAL CURVE
+C	Y = A(1) + A(2)*X + A(3)*X**2 + ...
+C
+C  USAGE
+C	CALL POLFIT(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR,ARR,IER)
+C
+C  DESCRIPTION OF PARAMETERS
+C	Y	- ARRAY OF DATA POINTS FOR DEPENDENT VARIABLE
+C	SIGMAY	- ARRAY OF STANDARD DEVIATIONS FOR Y-DATA POINTS
+C		  (OR - IN CASE MODE = 4 - WEIGHTS FOR POINTS)
+C	NPTS	- NO. OF PAIRS OF DATA POINTS
+C	NTERMS	- NO. OF COEFFICIENTS (DEGREE OF POLYNOMIAL + 1)
+C	MODE	- DETERMINES METHOD OF WEIGHTING LEAST SQUARES FIT
+C		  4 (USER DEFINED) WEIGHT(I) = SIGMAY(I)
+C		  3 (INSTRUMENTAL) WEIGHT(I) = 1./SIGMAY(I)**2
+C		  2 (NO WEIGHTING) WEIGHT(I) = 1.
+C		  1 (STATISTICAL)  WEIGHT(I) = 1./Y(I)
+C	A	- ARRAY OF COEFFICIENTS OF POLYNOMIAL
+C	CHISQR	- REDUCED CHISQUARE FOR FIT
+C	ARR	- DOUBLE PRECISION WORK BUFFER; MUST BE AT LEAST
+C		  400 WORDS LONG IN THE CALLING ROUTINE
+C	IER	- ERROR PARAMETER
+C		  -1 DET=0
+C		   0 O.K.
+C		  +1 NOT ENOUGH POINTS
+C
+C  SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
+C	DET(ARR,NORD) - EVALUATES THE DETERMINANT OF A SYMMETRIC
+C			TWO DIMENSIONAL MATRIX OF ORDER NORD
+C
+C  COMMENTS
+C	DIMENSION STATEMENT VALID FOR NTERMS UP TO 10
+C	FOR DETAILS OF ALGORITHM SEE "DATA REDUCTION AND ERROR
+C	ANALYSIS FOR THE PHYSICAL SCIENCES" - PH. R. BEVINGTON
+C	MC GRAW-HILL BOOK COMPANY
+C
+C	IN THIS SPECIAL VERSION THE ELEMENTS OF COEFFICIENT MATRIX
+C	ARR ARE NORMALIZED WITH RESPECT TO A VALUE COMPUTED VIA
+C	LOGARITHMIC INTERPOLATION BETWEEN THE SMALLEST AND LARGEST
+C	MATRIX ELEMENT
+C
+C  MODIFICATIONS BY A. SCHERMANN - INST. FOR ASTRONOMY,
+C	UNIV. OF VIENNA, AUSTRIA
+C
+C Modified to assume x-variable is index of Y array
+C-
+	SUBROUTINE POLFT1(Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR,ARR,IER)
+	DOUBLE PRECISION SUMX(19),SUMY(10),ARR(10,10),XTERM,YTERM,
+     1	CHISQ
+	DIMENSION Y(1),SIGMAY(1),A(1)
+	IER=0
+C
+C  CHECK DEGREE OF FREEDOM
+	FREE=NPTS-NTERMS
+	IF(FREE.GT.0.) GOTO 11
+	IER=1
+	GOTO 80
+C
+C  ACCUMULATE WEIGHTED SUMS
+11	NMAX=2*NTERMS-1
+	DO 13 N=1,NMAX
+13	SUMX(N)=0.
+	DO 15 J=1,NTERMS
+15	SUMY(J)=0.
+	CHISQ=0.
+21	DO 50 I=1,NPTS
+	XI=I
+	YI=Y(I)
+31	GOTO(32,37,38,39),MODE
+32	IF(YI)35,37,33
+33	WEIGHT=1./YI
+	GOTO 41
+35	WEIGHT=-1./YI
+	GOTO 41
+37	WEIGHT=1.
+	GOTO 41
+38	WEIGHT=1./SIGMAY(I)**2
+	GOTO 41
+39	WEIGHT=SIGMAY(I)
+41	XTERM=WEIGHT
+	DO 44 N=1,NMAX
+	SUMX(N)=SUMX(N)+XTERM
+44	XTERM=XTERM*XI
+45	YTERM=WEIGHT*YI
+	DO 48 N=1,NTERMS
+	SUMY(N)=SUMY(N)+YTERM
+48	YTERM=YTERM*XI
+49	CHISQ=CHISQ+WEIGHT*YI**2
+50	CONTINUE
+C
+C  GET LARGEST AND SMALLEST MATRIX ELEMENT (FOR NORMALIZATION)
+	XTERM=SUMX(1)
+	YTERM=XTERM
+	DO 100 I=2,NMAX
+	IF(SUMX(I).GT.XTERM) XTERM=SUMX(I)
+	IF(SUMX(I).LT.YTERM) YTERM=SUMX(I)
+100	CONTINUE
+	DO 110 I=1,NTERMS
+	IF(SUMY(I).GT.XTERM) XTERM=SUMY(I)
+	IF(SUMY(I).LT.YTERM) YTERM=SUMY(I)
+110	CONTINUE
+	IF(YTERM.LE.0.) YTERM=1.D0
+C
+C  LOGARITHMIC INTERPOLATION OF NORMALIZATION VALUE
+	XTERM=1.D1**((DLOG10(XTERM)+DLOG10(YTERM))/2.)
+C
+C  CONSTRUCT MATRICES AND CALCULATE COEFFICIENTS
+51	DO 54 J=1,NTERMS
+	DO 54 K=1,NTERMS
+	N=J+K-1
+54	ARR(J,K)=SUMX(N)/XTERM
+	DELTA=DET(ARR,NTERMS)
+	IF(DELTA.NE.0.) GOTO 61
+57	CHISQR=0.
+	DO 59 J=1,NTERMS
+59	A(J)=0.
+	IER=-1
+	GOTO 80
+61	DO 70 L=1,NTERMS
+62	DO 66 J=1,NTERMS
+	DO 65 K=1,NTERMS
+	N=J+K-1
+65	ARR(J,K)=SUMX(N)/XTERM
+66	ARR(J,L)=SUMY(J)/XTERM
+70	A(L)=DET(ARR,NTERMS)/DELTA
+C
+C  CALCULATE CHISQUARE
+71	DO 75 J=1,NTERMS
+	CHISQ=CHISQ-2.*A(J)*SUMY(J)
+	DO 75 K=1,NTERMS
+	N=J+K-1
+75	CHISQ=CHISQ+A(J)*A(K)*SUMX(N)
+77	CHISQR=CHISQ/FREE
+80	RETURN
+	END
+c-----------------------------------------------------------
+	FUNCTION POLYNO (COE,NPOL,IX)
+C
+C EVALUATE A POLYNOMIAL OF ORDER NPOL WITH COEFFICIENTS COE(1),...,
+C COE (NPOL+1)  FOR AN INDEX IX
+C
+	DIMENSION COE(*)
+C
+	IF(NPOL.GT.0) GO TO 10
+	POLYNO=COE(1)
+	RETURN
+C
+10	POLYNO=COE(NPOL+1)
+	DO 20 I=NPOL,1,-1
+20	    POLYNO=POLYNO*IX+COE(I)
+C
+	RETURN
+	END
+c-----------------------------------------------------------
+C+
+C
+C  FUNCTION DET
+C
+C  PURPOSE
+C	CALCULATE DETERMINANT OF A SQUARE MATRIX
+C
+C  USAGE
+C	DET = DET (ARR,NORDER)
+C
+C  DESCRIPTION OF PARAMETERS
+C	ARRAY	- MATRIX
+C	NORDER	- DEGREE OF MATRIX
+C
+C  SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
+C		NONE
+C
+C  COMMENTS
+C	THIS SUBPROGRAM DESTROYS THE INPUT MATRIX ARRAY
+C
+	FUNCTION DET(ARRAY,NORDER)
+	DOUBLE PRECISION ARRAY(10,10),SAVE
+	DET=1.
+	DO 90 K=1,NORDER
+C  INTERCHANGE COLUMNS, IF DIAGONAL ELEMENT IS ZERO
+	IF(ARRAY(K,K).NE.0.) GOTO 40
+	DO 10 J=K,NORDER
+	IF(ARRAY(K,J).NE.0.) GOTO 20
+10	CONTINUE
+	DET=0.
+	GOTO 100
+20	DO 30 I=K,NORDER
+	SAVE=ARRAY(I,J)
+	ARRAY(I,J)=ARRAY(I,K)
+30	ARRAY(I,K)=SAVE
+	DET=-DET
+C  SUBTRACT ROW K FROM LOWER ROWS TO GET DIAGONAL MATRIX
+40	DET=DET*ARRAY(K,K)
+	IF(K-NORDER.GE.0) GOTO 90
+	K1=K+1
+	DO 50 I=K1,NORDER
+	DO 50 J=K1,NORDER
+50	ARRAY(I,J)=ARRAY(I,J)-ARRAY(I,K)*ARRAY(K,J)/ARRAY(K,K)
+90	CONTINUE
+100	RETURN
+	END