C @(#)tsawdt.for 17.1.1.1 (ESO-DMD) 01/25/02 17:20:37 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 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence 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 (c) 1992 European Southern Observatory and Warsaw Observatory C.IDENT tsawdt.for C.AUTHOR Alex Schwarzenberg-Czerny, ESO and Warsaw Observatory C.KEYWORD MIDAS, time series analysis, WIDTH/TIME C.LANGUAGE FORTRAN 77 C.PURPOSE Evaluate width of a spectral line C.VERSION 0.0 June 1992 C.RETURNS None C.ENVIRON TSA context C C 010625 last modif C C----------------------------------------------------------------------------- C C PROGRAM TSAWDT C INCLUDE 'MID_REL_INCL:TSA_DEF.INC' INCLUDE 'MID_INCLUDE:ST_DEF.INC' C INTEGER NAXIS,NVALUES PARAMETER (NAXIS=2,NVALUES=5) INTEGER NSTEPS(NAXIS),IBUF(5) INTEGER IACTS,KNUL,KUN,IFLAG,ITYP,IAXIS INTEGER IWORK,IPER INTEGER*8 PPER,PWORK C DOUBLE PRECISION BANDW ! INPUT BAND WIDTH DOUBLE PRECISION BANDC ! INPUT BAND CENTRE DOUBLE PRECISION START(NAXIS),STEP(NAXIS) DOUBLE PRECISION WIDTH(NVALUES,2) C CHARACTER*60 IDENT CHARACTER*60 INAME ! NAME OF INPUT FRAME CHARACTER*64 CUNIT C C INCLUDE 'MID_REL_INCL:TSA_DAT.INC' INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C C Get parameters C CALL STSPRO ('tsawdt') CALL STKRDC ('IN_A', 1,1,60,IACTS,INAME,KUN,KNUL,ISTAT) CALL STKRDD ('WIDTH', 1,1, IACTS,BANDW,KUN,KNUL,ISTAT) CALL STKRDD ('CENTRE', 1,1, IACTS,BANDC,KUN,KNUL,ISTAT) C C Map input data C IFLAG=2 CALL STFINF (INAME,IFLAG,IBUF,ISTAT) ITYP=IBUF(2) IF (ITYP.NE.D_R8_FORMAT) THEN CALL STETER(1,'Input spectrum must be in DOUBLE PRECISION') ENDIF CALL STIGET(INAME,D_OLD_FORMAT,F_I_MODE,F_IMA_TYPE,NAXIS, $ IAXIS,NSTEPS,START,STEP,IDENT,CUNIT,PPER,IPER,ISTAT) C C Create work table C CALL STFCRE('ZZMIDWORK',D_R8_FORMAT,F_X_MODE,F_IMA_TYPE,NSTEPS(1), $ IWORK,ISTAT) CALL STFMAP(IWORK,F_X_MODE,1,NSTEPS(1),IACTS,PWORK,ISTAT) C C Evaluate frequency band C CALL TIMWIDTH(MADRID(PPER),MADRID(PWORK), $ NSTEPS(1),START(1),STEP(1),BANDW,BANDC,NVALUES,WIDTH) C C Return keyword values (NOT IMPLEMENTED) C CALL STKWRD('OUTPUTD',WIDTH(1,1), 1,NVALUES*2,KUN,ISTAT) CALL STKWRD('STARTTSA', BANDC, 1, 1,KUN,ISTAT) C CALL STKWRD('STEPTSA', STEP, 1,1,KUN,ISTAT) C CALL STKWRI('NSTEPS',NSTEPS,1,1,KUN,ISTAT) CALL STTPUT('Kywords START and OUTPUTD are set now.',ISTAT) C C Wind-up C CALL STSEPI C END C C C SUBROUTINE TIMWIDTH(PER,WORK, $ NSTEPS,START,STEP,BANDW,BANDC,NVALUES,WIDTH) C INCLUDE 'MID_REL_INCL:TSA_DEF.INC' DOUBLE PRECISION CRACOW INTEGER NSTEPS,NVALUES DOUBLE PRECISION PER(NSTEPS,2),WORK(NSTEPS) REAL*8 START,STEP,BANDW,BANDC DOUBLE PRECISION WIDTH(NVALUES,2) C INTEGER NPAR1 PARAMETER (NPAR1=4) REAL*8 FINISH DOUBLE PRECISION COND(NPAR1),COVAR(NPAR1,NPAR1) DOUBLE PRECISION VALPEAK,VALLINE,VALCONT,VALHALF,VAL DOUBLE PRECISION PEAKC,LINEC,PEAKW,LINEW,X INTEGER IPEAKC,IPEAKW INTEGER ISTART,NBAND INTEGER I,J,K,ITER,MXITER,IVAL,MODE CHARACTER*80 TEXT INCLUDE 'MID_REL_INCL:TSA_CONST.INC' DATA MXITER/10/ C C Get (trimmed) range of pixels C IF (STEP.EQ.ZERO) THEN CALL STETER(10,'Zero step: check descriptor STEP(1)') ENDIF FINISH=STEP*NSTEPS VAL=BANDC-START ISTART= MIN(NSTEPS,MAX(1,NINT( $ MIN(FINISH,MAX(ZERO,VAL-BANDW))/STEP)))-1 NBAND= MIN(NSTEPS,MAX(1,NINT( $ MIN(FINISH,MAX(ZERO,VAL+BANDW))/STEP))) NBAND=NBAND-ISTART IF (NBAND.LT.10) THEN CALL STETER(11,'Too narrow band') ENDIF C C Compute pixel histogramme and find continuum as median value C DO 1 I=1,NBAND WORK(I)=PER(ISTART+I,1) 1 CONTINUE CALL SORT(NBAND,WORK) VALPEAK=WORK(NBAND) VALCONT=WORK(NBAND/2+1) CALL STTPUT(' --- Line Parameters ---',ISTAT) CALL STTPUT( $ 'Component Level Frequency Width', $ ISTAT) WRITE(TEXT,'(A,1PE10.2,1PE13.5,11X,1PE10.2)') $ 'Continuum ',VALCONT,BANDC,BANDW CALL STTPUT(TEXT,ISTAT) C C Find approximate parameters of the line C VALHALF=(VALPEAK+VALCONT)*HALF IPEAKC=0 DO 2 I=1,NBAND IF (PER(ISTART+I,1).EQ.VALPEAK) THEN IPEAKC=ISTART+I ENDIF IF (WORK(I).LT.VALHALF) THEN IPEAKW=I ENDIF 2 CONTINUE PEAKW=(VALHALF-WORK(IPEAKW))/(WORK(IPEAKW+1)-WORK(IPEAKW)) PEAKW=(NBAND-(IPEAKW*(ONE-PEAKW)+(IPEAKW+1)*PEAKW))*STEP*HALF IF (IPEAKC.EQ.0) THEN CALL STETER(12,'Coding error: condition should not occure') ENDIF PEAKC=START+STEP*IPEAKC WRITE(TEXT,'(A,1PE10.2,1PE13.5,A,1PE7.1,1PE10.2,A,1PE7.1)') $ 'Peak ',VALPEAK,PEAKC,' +/-',STEP,PEAKW,' +/-',STEP CALL STTPUT(TEXT,ISTAT) BANDC=PEAKC C C Fit a line by nonlinear least squares iterations C LINEC=PEAKC LINEW=PEAKW VALLINE=VALPEAK ITER=0 19 CONTINUE ITER=ITER+1 IF (ITER.EQ.1) THEN VAL=LINEC-START ISTART= MIN(NSTEPS,MAX(1,NINT( $ MIN(FINISH,MAX(ZERO,VAL-TWO*LINEW))/STEP)))-1 I = MIN(NSTEPS,MAX(1,NINT( $ MIN(FINISH,MAX(ZERO,VAL+TWO*LINEW))/STEP))) NBAND= I-ISTART IF (NBAND.LT.4) THEN CALL STTPUT('Line very narrow: poor results',ISTAT) ISTART= MIN(NSTEPS,MAX(1,(ISTART+I)/2-2)) NBAND=4 ENDIF ENDIF DO 10 J=1,NPAR1 DO 11 K=1,NPAR1 COVAR(J,K)=ZERO 11 CONTINUE 10 CONTINUE MODE=1 DO 15 I=1,NBAND J=ISTART+I X=START+J*STEP-LINEC CALL PROFILE(LINEW,VALLINE-VALCONT,X,PER(J,1)-VALCONT, $ COND,NPAR1,MODE) DO 16 J=1,NPAR1 DO 17 K=1,NPAR1 COVAR(J,K)=COVAR(J,K)+COND(J)*COND(K) 17 CONTINUE 16 CONTINUE 15 CONTINUE IF (CRACOW(COVAR,NPAR1,NPAR1,NBAND).EQ.ZERO) THEN CALL STETER(15,'Singular conditions for profile') ENDIF C C Take care of initial convergence C VAL= EXP(-(ZERO+ITER-ONE)/ONE) VAL= 0.1*VAL+ONE*(ONE-VAL) LINEC= LINEC + COVAR(1,NPAR1)*VAL LINEW= LINEW *EXP(COVAR(2,NPAR1)*VAL) VALLINE= VALLINE*EXP(COVAR(3,NPAR1)*VAL) VAL=COVAR(1,NPAR1)*COVAR(1,NPAR1)+COVAR(2,NPAR1)*COVAR(2,NPAR1) IF (SQRT(VAL).GT.ONE/NSTEPS*0.1 .AND. ITER.LT.MXITER) GOTO 19 VAL=COVAR(NPAR1,NPAR1)/VALCONT/VALCONT IF (VAL.GT.ONE) THEN WRITE(TEXT,'(A,F12.2,A)') $ '*** Warning: Poor fit (RMS= ',SQRT(VAL),'*Cont.)' CALL STTPUT(TEXT,ISTAT) ENDIF C C Print report C IF (ITER.GE.MXITER) THEN CALL STTPUT('*** Warning: No convergence of fit',ISTAT) ELSE BANDC=LINEC ENDIF WRITE(TEXT,'(A,1PE10.2,1PE13.5,A,1PE7.1,1PE10.2,A,1PE7.1)') $ 'Fit ',VALLINE, $ LINEC,' +/-',SQRT(COVAR(1,1)*COVAR(NPAR1,NPAR1)), $ LINEW,' +/-',SQRT(COVAR(2,2)*COVAR(NPAR1,NPAR1)) CALL STTPUT(TEXT,ISTAT) C C Compute profile width at NVALUES levels C MODE=2 CALL STTPUT(' --- Fitted Line Profile ---',ISTAT) CALL STTPUT(' Width Level Fraction',ISTAT) DO 20 IVAL=1,NVALUES VAL=(IVAL-ONE)/(NVALUES-ONE) VAL=((ONE-VAL)*(VALLINE-VALCONT-VALCONT)+VAL*VALCONT) WIDTH(IVAL,2)=VAL+VALCONT CALL PROFILE(LINEW,VALLINE-VALCONT,WIDTH(IVAL,1),VAL, $ COND,NPAR1,MODE) WRITE(TEXT,'(4(1PE11.2))') WIDTH(IVAL,1),WIDTH(IVAL,2), $ (VALLINE-WIDTH(IVAL,2))/(VALLINE-VALCONT) CALL STTPUT(TEXT,ISTAT) 20 CONTINUE END C C C SUBROUTINE PROFILE(W,H,X,Y,COND,NPAR1,MODE) C INCLUDE 'MID_REL_INCL:TSA_DEF.INC' INTEGER NPAR1,MODE DOUBLE PRECISION COND(NPAR1) DOUBLE PRECISION W,H,X,Y C DOUBLE PRECISION PI,FI,DFIDX,H2,DFDFI,F,WEIGHT INTEGER I INCLUDE 'MID_REL_INCL:TSA_CONST.INC' C PI=TWO*TWO*ATAN(ONE) IF (MODE.EQ.1) THEN DFIDX=PI*HALF/W FI=MAX(-PI,MIN(PI,DFIDX*X)) H2=H*HALF F=(COS(FI)+ONE)*H2 DFDFI=-H2*SIN(FI) COND(1)=-DFDFI*DFIDX COND(2)=DFDFI*(-FI/W)*W COND(3)=F/H*F COND(NPAR1)=Y-F WEIGHT=F/H DO 1 I=1,NPAR1 COND(I)=COND(I)*WEIGHT 1 CONTINUE ELSEIF (MODE.EQ.2) THEN X=TWO*W/PI*ACOS(MIN(ONE,MAX(-ONE,(TWO*Y/H-ONE)))) ELSE CALL STETER(20,'Wrong mode') ENDIF END C C C