C @(#)tsawdt.for	19.1 (ESO-DMD) 02/25/03 13:33:26
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
