C @(#)spehpcntr.for 13.1.1.1 (ES0-DMD) 06/02/98 18:18:12 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: SPEHPCNTR C.PURPOSE: Interactive line centering. output stored in a table 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.INPUT/OUTPUT: Interactive input via the cursor C table columns : C :XSTART - first cursor input C :XEND - second cursor input C :XCEN - computed center C :YSTART - image value for first cursor input C :YEND - image value for second cursor input C :PEAK - maximum value C :YFIT - fitted maximum in gaussian method only C :FWHM - fwhm in gaussian method only C.VERSION: 840328 JDP Creation C.VERSION: 870115 JDP ?? C.VERSION: 880703 RHW ESO-FORTRAN Conversion C.VERSION: 900430 RHW Map plotted line separately C----------------------------------------------------------------- PROGRAM SPCNTR IMPLICIT NONE INTEGER MADRID INTEGER I1,I2,IAC,IACT,IAV,IFAIL,II,II1,IMETH,NPTS INTEGER IMODE,IPL1,IPL2,ISEQ,ISORT,ISTAT INTEGER L,NAXIS,NCOLS,NCOUNT,NNCOL,NPL,NL INTEGER NVAL, NSIZE(2) INTEGER*8 PNTRI, PNTRW, PNTRT INTEGER IMNOI, IMNOT, IMNOW INTEGER IMAGE(4) INTEGER NPIX(2),NN INTEGER KEY INTEGER APPFLG INTEGER STATUS,IC(8) INTEGER KUN, KNUL, TID,ACOL, AROW REAL SXY,X1,X2,XYS REAL XCUR,YCUR,PXVL REAL VALUE(8),W1(1024),W2(1024),ACOE(4) DOUBLE PRECISION START(2),STEP(2) DOUBLE PRECISION DXYS, DSXY CHARACTER NAME*20,INFO*20,METH*2 CHARACTER TABLE*64,TABU(8)*16,TABL(8)*16,IEA*1 CHARACTER CLINE*72,HEAD1*72,HEAD2*72,FORM(8)*6,PROMPT*72,IAPP*1 CHARACTER IDENT*72,CUNIT*48,COML*80 INCLUDE 'MID_INCLUDE:ST_DEF.INC/NOLIST' INCLUDE 'MID_INCLUDE:PLTDEC.INC/NOLIST' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC/NOLIST' C DATA NCOLS/6/ DATA TABU/' ',' ',' ',' ',' ',' ',' ',' '/ DATA TABL/'START','END','CENTER','INT_START','INT_END', 2 'INT_PEAK','INT_FIT','FWHM'/ DATA FORM/'F10.3','F10.3','F10.3','G15.5','G15.5','G15.5', 2 'G15.5','F10.3'/ DATA HEAD1/ 2 ' start end center pixel_value FWHM'/ DATA HEAD2/ 2 ' start end center pixel_value'/ DATA PROMPT/' Ready for cursor input ...'/ C 9000 FORMAT(3F13.3,G13.5,G13.5) C9010 FORMAT(4G13.5) C9020 FORMAT(A) C C *** start of executable code *************************************** CALL STSPRO('SPCNTR') ! initialize CALL STKRDC('PLCURSOR',1,1,20,IACT,INFO,KUN,KNUL,STATUS) L = INDEX(INFO//' ',' ') - 1 NAME = INFO(1:L) C CALL STKRDC('P3',1,1,1,IACT,IEA,KUN,KNUL,STATUS) CALL UPCAS(IEA,IEA) IF (IEA.EQ.'E') THEN IMODE = 1 ELSE IMODE = 0 END IF C CALL STKRDC('MID$CMND',1,1,80,IACT,COML,KUN,KNUL,STATUS) METH = COML(11:12) IF (METH.EQ.'GA') THEN IMETH = -1 NCOLS = 8 ELSE IF (METH.EQ.'GR') THEN IMETH = 0 ELSE IMETH = 1 END IF END IF C CALL STKRDC('P4',1,1,1,IACT,IAPP,KUN,KNUL,STATUS) CALL UPCAS(IAPP,IAPP) IF (IAPP.EQ.'A') THEN APPFLG = 1 ELSE APPFLG = 0 END IF ISEQ = 0 C C *** read image CALL STFOPN(NAME,D_R4_FORMAT,0,F_IMA_TYPE,IMNOI,ISTAT) CALL STDRDI(IMNOI,'NAXIS',1,1,IAC,NAXIS,KUN,KNUL,ISTAT) CALL STDRDI(IMNOI,'NPIX',1,NAXIS,IAC,NPIX,KUN,KNUL,ISTAT) CALL STDRDD(IMNOI,'START',1,NAXIS,IAC,START,KUN,KNUL,ISTAT) CALL STDRDD(IMNOI,'STEP',1,NAXIS,IAC,STEP,KUN,KNUL,ISTAT) CALL STDRDC(IMNOI,'IDENT',1,1,72,IAC,IDENT,KUN,KNUL,ISTAT) CALL STDRDC(IMNOI,'CUNIT',1,1,64,IAC,CUNIT,KUN,KNUL,ISTAT) C NPL = NPIX(1) NL = NPIX(2) II1 = 16*NAXIS C C *** allocate virtual memory for a single line of the frame CALL STKRDI('PLISTAT',11,4,IAC,IMAGE,KUN,KNUL,ISTAT) IF (IMAGE(3).EQ.IMAGE(4)) THEN II = (IMAGE(3)-1)*NPL + 1 CALL STFMAP(IMNOI,F_I_MODE,II,NPL,IAC,PNTRW,ISTAT) SXY = STEP(1) XYS = START(1) DSXY = STEP(1) DXYS = START(1) ELSE C ! here for the columns NN = NPIX(1)*NPIX(2) CALL STFMAP(IMNOI,F_I_MODE,1,NN,IAC,PNTRI,ISTAT) CALL STFCRE('MIDTRAN',D_R4_FORMAT,F_X_MODE,F_IMA_TYPE, 2 NVAL,IMNOT,ISTAT) CALL STFMAP(IMNOT,F_X_MODE,1,NN,IAV,PNTRT,ISTAT) NSIZE(1) = 128 NSIZE(2) = 256 CALL LINCOL(MADRID(PNTRI),NPIX,NSIZE,MADRID(PNTRT)) ! transpose frame CALL STFUNM(IMNOI,ISTAT) C IF (NVAL.LT.512) NVAL = 512 CALL STFCRE('PLOTWORK', D_R4_FORMAT,F_X_MODE, F_IMA_TYPE, 2 NVAL,IMNOW,ISTAT) ! create scratch frame CALL STFMAP(IMNOW,F_X_MODE,1,NVAL,IAC,PNTRW,ISTAT) II = (IMAGE(1)-1)*NL PNTRT = PNTRT + II CALL COPYF(MADRID(PNTRT),MADRID(PNTRW),NVAL) ! copy right column CALL STFUNM(IMNOT,ISTAT) SXY = STEP(2) XYS = START(2) DSXY = STEP(2) DXYS = START(2) ENDIF C TABU(1) = CUNIT(II1+1:II1+16) TABU(2) = TABU(1) TABU(3) = TABU(1) TABU(4) = CUNIT(1:16) TABU(5) = TABU(4) TABU(6) = TABU(4) TABU(7) = TABU(4) TABU(8) = TABU(1) C C *** get table name for storage of data CALL STKRDC('P2',1,1,64,IACT,TABLE,KUN,KNUL,STATUS) IF (APPFLG.EQ.0) THEN CALL TBTINI(TABLE,F_TRANS,F_O_MODE,10,200,TID,STATUS) DO 10 IAV = 1,NCOLS CALL TBCINI(TID,D_R4_FORMAT,1,FORM(IAV),TABU(IAV), + TABL(IAV),IC(IAV),STATUS) 10 CONTINUE ELSE CALL TBTOPN(TABLE,F_U_MODE,TID,STATUS) CALL TBIGET(TID,NNCOL,ISEQ,ISORT,ACOL,AROW,STATUS) DO 20 IAV = 1,NCOLS CALL TBLSER(TID,TABL(IAV),IC(IAV),STATUS) IF (IC(IAV).EQ.-1) THEN CALL STTPUT('*** FATAL: Wrong input table ',STATUS) CALL TBTCLO(TID,STATUS) CALL STSEPI END IF 20 CONTINUE END IF C C *** restore the graphics display C CALL GCINI CALL AGSSET('USER') LTYPE = 2 C C *** loop on input - pairs of points IF (IMETH.EQ.-1) THEN CALL STTPUT(HEAD1,STATUS) ELSE CALL STTPUT(HEAD2,STATUS) END IF CALL STTPUT(' ',STATUS) NCOUNT = 0 C 30 CONTINUE CALL AGVLOC(XCUR,YCUR,KEY,PXVL) IF ((KEY.EQ.13) .OR. (KEY.EQ.11)) THEN CALL STTPUT('*** WARNING: do NOT use the return key',ISTAT) GO TO 30 ELSE IF (KEY.EQ.32) THEN CALL TBTCLO(TID,STATUS) CALL PTCLOS CALL STSEPI ELSE CALL AGGPLM(XCUR,YCUR,1,4) X1 = XCUR VALUE(1) = XCUR IPL1 = NINT((X1-XYS)/SXY) + 1 END IF C 40 CONTINUE CALL AGVLOC(XCUR,YCUR,KEY,PXVL) IF ((KEY.EQ.13) .OR. (KEY.EQ.11)) THEN CALL STTPUT('*** WARNING: do NOT use the return key',ISTAT) GO TO 40 ELSE IF (KEY.EQ.32) THEN CALL TBTCLO(TID,STATUS) CALL PTCLOS CALL STSEPI ELSE CALL AGGPLM(XCUR,YCUR,1,4) X2 = XCUR VALUE(2) = XCUR IPL2 = NINT((X2-XYS)/SXY) + 1 ENDIF C C *** find center I1 = MIN(IPL1,IPL2) I2 = MAX(IPL1,IPL2) CALL FIND(MADRID(PNTRW),DXYS,DSXY,I1,I2,IMODE,IMETH, 2 VALUE(3),VALUE(6),IFAIL,W1,W2,ACOE,VALUE(4),VALUE(5)) VALUE(7) = ACOE(1) VALUE(8) = ACOE(3) C IF (IFAIL.EQ.0) THEN ISEQ = ISEQ + 1 CALL TBRWRR(TID,ISEQ,NCOLS,IC,VALUE,STATUS) IF (IMETH.EQ.-1) THEN WRITE (CLINE,9000) VALUE(1),VALUE(2),VALUE(3),VALUE(6), + VALUE(8) ELSE WRITE (CLINE,9000) VALUE(1),VALUE(2),VALUE(3),VALUE(6) END IF C CALL STTPUT(CLINE,STATUS) NCOUNT = NCOUNT + 1 IF (IMETH.EQ.-1) THEN NPTS = I2 - I1 + 1 IF( LTYPE .GT. 0 ) CALL AGGPLL(W1,W2,NPTS) END IF END IF C GO TO 30 END SUBROUTINE FIND(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 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(1),XBUF(1),YFIT(1),A(1) 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(1),Y(1),SIGMAY(1),DELTAA(1),YFIT(1) REAL A(1),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(1),A(1),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(1),SIGMAY(1),YFIT(1),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(1),A(1),DELTAA(1),DERIV(1) 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