C @(#)spec.for	19.1 (ES0-DMD) 02/25/03 14:29:03
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


