C @(#)feet.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
      SUBROUTINE FEET(MASK,ATEMP,IDEG,FLMSK,REJ,NMX,NMY,RXSTR,RYSTR,            
     +                RXSTP,RYSTP,NSX,NSY,SXSTR,SYSTR,SXSTP,SYSTP,ICONT,        
     +                RMEAN,SIGM1,SKEW1,RKURT,C,JER)                            
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++    
C.PURPOSE:   Calculate the coefficients of a 2D polinomial fitting the          
C            sky background,then give the first moments of the distribution     
C            of O-C residuals                                                   
C.ALGORITHM: Lleast squares 2-dim. polinomial fit                               
C.INPUT:     MASK  : frame masking bright objects in the field                  
C            ATEMP : temporary buffer containing the rebinned input image       
C            FLMSK : flag = .TRUE. if MASK exists                               
C            REJ   : value assigned to rejected pixels in frame ATEMP           
C            NMX,... usual descriptors of the frame MASK                        
C            NSX,... descriptors of the frame-buffer ATEMP                      
C            IDEG  : polinomial degree                                          
C.OUTPUT:    ICONT : number of processed 'good' points                          
C            RMEAN : mean value of residuals image-pol.fit for 'good' points    
C            SIGM1 : standard deviation of residuals                            
C            SKEW1 : 3rd moment of the distribution of residuals (skewness)     
C            RKURT : 4th moment of the distribution of residuals (kurtosis)     
C            C     : up to 21 coefficients of the polinomial surface            
C            JER   : error return ; 0=OK   1=interpolation failure              
C --------------------------------------------------------------------          
      IMPLICIT NONE
C
      INTEGER  NMX
      INTEGER  NMY
      REAL     MASK(NMX,NMY)
      REAL     ATEMP(128,128)
      INTEGER  IDEG
      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  ICONT
      REAL     RMEAN
      REAL     SIGM1
      REAL     SKEW1
      REAL     RKURT
      DOUBLE PRECISION C(21)
      INTEGER  JER
C
      INTEGER  I, I1, I2
      INTEGER  I1MSK, I2MSK
      INTEGER  ICX, ICY
      INTEGER  IF1, IF2
      INTEGER  IPERC
      INTEGER  IW1, IW2, IW3, IW4
      INTEGER  IX, IY
      INTEGER  J, J1, J2
      INTEGER  J1MSK, J2MSK
      INTEGER  NC, NC1
      INTEGER  NCV(5)                                                           
      REAL     ANP
      REAL     BUP(128)
      REAL     BLW(128)                                                
      REAL     CONT
      REAL     DCX, DCY
      REAL     D1, D2, D3, D4
      REAL     POLIN
      REAL     RES, RES2
      REAL     X1, X2
      REAL     Y1, Y2
      REAL     V
      REAL     VUP1, VUP2
      REAL     VLW1, VLW2
      DOUBLE PRECISION Z(64,10),A(21,22),B(21,22),CC(6)                   
      DOUBLE PRECISION W1,W2,W3,W4,DMEAN
      DOUBLE PRECISION XD,YD,SIGMA,SKEW,DKURT
C
      EXTERNAL POLIN
      EXTERNAL SHORT,LONG
C
      DATA NCV/3,6,10,15,21/                                                    
C                                                                               
C *** initialize the powers and arrays                                          
C                                                                               
      DO 20 I = 1,64                                                            
         Z(I,1) = I - 0.5                                                       
         DO 10 J = 2,10                                                         
            Z(I,J) = Z(I,J-1)*Z(I,1)                                            
   10    CONTINUE                                                               
   20 CONTINUE                                                                  
C                                                                               
      DO 40 I = 1,21                                                            
         DO 30 J = 1,22                                                         
            A(I,J) = 0.                                                         
   30    CONTINUE                                                               
   40 CONTINUE                                                                  
C                                                                               
      DO 50 I = 1,21                                                            
         C(I)   = 0.D0                                                          
   50 CONTINUE                                                                  
C                                                                               
C ***  initialize variables                                                     
C                                                                               
      ICX    = NSX/2                                                            
      ICY    = NSY/2                                                            
C                                                                               
C *** loop: load image values (testing mask and rejected pixels)                
C                                                                               
      IF1    = 1                                                                
      IF2    = 128* (NSY-1) + 1                                                 
C                                                                               
      DO 90 J1 = 1,ICY                                                          
         J2     = NSY - J1 + 1                                                  
         IY     = ICY - J1 + 1                                                  
         Y1     = SYSTR + SYSTP* (J1-1.)                                        
         Y2     = SYSTR + SYSTP* (J2-1.)                                        
         IF (FLMSK) THEN                        !   nearest pixel in mask frame
            J1MSK  = (Y1-RYSTR)/RYSTP + 1.5                                     
            J2MSK  = (Y2-RYSTR)/RYSTP + 1.5                                     
C                                                                               
C *** mask is extrapolated outside                                              
C                                                                               
            IF (J1MSK.LT.1) THEN                                                
               J1MSK  = 1                                                       
            ELSE IF (J1MSK.GT.NMY) THEN                                         
               J1MSK  = NMY                                                     
            END IF                                                              
                                                                                
            IF (J2MSK.LT.1) THEN                                                
               J2MSK  = 1                                                       
            ELSE IF (J2MSK.GT.NMY) THEN                                         
               J2MSK  = NMY                                                     
            END IF                                                              
         END IF                                                                 
C                                                                               
C *** read needed rows from temporary array                                     
C                                                                               
         CALL ROW(ATEMP,BUP,IF1,NSX)                                            
         CALL ROW(ATEMP,BLW,IF2,NSX)                                            
         IF1    = IF1 + 128                                                     
         IF2    = IF2 - 128                                                     
C                                                                               
         DO 70 IX = 1,ICX                                                       
            I1     = ICX - IX + 1                                               
            I2     = ICX + IX                                                   
            VUP1   = 1.                                                         
            VUP2   = 1.                                                         
            VLW1   = 1.                                                         
            VLW2   = 1.                                                         
C                                                                               
            IF (FLMSK) THEN                                                     
               X1     = SXSTR + SXSTP* (I1-1.)                                  
               X2     = SXSTR + SXSTP* (I2-1.)                                  
               I1MSK  = (X1-RXSTR)/RXSTP + 1.5                                  
               I2MSK  = (X2-RXSTR)/RXSTP + 1.5                                  
               IF (I1MSK.LT.1) THEN                                             
                  I1MSK  = 1                                                    
               ELSE IF (I1MSK.GT.NMX) THEN                                      
                  I1MSK  = NMX                                                  
               END IF                                                           
                                                                                
               IF (I2MSK.LT.1) THEN                                             
                  I2MSK  = 1                                                    
               ELSE IF (I2MSK.GT.NMX) THEN                                      
                  I2MSK  = NMX                                                  
               END IF                                                           
C                                                                               
C *** find out mask values  -  assumed floating p.                              
C                                                                               
               VUP1   = MASK(I1MSK,J1MSK)                                       
               VUP2   = MASK(I2MSK,J1MSK)                                       
               VLW1   = MASK(I1MSK,J2MSK)                                       
               VLW2   = MASK(I2MSK,J2MSK)                                       
            END IF                                                              
C                                                                               
C *** a pixel is 'bad' if rejected via mask frame OR temp. buffer               
C                                                                               
            IF (BUP(I1).EQ.REJ) VUP1   = -1.                                    
            IF (BUP(I2).EQ.REJ) VUP2   = -1.                                    
            IF (BLW(I1).EQ.REJ) VLW1   = -1.                                    
            IF (BLW(I2).EQ.REJ) VLW2   = -1.                                    
C                                                                               
            IF ((VUP1.GT.0.) .AND. (VUP2.GT.0.) .AND.                           
     +         (VLW1.GT.0.) .AND. (VLW2.GT.0)) THEN                             
               D1 = BUP(I1)
               D2 = BUP(I2)
               D3 = BLW(I1)
               D4 = BLW(I2)
               CALL SHORT(D1,D2,D3,D4,IX,IY,Z,A)            
                                                                                
            ELSE                                                                
               IW2    = 1                                                       
               IF (VUP2.LE.0.) THEN                                             
                  BUP(I2) = 0.                                                  
                  IW2    = 0                                                    
               END IF                                                           
                                                                                
               IW1    = 1                                                       
               IF (VUP1.LE.0.) THEN                                             
                  BUP(I1) = 0.                                                  
                  IW1    = 0                                                    
               END IF                                                           
                                                                                
               IW4    = 1                                                       
               IF (VLW2.LE.0.) THEN                                             
                  BLW(I2) = 0.                                                  
                  IW4    = 0                                                    
               END IF                                                           
                                                                                
               IW3    = 1                                                       
               IF (VLW1.LE.0.) THEN                                             
                  BLW(I1) = 0.                                                  
                  IW3    = 0                                                    
               END IF                                                           
                                                                                
               W1     = IW1 + IW2 + IW3 + IW4                                   
               IF (W1.EQ.0.0D0) GO TO 60                                        
               W2     = IW2 + IW4 - IW1 - IW3                                   
               W3     = IW1 + IW2 - IW3 - IW4                                   
               W4     = IW2 + IW3 - IW1 - IW4                                   
               D1 = BUP(I1)
               D2 = BUP(I2)
               D3 = BLW(I1)
               D4 = BLW(I2)
               CALL LONG(W1,W2,W3,W4,D1,D2,D3,D4,IX,IY,Z,A)                                             
            END IF                                                              
   60       CONTINUE                                                            
   70    CONTINUE                                                               
   90 CONTINUE                                                                  
C                                                                               
C *** fill matrix "a" using symmetries                                          
C                                                                               
      CALL FLMTX(A)                                         !   processed points
      ANP    = 4.0D0*A(1,1)                                                     
      DO 110 I = 1,20                                                           
         I1     = I + 1                                                         
         DO 100 J = I1,21                                                       
            A(J,I) = A(I,J)                                                     
  100    CONTINUE                                                               
  110 CONTINUE                                                                  
C                                                                               
C *** calculate coefficients   --->array c                                      
C                                                                               
      NC     = NCV(IDEG)                                                        
      NC1    = NC + 1                                                           
      DO 130 I = 1,NC                                                           
         DO 120 J = 1,NC                                                        
            B(I,J) = A(I,J)                                                     
  120    CONTINUE                                                               
         B(I,NC1) = A(I,22)                                                     
  130 CONTINUE                                                                  
C                                                                               
      CALL SYST(NC,B,C,JER)                               !interpolation failure
      IF (JER.NE.0) RETURN                                                      
C                                                                               
C *** calculate parameters of distribution of residuals                         
C     initialize variables                                                      
C                                                                               
      DCX    = (NSX-1.)/2.                                                      
      DCY    = (NSY-1.)/2.                                                      
      DMEAN  = 0.D0                                                             
      SIGMA  = 0.D0                                                             
      SKEW   = 0.D0                                                             
      DKURT  = 0.D0                                                             
      ICONT  = 0                                                                
C                                                                               
C *** loop : calculate parameters                                               
C                                                                               
      DO 170 J = 1,NSY                                                          
         IY     = J - 1                                                         
         YD     = DCY - IY                                                      
         CALL COEFY(YD,C,CC,IDEG)                                               
         Y1     = SYSTR + SYSTP*FLOAT(IY)                                       
         IF (FLMSK) THEN                                                        
            J1MSK  = (Y1-RYSTR)/RYSTP + 1.5                                     
            IF (J1MSK.LT.1) THEN                                                
               J1MSK  = 1                                                       
            ELSE IF (J1MSK.GT.NMY) THEN                                         
               J1MSK  = NMY                                                     
            END IF                                                              
         END IF                                                                 
C                                                                               
         DO 150 I = 1,NSX                                                       
            IX     = I - 1                                                      
            XD     = IX - DCX                                                   
            V      = 1.                                                         
            IF (FLMSK) THEN                                                     
               X1     = SXSTR + SXSTP*FLOAT(IX)                                 
               I1MSK  = (X1-RXSTR)/RXSTP + 1.5                                  
               IF (I1MSK.LT.1) THEN                                             
                  I1MSK  = 1                                                    
               ELSE IF (I1MSK.GT.NMX) THEN                                      
                  I1MSK  = NMX                                                  
               END IF                                                           
               V      = MASK(I1MSK,J1MSK)                                       
            END IF                                                              
                                                                                
            IF (ATEMP(I,J).EQ.REJ) V      = -1                                  
            IF (V.LE.0) GO TO 140                                               
            ICONT  = ICONT + 1                                                  
            RES    = ATEMP(I,J) - POLIN(XD,CC,IDEG)                             
            RES2   = RES*RES                                                    
            DMEAN  = DMEAN + RES                                                
            SIGMA  = SIGMA + RES2                                               
            SKEW   = SKEW + RES2*RES                                            
            DKURT  = DKURT + RES2*RES2                                          
  140       CONTINUE                                                            
  150    CONTINUE                                                               
  170 CONTINUE                                                                  
C                                                                               
      CONT   = FLOAT(ICONT)                                                     
      DMEAN  = DMEAN/DBLE(CONT)                                                 
      SIGMA  = DSQRT(SIGMA/DBLE(CONT))                                          
      SKEW   = SKEW/ (DBLE(CONT)* (SIGMA**3))                                   
      DKURT  = DKURT/ (DBLE(CONT)* (SIGMA**4))                                  
      RMEAN  = SNGL(DMEAN)                                                      
      SIGM1  = SNGL(SIGMA)                                                      
      SKEW1  = SNGL(SKEW)                                                       
      RKURT  = SNGL(DKURT)                                                      
      IPERC  = (CONT/ (NSX*NSY))*100.                                           
C                                                                               
      RETURN                                                                    
      END                                                                       
