C @(#)updat.for	19.1 (ES0-DMD) 02/25/03 13:30:49
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.IDENT:   UPDAT                                                                
C.PURPOSE: Update the mask frame and the buffer setting to a flag value points  
C          with residual greater than kappa*sigma                               
C.INPUT:   MASK  : INPUT_B is the frame masking bright objects in the image     
C          ATEMP : a buffer containing the rebinned input image - used for speed
C          FLMSK : flag = .TRUE. if MASK exists                                 
C          REJ   : value assigned to 'bad' pixels in the buffer ATEMP           
C          NMX,... usual descriptors of the frame MASK                          
C          IDEG  : polinomial degree                                            
C          RKAPPA: parameter for kappa-sigma clipping                           
C          SIGM1 : standard deviation of residuals image-fit.pol. for the       
C                  'good' pixels                                                
C          C     : up to 21 coefficients of the polinomial surface              
C          NSX,... usual descriptors for the temporary frame ATEMP              
C.OUTPUT:  none                                                                 
C ------------------------------------------------------------------            
      SUBROUTINE UPDAT(MASK,ATEMP,FLMSK,REJ,NMX,NMY,RXSTR,RYSTR,RXSTP,          
     +                 RYSTP,NSX,NSY,SXSTR,SYSTR,SXSTP,SYSTP,IDEG,              
     +                 RKAPPA,SIGM1,C)                                          
      IMPLICIT NONE
C
      INTEGER  NMX
      INTEGER  NMY
      REAL     MASK(NMX,NMY)
      REAL     ATEMP(128,128)
      LOGICAL  FLMSK
      REAL     REJ
      DOUBLE PRECISION RXSTR
      DOUBLE PRECISION RYSTR
      DOUBLE PRECISION RXSTP
      DOUBLE PRECISION RYSTP
      INTEGER  NSX
      INTEGER  NSY
      DOUBLE PRECISION SXSTR
      DOUBLE PRECISION SYSTR
      DOUBLE PRECISION SXSTP
      DOUBLE PRECISION SYSTP
      INTEGER  IDEG
      REAL     RKAPPA
      REAL     SIGM1
      DOUBLE PRECISION C(21)
C
      INTEGER  I
      INTEGER  I1MSK
      INTEGER  IX, IY
      INTEGER  J
      INTEGER  J1MSK
      INTEGER  NCV(5)
      REAL     DCX, DCY
      REAL     POLIN
      REAL     RES
      REAL     X1
      REAL     Y1
      REAL     V
      DOUBLE PRECISION CC(6)
      DOUBLE PRECISION XD,YD
      REAL     BTSIG
C
      EXTERNAL POLIN
C                                                                               
      DATA NCV/3,6,10,15,21/                                                    
C                                                                               
C *** calculate residuals and update mask and temporary array                   
C                                                                               
      DCX    = (NSX-1.)/2.                                                      
      DCY    = (NSY-1.)/2.                                                      
C                                                                               
C *** loop : update mask and set -rej- flag in temporary array                  
C                                                                               
      BTSIG  = RKAPPA*SIGM1                                                     
      DO 70 J = 1,NSY                                                           
         IY = J - 1                                                             
         Y1 = SYSTR + SYSTP*FLOAT(IY)                                           
         IF (FLMSK) THEN                                                        
             J1MSK  = (Y1-RYSTR)/RYSTP + 1.5                                    
         END IF                                                                 
                                                                                
         YD     = DCY - IY                                                      
         CALL COEFY(YD,C,CC,IDEG)                                               
C                                                                               
C *** check mask                                                                
C                                                                               
         IF ( .NOT. FLMSK) GO TO 30                           ! y outofbounds   
         IF (J1MSK.LT.1 .OR. J1MSK.GT.NMY) GO TO 30                             
C                                                                               
          DO 20 I = 1,NSX                                                       
             IX     = I - 1                                                     
             XD     = IX - DCX                                                  
             RES    = ATEMP(I,J) - POLIN(XD,CC,IDEG)                            
             X1     = SXSTR + SXSTP*FLOAT(IX)                                   
             I1MSK  = (X1-RXSTR)/RXSTP + 1.5                                    
             IF (I1MSK.LT.1 .OR. I1MSK.GT.NMX) THEN                             
                IF (ATEMP(I,J).EQ.REJ) THEN                                     
                   GO TO 10                                                     
                END IF                                                          
                IF (ABS(RES).GT.BTSIG) THEN                                     
                   ATEMP(I,J) = REJ                                             
                END IF                                                          
             ELSE                                                               
                V      = MASK(I1MSK,J1MSK)                                      
                IF (V.LE.0.) THEN                                               
                   GO TO 10                                                     
                END IF                                                          
                IF (ABS(RES).LE.BTSIG) THEN                                     
                   GO TO 10                                                     
                END IF                                                          
                MASK(I1MSK,J1MSK) = -1                                          
                ATEMP(I,J) = REJ                                                
             END IF                                                             
   10        CONTINUE                                                           
   20     CONTINUE                                                              
          GO TO 60                                                              
C                                                                               
   30     DO 50 I = 1,NSX                                                       
             IF (ATEMP(I,J).EQ.REJ) THEN                                        
                GO TO 40                                                        
             END IF                                                             
             IX     = I - 1                                                     
             XD     = IX - DCX                                                  
             RES    = ATEMP(I,J) - POLIN(XD,CC,IDEG)                            
             IF (ABS(RES).GT.BTSIG) THEN                                        
                ATEMP(I,J) = REJ                                                
             END IF                                                             
   40        CONTINUE                                                           
   50     CONTINUE                                                              
C                                                                               
   60     CONTINUE                                                              
   70 CONTINUE                                                                  
C                                                                               
      RETURN                                                                    
      END                                                                       
