C @(#)sbtrct.for	19.1 (ES0-DMD) 02/25/03 13:30:48
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 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
      REAL             SXSTR
      REAL             SYSTR
      REAL             SXSTP
      REAL             SYSTP
      DOUBLE PRECISION C(21)                                        
      INTEGER          IDEG
      LOGICAL          FLGSB                                              
      LOGICAL          FLGSU                                              
C
      DOUBLE PRECISION CC(6)                                        
      DOUBLE PRECISION XD,YD                                        
      INTEGER          I, IF, IY, IX
      INTEGER          J
      REAL             DCX, DCY
      REAL             X1
      REAL             XSRV
      REAL             Y1
      REAL             YSRV
C
      REAL             POLIN

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) FRAMED(I) = POLIN(XD,CC,IDEG)                         
            ELSE                                                                
               IF (FLGSU) FRAMED(I) = POLIN(XD,CC,IDEG)                         
            END IF                                                              
   10    CONTINUE                                                               
C                                                                               
         IF     = IF + NIX                                                      
   20 CONTINUE                                                                  
C                                                                               
      RETURN                                                                    
      END                                                                       
