diff options
Diffstat (limited to 'noao/onedspec/fortran/polft1.f')
| -rw-r--r-- | noao/onedspec/fortran/polft1.f | 205 |
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 |