C @(#)irscontsub.for	17.1.1.1 (ES0-DMD) 01/25/02 17:53:05
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
C.COPYRIGHT   (C) 1992 European Southern Observatory
C.IDENT       irscontsub.for
C.AUTHOR      E. Oliva,  Firenze-Arcetri
C.KEYWORDS    Spectroscopy, IRSPEC
C
C.PURPOSE     Fit row by row a poynom and subtract it from a
C             given image (IN_A). Regions to exclude are in dummy
C             image &w (input)
C
C.ALGORITHM   Standard fit of polynoms
C
C
C.INPUT/OUTPUT 'IN_A' : image from which to subtract the continuum
C               &w    : "selection image", fit only pixels &w(x,y)=1.0
C              'OUT_A': continuum subtracted image
C              'OUT_B': continuum image (optional)
C
C
C.VERSION     1.0    Creation     02.09.1992   E. Oliva
C
C-------------------------------------------------------
C
      PROGRAM CONTSUB
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      CHARACTER*64 CUNIT
      CHARACTER*72 IDENT
      CHARACTER*60 FRAMI,FRAMO,FRAMC,FRAMW
      INTEGER NPIX(2),MMPIX(2)
      REAL CUTS(2)
      DOUBLE PRECISION START(2),STEP(2)
c
c pointers
c
      INTEGER*8 IPI,IPO,IPC,IPW
      COMMON /VMR/ MADRID(1)

      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'

      DATA MAXDIM/2/

      IRET=1
      CALL STSPRO('CONTSUB')
C
C GET INPUT FRAME AND MAP IT
C CHECK ALSO WHETHER IT IS A 2D FRAME.....
C
      CALL STKRDC('IN_A',1,1,60,IRET,FRAMI,KUNIT,KNUL,ISTAT)
      CALL CLNFRA(FRAMI,FRAMI,0)
      CALL STIGET(FRAMI,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,MAXDIM,
     ,            NAXIS,NPIX,START,STEP,
     ,            IDENT,CUNIT,IPI,NIN,ISTAT)
C      IF(NAXIS.NE.2) 
C     +     CALL STETER(1,'Input frame must be 2-D')
      NX=NPIX(1)
      NY=NPIX(2)
C
C Get "weight" frame and map it
C
      CALL CLNFRA('&csw',FRAMW,0)
      CALL STIGET(FRAMW,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,MAXDIM,
     ,            NAXIS,NPIX,START,STEP,
     ,            IDENT,CUNIT,IPW,NW,ISTAT)
C
C GET OUTPUT FRAME AND MAP IT
C
      CALL STKRDC('OUT_A',1,1,60,IRET,FRAMO,KUNIT,KNUL,ISTAT)
      CALL CLNFRA(FRAMO,FRAMO,0)
      CALL STIPUT(FRAMO,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     ,            NAXIS,NPIX,START,STEP,
     ,            IDENT,CUNIT,IPO,NOU,ISTAT)
C
C GET FLAG TO check whether a continuum frame is requested
C
      CALL STKRDI('cont',1,1,IRET,ICONT,KUNIT,KNUL,ISTAT)
      IF(ICONT.GT.0) THEN
C
C GET continuum FRAME AND MAP IT
C
        CALL STKRDC('OUT_B',1,1,60,IRET,FRAMC,KUNIT,KNUL,ISTAT)
        CALL CLNFRA(FRAMC,FRAMC,0)
        CALL STIPUT(FRAMC,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     ,              NAXIS,NPIX,START,STEP,
     ,              'continuum frame',CUNIT,IPC,NC,ISTAT)
C
C GET FLAG TO DEFINE degree of interp. polinom
C
      ENDIF
      CALL STKRDI('deg',1,1,IRET,IDEG,KUNIT,KNUL,ISTAT)
C
C CALL ROUTINE WHICH subtracts row by row
C
      IF(ICONT.GT.0) THEN
        CALL SUBC(MADRID(IPI),MADRID(IPO),MADRID(IPW),MADRID(IPC),
     ,            NX,NY,IDEG,ICONT)
      ELSE
        CALL SUBC(MADRID(IPI),MADRID(IPO),MADRID(IPW),DUMMY,
     ,            NX,NY,IDEG,ICONT)
      ENDIF
C
C
C Define and write cuts of output file
C
      CALL MNMX(MADRID(IPO),NPIX,CUTS,MMPIX)
      CALL STDWRR(NOU,'LHCUTS',CUTS,3,2,KUNIT,ISTAT)
C
C RELEASE FILES, UPDATE KEYWORDS AND EXIT
C
      CALL STSEPI
      END
C
C
C
      SUBROUTINE SUBC(A,B,W2D,CC,NX,NY,IDEG,ICONT)
C
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
C
      DIMENSION A(NX,NY),B(NX,NY),CC(NX,NY),W2D(NX,NY)
      REAL*4 DUM
      PARAMETER (MAXPAR=20)
      DOUBLE PRECISION G(MAXPAR,MAXPAR),COEF(MAXPAR)
      DOUBLE PRECISION X(10000),Y(10000)
      DOUBLE PRECISION POL


c
c Define X(I) once forever
c
      DO I=1,NX
        X(I)=I
      ENDDO
C
      NCOEF=IDEG+1
      DO J=1,NY
C
C COPY J-TH row of input IMAGE (A) into  Y 
C same for weight-image (W)
C
       LL = 0
        DO I=1,NX
          IF (W2D(I,J).LE.1.1 .AND. W2D(I,J).GE.0.9) THEN
             NPT = NPT+1
             Y(I)=A(I,J)
C
C determine coefficients of pol. fitting
C
             CALL TDSET2(LL+1,X(I),0.D0,Y(I),IDEG,0,G,COEF,NCOEF,MAXPAR)
             IF (LL.NE.0) THEN
                KMIN = MIN(LL,NCOEF+1)
                DO K1 = 1,KMIN
                 CALL TDHHTR(K1,LL+1,G,COEF,NCOEF,MAXPAR)
                ENDDO
             ENDIF
             LL = MIN(LL+1,NCOEF+1)
          ENDIF
        ENDDO 
        CALL TDSOLV(G,COEF,NCOEF,MAXPAR)
C
C compute output image row as input-pol_fit
C
        DO I=1,NX
          POL = COEF(1)
          DUM = 1.0
          DO K = 2,NCOEF
             DUM = DUM*X(I)
             POL = POL+COEF(K)*DUM
          ENDDO
          VINT=POL
          B(I,J)=A(I,J)-VINT
          IF(ICONT.GT.0) CC(I,J)=VINT
        ENDDO
      ENDDO
C
      RETURN
      END