C===========================================================================
C Copyright (C) 1995-2005 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.IDENTIFICATION: BKGRND
C.LANGUAGE:       F77+ESOext
C.AUTHOR:         Capaccioli/Held
C.KEYWORDS:       polinomial fit
C.PURPOSE:        Calculate the coefficients of a 2D polinomial approximation 
C                 of the sky background. Optionally create frames with the 
C                 background subtracted image and the polinomial surface.
C.ALGORITHM:      Coefficients are calculated by an iterative least-squares 
C                 fitting of the pixels in the rebinned image frame. If a 
C                 predefined 'masking frame' is given, only not masked pixels 
C                 are taken in account.
C.INPUT/OUTPUT:   the following keywords are used:
C                 IN_A/C/1/60   name of image data frame
C                 IN_B/C/1/60   name of masking frame
C                 OUT_A/C/1/60  name of output corrected image
C                 OUT_B/C/1/60  name of output background surface
C                 INPUTI/I/1/2  degree, maximum n. of iterations
C                 INPUTR/R/1/1  clipping par. for the 1st iteration
C                 INPUTR/R/2/1  clipping par. for following iter.
C                 INPUTR/R/3/1  upper limit of acceptable skewness
C                 INPUTC/C/1/8  name of the keyword storing coeff.
C.VERSION         840517
C.VERSION:        871124  RHW ESO-FORTRAN Conversion
C 
C 051107	last modif
C 
C ---------------------------------------------------------------------
      PROGRAM BCKGRN

      IMPLICIT NONE
C
      CHARACTER*60 FRAMEA,FRAMEB,FRAMEC,FRAMED
      CHARACTER*8  INPUTC
      CHARACTER    CUNITA*64,IDENTA*72,CUNITB*64,IDENTB*72
      CHARACTER    HISTC*80,HISTD*80
      CHARACTER    MLINE1*40,MLINE2*40,LINE*80

      INTEGER      MADRID,KUN,KNUL,IMFA,IMFB,IMFD,IMFC
      INTEGER      IAV,ISTAT
      INTEGER      JER,MODE,NAXISA,NAXISB,NLOOP
      INTEGER      INPUTI(2)
      INTEGER      NPIXA(3),NPIXB(3),NPIXT(3)
      INTEGER      ICONT,ICONTP
      INTEGER*8      IPNTRA,IPNTRB,IPNTRC,IPNTRD
C
      REAL         AINCR,CONV
      REAL         INPUTR(3)
      REAL         ATEMP(128,128)
      REAL         RKAPPA,REJ,RMEAN,SIGM1,SKEW1,RKURT
      DOUBLE PRECISION STEPA(3),STARTA(3),STEPB(3),STARTB(3)
      DOUBLE PRECISION STEPT(3),STARTT(3)
      DOUBLE PRECISION C(21),DSTORE(7)
C
      LOGICAL      FLMSK,FLGSB,FLGSU,CASEA,CASEB

      INCLUDE      'MID_INCLUDE:ST_DEF.INC/NOLIST'
      COMMON /VMR/MADRID(1)
      INCLUDE      'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
      DATA         MLINE1/'    Iter  Points  Mean (10^6)   Sigma'/
      DATA         MLINE2/'  Skewness  Kurtosis'/
      DATA         LINE/' '/
      DATA         CASEA,CASEB/.FALSE.,.FALSE./
      DATA         REJ/-99999./                      !  minimum relative change
      DATA         CONV/.002/                        !  needed if no mask
      DATA         NPIXB/1,1,1/
C
C *** set up MIDAS environment + enable automatic error abort
      CALL STSPRO('BKGRND')
      CALL STKRDI('MID$MODE',1,1,IAV,MODE,KUN,KNUL,ISTAT)
C
C *** get names of input frame + masking frame
      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUN,KNUL,ISTAT)
      CALL STKRDC('IN_B',1,1,60,IAV,FRAMEB,KUN,KNUL,ISTAT)
C
C *** setup history update
      HISTC = 'background corrected frame obtained from '//FRAMEA
      HISTD = 'background pol. surface using dimensions from '//
     +         FRAMEA
C
C *** get output frames
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEC,KUN,KNUL,ISTAT)
      CALL STKRDC('OUT_B',1,1,60,IAV,FRAMED,KUN,KNUL,ISTAT)
C
C *** read flag keyword INPUTC
      CALL STKRDC('INPUTC',1,1,12,IAV,INPUTC,KUN,KNUL,ISTAT)
C
C *** find out which case
      IF (INPUTC(1:1).EQ.'?') THEN                     ! read useful parameters
         CALL STKRDI('INPUTI',1,2,IAV,INPUTI,KUN,KNUL,ISTAT)
         CALL STKRDR('INPUTR',1,3,IAV,INPUTR,KUN,KNUL,ISTAT)
         CASEA = .TRUE.
      ELSE              ! only create new images reading coeff. from a keyword
         CASEB = .TRUE.
      END IF
C
C *** get input frame and map it
      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     2            2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
     3            IPNTRA,IMFA,ISTAT)
      IF (NAXISA.GT.2) THEN
         CALL STTPUT('*** FATAL: maximum 2 dim. currently supported...',
     +               ISTAT)
         CALL STSEPI
      END IF
C
C *** get masking frame (if exists) and map it
      IF (FRAMEB(1:1).NE.'+') THEN
         FLMSK  = .TRUE.
         CALL STIGET(FRAMEB,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
     2               2,NAXISB,NPIXB,STARTB,STEPB,IDENTB,CUNITB,
     3               IPNTRB,IMFB,ISTAT)
      ELSE
         FLMSK  = .FALSE.
      END IF
C
C *** get output frame for corrected image, map it and set descriptors
      IF (FRAMEC(1:1).NE.'+') THEN
          FLGSB  = .TRUE.
          CALL STIPUT(FRAMEC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     2                NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
     3                IPNTRC,IMFC,ISTAT)
          CALL STDCOP(IMFA,IMFC,4,'HISTORY',ISTAT)
          CALL STDCOP(IMFA,IMFC,4,'LHCUTS',ISTAT)
C         CALL STDWRC(IMFC,'HISTORY',1,HISTC,-1,80,ISTAT)

      ELSE
          FLGSB  = .FALSE.
      END IF
C
C *** get output frame for polinomial surface , map it and set descriptors
      IF (FRAMED(1:1).NE.'+') THEN
          FLGSU  = .TRUE.
          CALL STIPUT(FRAMED,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     2                NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
     3                IPNTRD,IMFD,ISTAT)
          CALL STDCOP(IMFA,IMFD,4,'HISTORY',ISTAT)
          CALL STDCOP(IMFA,IMFD,4,'LHCUTS',ISTAT)
C         CALL STDWRC(IMFD,'HISTORY',1,HISTD,-1,80,ISTAT)

      ELSE
          FLGSU  = .FALSE.
      END IF
C
C *** check parameters values
      IF (CASEA) THEN
          IF (INPUTI(1).LE.0 .OR. INPUTI(1).GT.5) THEN
             CALL STTPUT('*** FATAL: Invalid degree: should be 1 to 5',
     +                   ISTAT)
             CALL STSEPI
          END IF

          IF (INPUTI(2).LE.0) THEN
             CALL STTPUT('*** FATAL: Invalid number of iterations: '//
     +                    'should be greater than zero',ISTAT)
             CALL STSEPI
          END IF

          IF (INPUTI(2).GT.20) THEN
             INPUTI(2) = 20
             CALL STTPUT('*** WARNING: Maximum number of iterations '//
     +                    'set to 20 ',ISTAT)
          END IF

          IF (INPUTR(1).LE.0 .OR. INPUTR(2).LE.0) THEN
             CALL STTPUT('*** FATAL: Invalid kappa: should be '//
     +                   'greater than zero',ISTAT)
             CALL STSEPI
          END IF

          IF (INPUTR(3).LE.0) THEN
             CALL STTPUT('*** FATAL: Invalid skewness: should be '//
     +                   'greater than zero',ISTAT)
            CALL STSEPI
         END IF
      END IF
C
C *** if case b jump to section creating output frames
      IF (CASEB) GO TO 20
C
C *** rebin input image to run faster
      CALL SQZTMP(MADRID(IPNTRA),ATEMP,NPIXA(1),NPIXA(2),STARTA(1),
     +                 STARTA(2),STEPA(1),STEPA(2),NPIXT(1),NPIXT(2),
     +                 STARTT(1),STARTT(2),STEPT(1),STEPT(2))
C
C *** write scaling parameters and degree in keyword INPUTC
      DSTORE(1) = NPIXT(1)
      DSTORE(2) = NPIXT(2)
      DSTORE(3) = STARTT(1)
      DSTORE(4) = STARTT(2)
      DSTORE(5) = STEPT(1)
      DSTORE(6) = STEPT(2)
      DSTORE(7) = INPUTI(1)
      CALL STKWRD('BG_COEF',DSTORE,22,7,KUN,ISTAT)
C
C *** main loop : iterations
      RKAPPA = INPUTR(1)
      DO 10 NLOOP = 1,INPUTI(2)
         IF (NLOOP.NE.1) RKAPPA = INPUTR(2)
C
C *** calculate coefficients for this loop and param. of residuals
         CALL FEET(MADRID(IPNTRB),ATEMP,INPUTI(1),FLMSK,REJ,NPIXB(1),
     +             NPIXB(2),STARTB(1),STARTB(2),STEPB(1),STEPB(2),
     +             NPIXT(1),NPIXT(2),STARTT(1),STARTT(2),STEPT(1),
     +             STEPT(2),ICONT,RMEAN,SIGM1,SKEW1,RKURT,C,JER)
C
C *** check interpolation failure
         IF (JER.NE.0) THEN
            CALL STTPUT('*** FATAL: Interpolation failed',ISTAT)
            CALL STSEPI
         END IF
C
C *** write coefficients in keyword INPUTC
         CALL STKWRD('BG_COEF',C,1,21,KUN,ISTAT)
C
C *** write significant parameters on terminal
         IF (NLOOP.EQ.1) THEN
            CALL STTPUT(MLINE1//MLINE2,ISTAT)
         ENDIF
         WRITE(LINE(1:16),'(2I8)') NLOOP,ICONT
         CALL VFL(RMEAN*1.E6,LINE(19:26))
         CALL VFL(SIGM1,LINE(29:36))
         CALL VFL(SKEW1,LINE(39:46))
         CALL VFL(RKURT,LINE(49:56))
         CALL STTPUT(LINE,ISTAT)
C
C *** if skewness low enough , stop iteration
         IF (ABS(SKEW1).LE.INPUTR(3)) THEN
            CALL STSEPI
         ENDIF
C
C *** also if no significant change occurred
         IF (NLOOP.NE.1) THEN
            AINCR  = FLOAT(ICONTP-ICONT)/ICONTP
            IF (AINCR.LT.CONV) THEN
               CALL STTPUT('*** WARNING: iteration stopped: not '//
     +                     '  enough pixels were rejected',ISTAT)
               GO TO 20
            END IF
         END IF  
         ICONTP = ICONT                         !   save n. of processed points
C
C *** if last iteration , skip mask update
         IF (NLOOP.EQ.INPUTI(2)) GO TO 20
C
C *** if not enough pixels , stop iteration
         IF (ICONT.LT.50) THEN
            CALL STTPUT('*** WARNING: Iteration stopped: '//
     2                  'not enough points ',ISTAT)
             GO TO 20
          END IF
C
C *** now update masking frame and temporary buffer for next iteration
          CALL UPDAT(MADRID(IPNTRB),ATEMP,FLMSK,REJ,NPIXB(1),NPIXB(2),
     +               STARTB(1),STARTB(2),STEPB(1),STEPB(2),NPIXT(1),
     +               NPIXT(2),STARTT(1),STARTT(2),STEPT(1),STEPT(2),
     +               INPUTI(1),RKAPPA,SIGM1,C)
C
   10 CONTINUE                                            !   end main loop
C  
   20 CONTINUE                               !   this is the skip out point
C
C *** read coefficients from keyword (case b)
      IF (CASEB) THEN
         CALL STKRDD('BG_COEF',1,21,C,IAV,KUN,ISTAT)
         CALL STKRDD('BG_COEF',22,7,DSTORE,IAV,KUN,ISTAT)
         NPIXT(1)  = DSTORE(1)
         NPIXT(2)  = DSTORE(2)
         STARTT(1) = DSTORE(3)
         STARTT(2) = DSTORE(4)
         STEPT(1)  = DSTORE(5)
         STEPT(2)  = DSTORE(6)
         INPUTI(1) = DSTORE(7)
      END IF
C
C *** create back.corrected and/or pol. surf. image frame (testing flags)
      IF (.NOT.(FLGSB.OR.FLGSU)) THEN
         CALL STSEPI
      ELSE
         CALL SBTRCT(MADRID(IPNTRA),MADRID(IPNTRC),MADRID(IPNTRD),
     +               NPIXA(1),NPIXA(2),STARTA(1),STARTA(2),STEPA(1),
     +               STEPA(2),NPIXT(1),NPIXT(2),STARTT(1),STARTT(2),
     +               STEPT(1),STEPT(2),C,INPUTI(1),FLGSB,FLGSU)
         CALL STSEPI
      ENDIF
C 
      END


      SUBROUTINE SBTRCT(FRAMEA,FRAMEC,FRAMED,NIX,NIY,DXSTR,DYSTR,               
     +                  DXSTP,DYSTP,NSX,NSY,SXSTR,SYSTR,SXSTP,SYSTP,C,          
     +                  IDEG,FLGSB,FLGSU)                                       
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++:           
C.PURPOSE: Create new image frames with the background subtracted image and/or  
C          the polinomial surface ,testing option flags FLGSB and FLGSU:        
C.INPUT:   FRAMEA : INPUT_A is the input image                                  
C          FLGSB  : flag = .TRUE. if you want backgrond subtraction             
C          FLGSU  : flag = .TRUE. if you want a background surface image        
C.OUTPUT:  FRAMEC : OUTPUT_C is the background subtracted image                 
C          FRAMED : OUTPUT_D is the polinomial surface                          
C   ------------------------------------------------------------------          
      IMPLICIT         NONE
      REAL             FRAMEA(1)                        
      REAL             FRAMEC(1)                        
      REAL             FRAMED(1)                        
      INTEGER          NIX
      INTEGER          NIY
      DOUBLE PRECISION DXSTR
      DOUBLE PRECISION DYSTR
      DOUBLE PRECISION DXSTP
      DOUBLE PRECISION DYSTP
      INTEGER          NSX
      INTEGER          NSY
      DOUBLE PRECISION SXSTR
      DOUBLE PRECISION SYSTR
      DOUBLE PRECISION SXSTP
      DOUBLE PRECISION SYSTP
      DOUBLE PRECISION C(21)                                        
      INTEGER          IDEG
      LOGICAL          FLGSB                                        
      LOGICAL          FLGSU                                              
C
      INTEGER          I, IX, IY, IF
      INTEGER          J
      REAL             DCX, DCY
      REAL             POLIN
      REAL             Y1, YSRV
      REAL             X1, XSRV
      DOUBLE PRECISION CC(6)
      DOUBLE PRECISION XD,YD                                        
C                                                                               
C *** initialize variables                                                      
C                                                                               
      DCX    = (NSX-1.)/2.                                                      
      DCY    = (NSY-1.)/2.                                                      
C                                                                               
C *** fill image(S)                                                             
C                                                                               
      IF     = 1                                                                
      DO 20 J = 1,NIY                                                           
         IY     = J - 1                                                         
         Y1     = DYSTR + DYSTP*FLOAT(IY)                                       
C                                                                               
C *** find pixel coord. in temporary array used by fitting routine              
C                                                                               
         YSRV   = (Y1-SYSTR)/SYSTP                                              
         YD     = DCY - YSRV                                                    
C                                                                               
C *** initialize calculation of coefficients                                    
C                                                                               
         CALL COEFY(YD,C,CC,IDEG)                                               
C                                                                               
         DO 10 I = IF,IF + NIX - 1                                              
            IX     = I - IF                                                     
            X1     = DXSTR + DXSTP*FLOAT(IX)                                    
            XSRV   = (X1-SXSTR)/SXSTP                                           
            XD     = XSRV - DCX                                                 
            IF (FLGSB) THEN                                                     
               FRAMEC(I) = FRAMEA(I) - POLIN(XD,CC,IDEG)                        
               IF (FLGSU) THEN
                  FRAMED(I) = POLIN(XD,CC,IDEG)                         
               ENDIF
            ELSE                                                                
               IF (FLGSU) THEN
                  FRAMED(I) = POLIN(XD,CC,IDEG)                         
               ENDIF
            END IF                                                              
   10    CONTINUE                                                               
C                                                                               
         IF     = IF + NIX                                                      
   20 CONTINUE                                                                  
C                                                                               
      RETURN                                                                    
      END                                                                       
