C @(#)spec.for 17.1.1.1 (ES0-DMD) 01/25/02 17:55:52 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C.IDENTIFICATION: SPEC C.PURPOSE: Interactive line routines C.LANGUAGE: F77+ESOext C.AUTHOR: J.D.Ponz C.KEYWORDS: Line, centre C.ALGORITHM: Line center is found by: C maximum/minimum - position of the maximum/minimum value C gaussian - center of the fitted gaussian C gravity - gravity center of the 2 highest pixels C with respect to the third one C----------------------------------------------------------------- SUBROUTINE FIND(X,N1,N2,XSTR,XSTP,NL,IPL1,IPL2,IMODE,IMETH, + XCENTR,PEAK,IFAIL,W1,W2,ACOE,Y1,Y2) C C INTERMEDIATE ROUTINE C FIND CENTER C IMPLICIT NONE INTEGER N1,N2,NL,IPL1,IPL2,IMODE,IMETH,IFAIL,NPIX REAL XSTR,XSTP,XCENTR,PEAK,CHICHK,XGO REAL X(N1,N2),W1(1),W2(1),ACOE(4),Y1,Y2 C CHICHK = 0.005 XGO = XSTR + (IPL1-1)*XSTP NPIX = IPL2 - IPL1 + 1 Y1 = X(IPL1,NL) Y2 = X(IPL2,NL) IF (IMETH) 10,20,30 C C GAUSSIAN CENTER C 10 CALL SGAUS(X(IPL1,NL),W1,W2,IMODE,NPIX,IFAIL,XGO,XSTP,XCENTR, + CHICHK,0,PEAK,ACOE) RETURN C C GRAVITY C 20 CALL GRAVT(X(IPL1,NL),NPIX,IMODE,IFAIL,XGO,XSTP,XCENTR,PEAK) RETURN C C MAXUMUM C C 30 CALL CNTRH(X(IPL1,NL),NPIX,IMODE,IFAIL,XGO,XSTP,XCENTR,PEAK) ! Min/Max 30 CALL CNTRW(X(IPL1,NL),NPIX,IMODE,IFAIL,XGO,XSTP,XCENTR,PEAK) ! Moment RETURN END SUBROUTINE SFIND(X,XSTR,XSTP,IPL1,IPL2,IMODE,IMETH, + XCENTR,PEAK,IFAIL,W1,W2,ACOE,Y1,Y2) C++++ C. PURPOSE: Find center C--- IMPLICIT NONE REAL X(1) DOUBLE PRECISION XSTR,XSTP INTEGER IPL1,IPL2,IMODE,IMETH REAL XCENTR,PEAK,CHICHK,XSP,XGO INTEGER IFAIL,NPX REAL W1(1),W2(1),ACOE(4),Y1,Y2 C CHICHK = 0.005 XSP = REAL(XSTP) XGO = REAL(XSTR) + (IPL1-1)*XSP NPX = IPL2 - IPL1 + 1 Y1 = X(IPL1) Y2 = X(IPL2) IF (IMETH.LT.0) THEN ! gaussian center CALL SGAUS(X(IPL1),W1,W2,IMODE,NPX,IFAIL,XGO,XSP, + XCENTR,CHICHK,0,PEAK,ACOE) ELSE IF (IMETH.EQ.0) THEN CALL GRAVT(X(IPL1),NPX,IMODE,IFAIL,XGO,XSP,XCENTR,PEAK) ELSE CALL CNTRH(X(IPL1),NPX,IMODE,IFAIL,XGO,XSP,XCENTR,PEAK) ENDIF RETURN END SUBROUTINE GRAVT(WINDOW,NPIX,IWID,IWERR,XGO,DXSTP,XCNTR,AM) C C LINE CENTERING BY FINDING THE CENTER OF GRAVITY OF 2 HIGHEST C POINTS OF 3 POINTS C THE 2 HIGHEST POINTS HAVE VALUES RELATIVES TO THE THIRD C ONLY APPLICABLE IN EMISSION MODE C PARAMETERS : C WINDOW BUFFER WITH EXPECTED LINE CENTER IN MIDDLE C NPIX NUMBER OF SAMPLES IN WINDOW C IWID CODE 0 ABSORPTION, 1 EMISSION C IWERR ERROR RETURN, 0 OK, 1 ERROR C XGO XCOORDINATE FIRST ELEMENT IN WINDOW C DXSTP STEP IN X C XCNTR RETURN X CENTER CALCULATED C AM INTENSITY CENTER PIXEL C IMPLICIT NONE C REAL AM, WINDOW(1), AL, AH, DR INTEGER NPIX, IWID, IWERR, I, K C modified MP 900119 , C DOUBLE PRECISION XGO, DXSTP, XCNTR, A, B, XSHIFT REAL XGO, DXSTP, XCNTR, A, B, XSHIFT C IF (IWID.NE.1) GO TO 20 C C FIND MAXIMUM C AM = WINDOW(1) I = 1 DO 10 K = 2,NPIX IF (WINDOW(K).GT.AM) THEN AM = WINDOW(K) I = K END IF 10 CONTINUE C C CHECK NOT BOUNDARY C IF (I.EQ.1 .OR. I.EQ.NPIX) GO TO 20 XCNTR = XGO + (I-1)*DXSTP C C FIND LOWEST OF THE FLANKING PIXELS C AL = WINDOW(I-1) AH = WINDOW(I+1) DR = 1. IF (AL.GE.AH) THEN AL = WINDOW(I+1) AH = WINDOW(I-1) DR = -1. END IF AM = WINDOW(I) A = AM - AL B = AH - AL XSHIFT = (B/ (A+B))*DXSTP XCNTR = XCNTR + XSHIFT*DR IWERR = 0 RETURN C C ERROR RETURN C 20 IWERR = 1 RETURN END SUBROUTINE CNTRW (WINDOW,NPIX,IWID,IWERR,XGO,DXSTP,XCNTR,AM) C C CENTER OF LINE ESTIMATED BY FIRST MOMENT OF THE VALUES C C PARAMETERS : C WINDOW BUFFER WITH EXPECTED LINE CENTER IN MIDDLE C NPIX NUMBER OF SAMPLES IN WINDOW C IWID CODE 0 ABSORPTION, 1 EMISSION C IWERR ERROR RETURN, 0 OK, 1 ERROR C XGO XCOORDINATE FIRST ELEMENT IN WINDOW C DXSTP STEP IN X C XCNTR RETURN X CENTER CALCULATED C AM INTENSITY CENTER PIXEL C IMPLICIT NONE C REAL WINDOW(1), AM INTEGER NPIX, IWID, IWERR REAL XGO, DXSTP, XCNTR C INTEGER IPN, I REAL SUMY, SUMXY,XMOM C IWERR = 0 AM = WINDOW(1) IPN = 1 SUMY = 0. SUMXY = 0. DO 10 I = 1,NPIX SUMXY = SUMXY + I*WINDOW(I) SUMY = SUMY + WINDOW(I) 10 CONTINUE XMOM = SUMXY / SUMY C CHECK RESULT C IF (XMOM.LE.1 .OR. XMOM.GE.NPIX) GO TO 30 XCNTR = XGO + (XMOM-1)*DXSTP RETURN C C ERROR RETURN C 30 IWERR = 1 RETURN END SUBROUTINE CNTRH(WINDOW,NPIX,IWID,IWERR,XGO,DXSTP,XCNTR,AM) C C CENTER OF LINE FOUND BY SIMPLE MAXIMUM WITHIN A WINDOW C C PARAMETERS : C WINDOW BUFFER WITH EXPECTED LINE CENTER IN MIDDLE C NPIX NUMBER OF SAMPLES IN WINDOW C IWID CODE 0 ABSORPTION, 1 EMISSION C IWERR ERROR RETURN, 0 OK, 1 ERROR C XGO XCOORDINATE FIRST ELEMENT IN WINDOW C DXSTP STEP IN X C XCNTR RETURN X CENTER CALCULATED C AM INTENSITY CENTER PIXEL C IMPLICIT NONE C REAL WINDOW(1), AM INTEGER NPIX, IWID, IWERR REAL XGO, DXSTP, XCNTR C INTEGER IPN, I C IWERR = 0 AM = WINDOW(1) IPN = 1 IF (IWID.NE.1) THEN C C FIND MINIMUM C DO 10 I = 2,NPIX IF (WINDOW(I).LT.AM) THEN AM = WINDOW(I) IPN = I END IF 10 CONTINUE ELSE C C FIND MAXIMUM C DO 20 I = 2,NPIX IF (WINDOW(I).GT.AM) THEN AM = WINDOW(I) IPN = I END IF 20 CONTINUE END IF C C CHECK RESULT C IF (IPN.EQ.1 .OR. IPN.EQ.NPIX) GO TO 30 XCNTR = XGO + (IPN-1)*DXSTP RETURN C C ERROR RETURN C 30 IWERR = 1 RETURN END SUBROUTINE SGAUS(WINDOW,XBUF,YFIT,IWID,NPIX,IWERR,XGO,DXSTP, 2 XCNTR,CHICHK,IOUTPT,AM,A) C C LINE CENTER FINDING WITH GAUSSIAN FIT C C PARAMETERS : C WINDOW BUFFER WITH EXPECTED LINE CENTER IN MIDDLE C XBUF WORK BUFFER C YFIT WORK BUFFER C IWID CODE 0 ABSORPTION, 1 EMISSION C NPIX NUMBER OF SAMPLES IN WINDOW C IWERR ERROR RETURN, 0 OK, 1 ERROR C XGO XCOORDINATE FIRST ELEMENT IN WINDOW C DXSTP STEP IN X C XCNTR RETURN X CENTER CALCULATED C CHICHK CHECK CHISQUARE C IOUTPT DISPLAY INTERMEDIATE RESULTS 0 NO, 1 YES C AM INTENSITY CENTER PIXEL C A ARRAY(4) COEFFICIENTS C INTEGER I,K1,K2,NPIX,IWID,IWERR,IOUTPT INTEGER ITCNT,ITCHK,JER,IBRK REAL WINDOW(*),XBUF(*),YFIT(*),A(*) REAL DELTAA(4),SIGMAA(4),XCNTR,XGO,CHICHK REAL DXSTP,CHISQR,FLAMDA,OLDCH,DIF,HALF,WIDTH REAL AM,X,STOP,TOPM CHARACTER*70 IHDD DATA IHDD/' CHISQ FWHM CENTER'/ C C BUFFER WITH X VALUES C X = XGO DO 10 I = 1,NPIX XBUF(I) = XGO + (I-1)*DXSTP 10 CONTINUE C C FIND INITIAL GUESS C CALL CNTRH(WINDOW,NPIX,IWID,IWERR,XGO,DXSTP,A(2),AM) IF (IWERR.NE.0) GO TO 120 A(4) = (WINDOW(1)+WINDOW(NPIX))/2. A(1) = AM - A(4) C C CALCULATE HALF WIDTH C HALF = A(1)/2. + A(4) C C CALCULATION OF FWHM DEPENDS ON SIGN OF THE LINE C IF (IWID.NE.1) THEN DO 20 I = 1,NPIX IF (WINDOW(I).LT.HALF) GO TO 30 20 CONTINUE 30 K1 = I DO 40 I = K1,NPIX IF (WINDOW(I).GT.HALF) GO TO 50 40 CONTINUE 50 K2 = I ELSE DO 60 I = 1,NPIX IF (WINDOW(I).GT.HALF) GO TO 70 60 CONTINUE 70 K1 = I DO 80 I = K1,NPIX IF (WINDOW(I).LT.HALF) GO TO 90 80 CONTINUE 90 K2 = I END IF WIDTH = ABS(FLOAT(K2-K1)*DXSTP) A(3) = WIDTH/2.354 OLDCH = 9.E16 C C OPTIONAL OUTPUT C C IF (IOUTPT .NE. 0) CALL WRUSER(0,IHDD,' ',' ',ISTAT) C C CALL FITTING ROUTINE C ITERATE 50 TIMES C ITCNT = 0 ITCHK = 50 100 FLAMDA = 0.001 IBRK = 0 C CALL CURFI(XBUF,WINDOW,SIGMAA,NPIX,4,0,A,DELTAA,FLAMDA,YFIT, + CHISQR,JER) IF (JER.NE.0) GO TO 120 C C CALC AND DISPLAY INTERMEDIATE RESULTS C IF (IBRK.EQ.-1) GO TO 110 C IF (IOUTPT.NE.0) THEN C TYPE *,NPIX C WRITE(LABL, 100) CHISQR,A(3)*2.345,A(2) C CALL STTPUT(LABL,ISTAT) C ENDIF C C SEE IF CHISQR HAS CHANGED SIGNIFICANTLY C DIF = OLDCH - CHISQR STOP = DIF/CHISQR OLDCH = CHISQR ITCNT = ITCNT + 1 IF (ITCNT.GT.ITCHK) GO TO 120 IF (STOP.GT.CHICHK) GO TO 100 110 XCNTR = A(2) C C LAST CHECK THAT WE ARE WITHIN THE WINDOW C TOPM = XBUF(NPIX) IF (DXSTP.LT.0.) THEN IF ((XCNTR.GT.XGO) .OR. (XCNTR.LT.TOPM)) GO TO 120 ELSE IF ((XCNTR.LT.XGO) .OR. (XCNTR.GT.TOPM)) GO TO 120 END IF IWERR = 0 A(3) = A(3)*2.345 RETURN C C ERROR RETURN C 120 IWERR = 1 RETURN C9000 FORMAT (E10.3,E10.3,F9.3) END SUBROUTINE CURFI(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,DELTAA,FLAMDA,YFIT, + CHISQR,IERR) C C LEAST SQUARES FIT TO A NON-LINEAR FUNCTION C C PARAMETERS: C X ARRAY OF DATA POINTS IND. VAR. C Y ARRAY OF DATA POINTS DEP. VAR. C SIGMAY ARRAY OF STD DEV FOR Y DATA POINTS C NPTS NUMBER OF PAIRS OF DATA POINTS C MODE DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT C +1 (INSTRUMENTAL) WEIGHT = 1./SIGMAY(I)**2 C 0 (NO WEIGHTING) WEIGHT = 1. C -1 (STATISTICAL ) WEIGHT = 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 IERR ERROR RETURN 0 OK, 1 NOT CONVERGING C C DIMENSION STATEMENT VALID FOR NTERMS UP TO 10 C SET FLAMDA TO 0.001 AT THE BEGINNING OF THE SEARCH C INTEGER I,J,K,NPTS,NFREE ,NTERMS,IERR INTEGER MODE ,ICNT REAL X(*),Y(*),SIGMAY(*),DELTAA(*),YFIT(*) REAL A(*),DET,FLAMDA,CHISF,FUNCT REAL WEIGHT(400),ALPHA(10,10),BETA(10) REAL DERIV(10),B(10),CHISQR ,CHISQ1 DOUBLE PRECISION ARRAY(10,10) EXTERNAL FUNCT C ICNT = 0 IERR = 1 NFREE = NPTS - NTERMS IF (NFREE) 20,20,30 20 CHISQR = 0. RETURN C C EVALUATE WEIGHTS C 30 CONTINUE IERR = 0 DO 100 I = 1,NPTS IF (MODE) 40,70,80 40 IF (Y(I)) 60,70,50 50 WEIGHT(I) = 1./Y(I) GO TO 90 60 WEIGHT(I) = 1./ (-Y(I)) GO TO 90 70 WEIGHT(I) = 1. GO TO 90 80 WEIGHT(I) = 1./SIGMAY(I)**2 90 CONTINUE 100 CONTINUE C C EVALUATE ALPHA AND BETA MATRICES C DO 120 J = 1,NTERMS BETA(J) = 0. DO 110 K = 1,J ALPHA(J,K) = 0. 110 CONTINUE 120 CONTINUE DO 150 I = 1,NPTS CALL FDERI(X,I,A,DELTAA,NTERMS,DERIV) DO 140 J = 1,NTERMS BETA(J) = BETA(J) + WEIGHT(I)* (Y(I)-FUNCT(X,I,A))* + DERIV(J) DO 130 K = 1,J ALPHA(J,K) = ALPHA(J,K) + WEIGHT(I)*DERIV(J)*DERIV(K) 130 CONTINUE 140 CONTINUE 150 CONTINUE DO 170 J = 1,NTERMS DO 160 K = 1,J ALPHA(K,J) = ALPHA(J,K) 160 CONTINUE 170 CONTINUE C C EVALUATE CHI SQUARE AT STARTING POINT C DO 180 I = 1,NPTS YFIT(I) = FUNCT(X,I,A) 180 CONTINUE CHISQ1 = CHISF(Y,SIGMAY,NPTS,NFREE,MODE,YFIT) C C INVERT MATRIX C 190 DO 210 J = 1,NTERMS DO 200 K = 1,NTERMS IF (ABS(ALPHA(J,J)).LT.1.E-10.OR. . ABS(ALPHA(K,K)).LT.1.E-10) GOTO 987 ARRAY(J,K) = ALPHA(J,K)/SQRT(ALPHA(J,J)*ALPHA(K,K)) 200 CONTINUE ARRAY(J,J) = 1. + FLAMDA 210 CONTINUE CALL INVMAT(ARRAY,NTERMS,DET) DO 230 J = 1,NTERMS B(J) = A(J) DO 220 K = 1,NTERMS B(J) = B(J) + BETA(K)*ARRAY(J,K)/ + SQRT(ALPHA(J,J)*ALPHA(K,K)) 220 CONTINUE 230 CONTINUE C C IF CHI SQUARE INCREASED, INCREASE FLAMDA AND TRY AGAIN C DO 240 I = 1,NPTS YFIT(I) = FUNCT(X,I,B) 240 CONTINUE CHISQR = CHISF(Y,SIGMAY,NPTS,NFREE,MODE,YFIT) IF (CHISQ1-CHISQR) 250,270,270 C C CHECK LOOPS C 250 ICNT = ICNT + 1 IF (ICNT.LT.60) GO TO 260 987 IERR = 1 RETURN 260 FLAMDA = 10.*FLAMDA GO TO 190 C C EVALUATE PARAMETERS C 270 DO 280 J = 1,NTERMS A(J) = B(J) 280 CONTINUE FLAMDA = FLAMDA/10. RETURN END REAL FUNCTION FUNCT(X,I,A) C C EVALUATE TERMS OF FUNCTION FOR NON-LINEAR LEAST-SQUARES SEARCH C WITH FORM OF GAUSSIAN PEAK C C PARAMETERS: C X ARRAY OF DATA POINTS C I INDEX OF DATA POINTS C A ARRAY OF PARAMETERS C REAL X(*),A(*),Z,XI,Z2 INTEGER I XI = X(I) FUNCT = A(4) Z = (XI-A(2))/A(3) Z2 = Z*Z IF (Z2-50.) 40,50,50 40 FUNCT = FUNCT + A(1)*EXP(-Z2*0.5) 50 RETURN END FUNCTION CHISF(Y,SIGMAY,NPTS,NFREE,MODE,YFIT) C C EVALUATE REDUCED CHI SQUARE FOR FIT TO DATA C CHISF = SUM((Y-YFIT)**2/SIGMA**2)/NFREE C C PARAMETERS: C Y ARRAY OF DATA POINTS C SIGMAY ARRAY OF STANDARD DEVIATIONS C NPTS NUMBER OF DATA POINTS C NFREE NUMBER OF DEGREES OF FREEDOM C MODE DETERMINES METHOD OF WEIGHTING LEAST-SQUARES FIT C +1 (INSTRUMENTAL) WEIGHT = 1./SIGMAY(I)**2 C 0 (NO WEIGHTING) WEIGHT = 1. C -1 (STATISTICAL ) WEIGHT = 1./Y(I) C YFIT ARRAY OF CALCULATED VALUES OF Y C INTEGER NFREE,MODE,NPTS,I REAL Y(*),SIGMAY(*),YFIT(*),CHISF DOUBLE PRECISION CHISQ,WEIGHT C CHISQ = 0. IF (NFREE) 30,30,40 30 CHISF = 0. RETURN C C ACCUMULATE CHI SQUARE C 40 DO 110 I = 1,NPTS IF (MODE) 50,80,90 50 IF (Y(I)) 70,80,60 60 WEIGHT = 1./Y(I) GO TO 100 70 WEIGHT = 1./ (-Y(I)) GO TO 100 80 WEIGHT = 1. GO TO 100 90 WEIGHT = 1./SIGMAY(I)**2 100 CHISQ = CHISQ + WEIGHT* (Y(I)-YFIT(I))**2 110 CONTINUE C C DIVIDE BY NUMBER OF DEGREES OF FREEDOM C CHISF = CHISQ/NFREE RETURN END SUBROUTINE FDERI(X,I,A,DELTAA,NTERMS,DERIV) C C EVALUATES DERIVATIVES OF FUNCTION FOR LEAST SQUARES SEARCH C WITH FORM OF GAUSSIAN PEAK C C PARAMETERS: C X ARRAY OF DATA POINTS OF INDEPENDENT VARIABLE C I INDEX OF DATA POINTS C A ARRAY OF PARAMETERS C DELTAA ARRAY OF PARAMETER INCREMENTS C NTERMS NUMBER OF PARAMETERS C DERIV DERIVATIVES OF FUNCTION C INTEGER I,J,NTERMS REAL X(*),A(*),DELTAA(*),DERIV(*) REAL Z,Z2,XI XI = X(I) Z = (XI-A(2))/A(3) Z2 = Z*Z IF (Z2-50.) 40,20,20 20 DO 30 J = 1,3 DERIV(J) = 0. 30 CONTINUE GO TO 50 C C ANALYTICAL EXPRESSION FOR DERIVATIVES C 40 DERIV(1) = EXP(-Z2/2.) DERIV(2) = A(1)*DERIV(1)*Z/A(3) DERIV(3) = DERIV(2)*Z 50 DERIV(4) = 1. RETURN END SUBROUTINE INVMAT(ARR,N,DET) C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 0:26 - 4 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C INVERTS A SYMMETRIC MATRIX AND CALCULATES ITS DETERMINANT C USING A GAUSS-JORDAN ELIMINATION METHOD. THE DEGREE N OF C THE MATRIX MUST NOT BE GREATER THAN 10. C THE INVERSION IS IN DOUBLE PRECISION, AND C THE ARRAY ARR MUST BE D-P. THE DETERMINANT IS ALLWAYS REAL. C C PARAMS : C ARR : ARRAY CONTAINING THE (N*N) MATRIX TO BE INVERTED. C THE INVERTED MATRIX IS RETURNED IN ARR. C N : DEGREE OF MATRIX C DET : DETERMINANT OF INPUT-MATRIX C C IMPLICIT NONE REAL DET DOUBLE PRECISION ARR(10,10),AMAX,SAVE,DDET INTEGER IK(10),JK(10),I,J,K,N,L C DET = 1.0 DDET = 1.0D0 DO 210 K = 1,N AMAX = 0.0D0 C C FIND LARGEST ELEMENT ARR(I,J) IN REST OF MATRIX C 10 DO 40 I = K,N DO 30 J = K,N IF (DABS(AMAX)-DABS(ARR(I,J))) 20,20,30 20 AMAX = ARR(I,J) IK(K) = I JK(K) = J 30 CONTINUE 40 CONTINUE C ---------------------------------------------------- C INTERCHANGE ROWS AND COLUMNS TO PUT AMAX IN ARR(K,K) C ---------------------------------------------------- IF (AMAX) 60,50,60 50 DET = 0. GO TO 280 60 I = IK(K) IF (I-K) 10,90,70 70 DO 80 J = 1,N SAVE = ARR(K,J) ARR(K,J) = ARR(I,J) ARR(I,J) = -SAVE 80 CONTINUE C 90 J = JK(K) IF (J-K) 10,120,100 100 DO 110 I = 1,N SAVE = ARR(I,K) ARR(I,K) = ARR(I,J) ARR(I,J) = -SAVE 110 CONTINUE C ------------------------------------- C ACCUMULATE ELEMENTS OF INVERSE MATRIX C ------------------------------------- 120 DO 140 I = 1,N IF (I-K) 130,140,130 130 ARR(I,K) = -ARR(I,K)/AMAX 140 CONTINUE C DO 180 I = 1,N DO 170 J = 1,N IF (I-K) 150,170,150 150 IF (J-K) 160,170,160 160 ARR(I,J) = ARR(I,J) + ARR(I,K)*ARR(K,J) 170 CONTINUE 180 CONTINUE C DO 200 J = 1,N IF (J-K) 190,200,190 190 ARR(K,J) = ARR(K,J)/AMAX 200 CONTINUE ARR(K,K) = 1.0/AMAX DET = DET*AMAX DDET = DDET*AMAX C 210 CONTINUE C -------------------------- C RESTORE ORDERING OF MATRIX C -------------------------- DO 270 L = 1,N K = N - L + 1 J = IK(K) IF (J-K) 240,240,220 220 DO 230 I = 1,N SAVE = ARR(I,K) ARR(I,K) = -ARR(I,J) ARR(I,J) = SAVE 230 CONTINUE C 240 I = JK(K) IF (I-K) 270,270,250 250 DO 260 J = 1,N SAVE = ARR(K,J) ARR(K,J) = -ARR(I,J) ARR(I,J) = SAVE 260 CONTINUE C 270 CONTINUE DET = DDET RETURN C ------------------- C DETERMINANT IS ZERO C ------------------- 280 DET = 0.0 RETURN END