C @(#)astrm.for	19.1 (ES0-DMD) 02/25/03 13:30:45
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: ASTRM                                                         
C.PURPOSE:        Routine for main program star written by P.B.Stetson          
C                 performs the gaussian astrometric fitting                     
C.LANGUAGE:       F77+ESOext                                                    
C.AUTHOR:         P. Stetsan                                                    
C.ALGORITHM:      see Stetson,  P.B.,  1979  Astron. J. 84  1149                
C.INPUT/OUTPUT:   ->  KRX    input array                                        
C                 ->  JX     low  side local minimun                            
C                 ->  LX     high side local minimum                            
C                 <-> XI     input/ouput estimates of the image centre          
C                 <-  S1     error estimate of the output image centre          
C                            0. fitting procedure failed or if                  
C                            KOUNT .NE. 0  then KOUNT too small                 
C                 <-> SX     input/output estimate of the image width           
C                 <-  DO     output estimate of the image height                
C                 <-  B      output estimate of the background marginal iensity.
C                 <-> KOUNT  number of iterations required/ status              
C                            0 fitting procedure failed                         
C.COMMENTS:       If the solution fails to converge, XI will come back with     
C                 the input value.... with S1 = 0. just in case you want to try 
C                 doing the photometry with the estimated centre from SERCH.    
C.KEYWORDS :      photometry,astrometry                                         
C.VERSION:        8212?? P. Stetsan                                             
C.VERSION:        871123 R.H. Warmels  ESO Fortran Conversion                   
C.VERSION:        ?????? H. Waldhausen Vax implementation                       
C ---------------------------------------------------------------------         
      SUBROUTINE ASTRM(KRX,JX,LX,XI,S1,SX,DO,B,KOUNT)                           
C                                                                               
      IMPLICIT         NONE
      DOUBLE PRECISION KRX(1)
      INTEGER          JX
      INTEGER          LX
      REAL             XI
      REAL             S1
      REAL             SX
      REAL             DO
      REAL             B
      INTEGER          KOUNT
C
      INTEGER          I
      INTEGER          J
      INTEGER          KNT
      INTEGER          L
      REAL             D, DX, DELX
      REAL             EX
      REAL             F, FF, FOLD
      REAL             S, S2, SXSI
      REAL             XO
      REAL             Y
      DOUBLE PRECISION CTF(5)                         
      DOUBLE PRECISION CTC(5,5)                         
      DOUBLE PRECISION V(5)                         
      DOUBLE PRECISION DABS                         
      REAL             C(3)                                                     
C                                                                               
C                                                                               
C..  this part of the program uses the estimates of the image centre            
C..  and width from SERCH to derive initial estimates of the image              
C..  height and background density by linear least-squares.                     
C                                                                               
      IF ((LX-JX).GT.4) THEN                                                    
         XO     = XI                                                            
         DELX   = LX - JX + 1                                                   
         XO     = XI                                                            
         DELX   = LX - JX + 1                                                   
         S1     = 0.0                                                           
         S2     = 0.0                                                           
         DO     = 0.0                                                           
         B      = 0.0                                                           
         SXSI   = 0.5/(SX*SX)                                                   
C                                                                               
         DO 10 I = JX,LX                                                        
            DX     = I - XO                                                     
            EX     = EXP(-DX*DX*SXSI)                                           
            DO     = DO + EX*KRX(I)                                             
            B      = B + EX*EX                                                  
            S1     = S1 + KRX(I)                                                
            S2     = S2 + EX                                                    
   10    CONTINUE                                                               
C                                                                               
         F      = LX - JX + 1                                                   
         DO     = (DO-S1*S2/F)/ (B-S2*S2/F)                                     
         B      = (S1-DO*S2)/F                                                  
         IF (DO.LT.20.0) D  = 20.0                                              
         IF (SX.LT.2.0)  SX = 2.0                                               
         FOLD   = 1.0E38                                                        
C                                                                               
C..  main iteration loop                                                        
         DO 90 KNT = 1,KOUNT                           !   iteration loop       
            S1     = 0.0                                                        
            S2     = 0.0                                                        
            FF     = 0.0                                                        
            SXSI   = 1.0/ (SX*SX)                                               
            DO 30 I = 1,4                                                       
               CTF(I) = 0.0D0                                                   
               DO 20 J = I,4                                                    
                  CTC(I,J) = 0.0D0                                              
   20          CONTINUE                                                         
   30       CONTINUE                                                            
C                                                                               
C..  set normal equations                                                       
            DO 60 L = JX,LX                              !   equation loop      
               DX     = L - XO                                                  
               EX     = 0.5*DX*DX*SXSI                                          
               C(2)   = EXP(-EX)                                                
               S      = DO*C(2)                                                 
               Y      = S + B                                                   
               F      = Y - KRX(L)                                              
               C(1)   = DX*SXSI*S                                               
               EX     = C(1)*C(1)                                               
               S1     = S1 + EX*F*F                                             
               S2     = S2 + EX                                                 
               FF     = FF + F*F                                                
               C(3)   = C(1)*DX/SX                                              
               DO 50 I = 1,3                                                    
                  CTF(I) = CTF(I) + F*C(I)                                      
                  CTC(I,4) = CTC(I,4) + C(I)                                    
                  DO 40 J = I,3                                                 
                     CTC(I,J) = CTC(I,J) + C(I)*C(J)                            
   40             CONTINUE                                                      
   50          CONTINUE                                                         
               CTF(4) = CTF(4) + F  !   equation loop                           
   60       CONTINUE                                                            
C                                                                               
            CTC(4,4) = FLOAT(LX-JX+1)                                           
            DO 80 I = 2,4                                                       
               L      = I - 1                                                   
               DO 70 J = 1,L                                                    
                  CTC(I,J) = CTC(J,I)                                           
   70          CONTINUE                                                         
   80       CONTINUE                                                            
C                                                                               
C..  solve normal equations                                                     
            CALL INVRS(CTC,4)                                                   
            CALL VMUL(CTC,4,CTF,V)                                              
C                                                                               
C..  if DA is the correction to the estimate of A, then replace DA by           
C..  DA/(1+(DA/X)) = DA' prevents an ill-determined solution from               
C..  jumping around too much in subsequent iterations. DA' is less              
C..  than X, and DA' goes to DA for DA much less than X.                        
C..  If X = A/R for some R greater than one, then an intrinsically              
C..  positive parameter will never be permitted to pass through zero.           
C..  This can be helpful for parameters like image width or height.             
C                                                                               
            XO     = XO - V(1)/ (1.0+DABS(V(1)/5.0))                            
            DO     = DO - V(2)/ (1.0+2.0*DABS(V(2)/DO))                         
            SX     = SX - V(3)/ (1.0+2.0*DABS(V(3)/SX))                         
            B      = B - V(4)/ (1.0+2.0*DABS(V(4)/B))                           
C                                                                               
C..  failure                                                                    
            IF ((DO.LT.20.0) .OR. (SX.LT.1.0) .OR.                              
     +         (SX.GT.0.5*DELX)) THEN                                           
               KOUNT  = 0                                                       
               S1     = 0.0                                                     
               RETURN                                                           
            END IF                                                              
C                                                                               
C..  success                                                                    
            IF ((DABS(V(1)).LT.0.01) .AND.                                      
     +         ((FOLD-FF)/FOLD.LT.0.05)) THEN                                   
               XI     = XO                                                      
               S1     = S1*CTC(1,1)*DELX/ (S2* (DELX-4.0))                      
               RETURN                                                           
            END IF                                                              
C                                                                               
            FOLD   = FF                                                         
   90    CONTINUE                                         !   iteration loop    
C                                                                               
C..  number of iterations too small                                             
         S1     = 0.0                                                           
C                                                                               
C..  undersampled                                                               
      ELSE                                                                      
         KOUNT  = 0                                                             
         S1     = 0.0                                                           
      END IF                                                                    
C..  exit                                                                       
      RETURN                                                                    
      END                                                                       
