C===========================================================================
C Copyright (C) 1995,2003 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:      Copyright (c) 1987,2003 European Southern Observatory,
C                            all rights reserved
C.IDENTIFICATION: CENTLINE
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 
C.OJO:	to increase the size of the interval, increment all the
C arrays which currently have dimension = 2048 (031112)
C SW1, SW2, W1, W2, WEIGHT
C 
C.VERSION: 840328  JDP  Creation   
C.VERSION: 880703  RHW  ESO-FORTRAN Conversion 
C.VERSION: 900430  RHW  Change to standard graphic interfaces 
C 
C 031112	last modif
C-----------------------------------------------------------------
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   IMNOI, IMNOT, IMNOW
      INTEGER   IMAGE(4)
      INTEGER   NPIX(2),NN, INP
      INTEGER   KEY
      INTEGER   APPFLG
      INTEGER   STATUS,IC(8)
      INTEGER   KUN, KNUL, TID,ACOL, AROW
      INTEGER   ACCESS, PLMODE
      INTEGER   LTYPE, STYPE, BIN

      INTEGER*8 PNTRI, PNTRW, PNTRT

      REAL      SXY,X1,X2,XYS,YOFF
      REAL      XCUR,YCUR,RRBUF(5)
      REAL      RMAGE(4), SVALUE(8), SW1(2048), SW2(2048)

      DOUBLE PRECISION VALUE(8),W1(2048),W2(2048),ACOE(4)
      DOUBLE PRECISION START(2),STEP(2)
      DOUBLE PRECISION DXYS, DSXY

      CHARACTER NAME*60,INFO*60,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' 
      COMMON    /VMR/MADRID(1)
      INCLUDE   'MID_INCLUDE:ST_DAT.INC/NOLIST' 
C
      DATA      NCOLS/6/
      DATA      STYPE/5/
      DATA      BIN/0/
      DATA      YOFF/0.0/
      DATA      ACCESS/1/
      DATA      PLMODE/-1/

      DATA      TABU/' ',' ',' ',' ',' ',' ',' ',' '/
      DATA      TABL/'START','END','CENTER','INT_START','INT_END',
     +               'INT_PEAK','INT_FIT','FWHM'/
      DATA      FORM/'F10.3','F10.3','F10.3','G15.5','G15.5','G15.5',
     +               'G15.5','F10.3'/
      DATA      HEAD1/
     +  '     start         end        center     pixel_value    FWHM'/
      DATA      HEAD2/
     +  '     start         end        center     pixel_value'/
      DATA      PROMPT/' Ready for cursor input ...'/
C
C
      CALL STSPRO('CLINE')                                       !  initialize
      CALL STKRDC('PLCDATA',1,1,60,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
      NVAL    = NPIX(2)
C
C *** allocate virtual memory for a single line of the frame
      CALL PTKRDR('PIXEL',4,IAC,RMAGE)
      IMAGE(1) = RMAGE(1)
      IMAGE(2) = RMAGE(2)
      IMAGE(3) = RMAGE(3)
      IMAGE(4) = RMAGE(4)
      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,
     +               NN,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,
     +               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
      CALL PTOPEN(' ','none',ACCESS,PLMODE)
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 PTGCUR(XCUR,YCUR,KEY,STATUS)
      IF (STATUS.EQ.1) THEN
         GOTO 30
      ELSE IF (KEY.EQ.32) THEN
         CALL TBTCLO(TID,STATUS)
         CALL PTCLOS()
         CALL STSEPI
      ELSE
         LTYPE = 0
         BIN   = 0
         CALL PTDATA(STYPE,LTYPE,BIN,XCUR,YCUR,YOFF,1)
         X1       = XCUR
         VALUE(1) = DBLE(XCUR)
         IPL1     = NINT((X1-XYS)/SXY) + 1
      END IF
C
   40 CONTINUE
      CALL PTGCUR(XCUR,YCUR,KEY,STATUS)
      IF (STATUS.EQ.1) THEN
         GOTO 40
      ELSE IF (KEY.EQ.32) THEN
         CALL TBTCLO(TID,STATUS)
         CALL PTCLOS()
         CALL STSEPI
      ELSE
         LTYPE = 0
         CALL PTDATA(STYPE,LTYPE,BIN,XCUR,YCUR,YOFF,1)
         X2       = XCUR
         VALUE(2) = DBLE(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,
     +          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
         SVALUE(1) = VALUE(1)
         SVALUE(2) = VALUE(2)
         SVALUE(3) = VALUE(3)
         SVALUE(4) = VALUE(4)
         SVALUE(5) = VALUE(5)
         SVALUE(6) = VALUE(6)
         SVALUE(7) = VALUE(7)
         SVALUE(8) = VALUE(8)
C
         RRBUF(1)  = SVALUE(1)
         RRBUF(2)  = SVALUE(2)
         RRBUF(3)  = SVALUE(3)
         RRBUF(4)  = SVALUE(6)
         RRBUF(5)  = SVALUE(8)
         CALL TBRWRR(TID,ISEQ,NCOLS,IC,SVALUE,STATUS)
         IF (IMETH.EQ.-1) THEN
            WRITE (CLINE,9000) SVALUE(1), SVALUE(2),
     +                         SVALUE(3), SVALUE(6), SVALUE(8)
            CALL STKWRR('OUTPUTR',RRBUF,1,5,KUN,STATUS)
         ELSE
            WRITE (CLINE,9000) SVALUE(1), SVALUE(2), 
     +                         SVALUE(3), SVALUE(6)
            CALL STKWRR('OUTPUTR',RRBUF,1,4,KUN,STATUS)
         ENDIF
C
         CALL STTPUT(CLINE,STATUS)
         NCOUNT = NCOUNT + 1
         IF (IMETH.EQ.-1) THEN
            NPTS = I2 - I1 + 1
            LTYPE = 1
            DO INP = 1, NPTS
               SW1(INP) = SNGL(W1(INP))
               SW2(INP) = SNGL(W2(INP))
            ENDDO
            CALL PTDATA(STYPE,LTYPE,BIN,SW1,SW2,YOFF,NPTS)
         END IF
      END IF
C
      GO TO 30
C 
9000  FORMAT(3F13.3,G13.5,G13.5)
      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
      INTEGER          IPL1,IPL2,IMODE,IMETH
      INTEGER          IFAIL,NPX

      REAL             X(*)
      DOUBLE PRECISION XSTR,XSTP
      DOUBLE PRECISION XCENTR,PEAK,CHICHK,XSP,XGO
      DOUBLE PRECISION W1(*),W2(*),ACOE(*),Y1,Y2
C
      CHICHK = 0.005
      XSP    = XSTP
      XGO    = XSTR + (IPL1-1)*XSP
      NPX    = IPL2 - IPL1 + 1
      Y1     = DBLE(X(IPL1))
      Y2     = DBLE(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
C 
      RETURN
      END

C 
      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
      INTEGER  NPIX, IWID, IWERR, I, K

      REAL             WINDOW(*)
      DOUBLE PRECISION AM, AL, AH, DR
      DOUBLE PRECISION 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
      INTEGER           NPIX, IWID, IWERR
      REAL              WINDOW(*)
      DOUBLE PRECISION  AM, 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(*)
      DOUBLE PRECISION XBUF(*),YFIT(*),A(*)          
      DOUBLE PRECISION DELTAA(4),SIGMAA(4),XCNTR,XGO,CHICHK
      DOUBLE PRECISION DXSTP,CHISQR,FLAMDA,OLDCH,DIF,HALF,WIDTH
      DOUBLE PRECISION 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 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                                                                              
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                                                                   
                                                                               
      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, ISTAT                               

      REAL             Y(*)
      DOUBLE PRECISION X(*), SIGMAY(*), DELTAA(*), YFIT(*) 
      DOUBLE PRECISION A(*), DET, FLAMDA, CHISF, FUNCT 
      DOUBLE PRECISION WEIGHT(2048), ALPHA(10,10), BETA(10) 
      DOUBLE PRECISION 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-30 .OR.
     +                ABS(ALPHA(K,K)).LT.1.E-30) THEN
                 CALL STTPUT
     +                ('*** WARNING: Insufficient accuracy: NO RESULT',
     +                ISTAT)
                 CALL STTPUT
     +                ('              Scale your input data first', 
     +                ISTAT)
                 GOTO 987
              ENDIF
              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                                                                      

      DOUBLE PRECISION 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
      INTEGER I
      DOUBLE PRECISION X(*),A(*),Z,XI,Z2

      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

      DOUBLE PRECISION 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)
CC      DOUBLE PRECISION SIGMAY(1),YFIT(1),CHISF
      DOUBLE PRECISION SIGMAY(1),YFIT(1)
      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

      DOUBLE PRECISION X(*),A(*),DELTAA(*),DERIV(*) 
      DOUBLE PRECISION 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                                         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 double. 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
                                         
      DOUBLE PRECISION 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
