C @(#)rfotaddst.for	19.1 (ES0-DMD) 02/25/03 13:30:13
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
       PROGRAM ADDST
C+++
C.IDENTIFICATION: RFOTADDST
C.PURPOSE: Create an artifial image identical to the original with inserted
C          subframes
C.AUTHOR:  R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola
C          Osservatorio Astronomico di Roma
C.VERSION: 880922 RXB First version in MIDAS
C          890217 RHW Converted to ST intercafes, standard code
C          890808 RHW MIDAS table file system implemented 
C.VERSION  900122 RHW IMPLICIT NONE added; all variables defiend	
C----
      IMPLICIT     NONE
      INCLUDE     'MID_REL_INCL:RFOTDECL.INC'
C
      INTEGER      MXDIM
      INTEGER      MXCMP
      INTEGER      MXWND
      PARAMETER    (MXDIM=1024)
      PARAMETER    (MXCMP=500)
      PARAMETER    (MXWND=100)
C
      INTEGER      ICOL(12)
      INTEGER      IDUM2, IDUM3
      INTEGER      IDUM6, IDUM7, IDUM8, IDUM9, IDUM14, IDUM15
      INTEGER      IXA(MXCMP), IYA(MXCMP)
      INTEGER      ISTAT, IAV
      INTEGER      I, IY, IDX, IDY, IGRP, IRN, IROW
      INTEGER      IBO, IBA
      INTEGER      IMFO
      INTEGER*8    IPOU
      INTEGER      IX1, IY1
      INTEGER      IDIFX, IDIFY, INDX, INDY
      INTEGER      JX, JY
      INTEGER      J, JBO
      INTEGER      KVI(3)
      INTEGER      KUN, KNUL
      INTEGER      K, KW
      INTEGER      LX0, LY0, LX, LY, IX0, IY0
      INTEGER      MADRID(1)
      INTEGER      NPIX(3)
      INTEGER      NAXIS, NPL, NL
      INTEGER      NCOMP, NCO, NRO, NSC, NSTOT
      INTEGER      NSA, ICROW, NUMST
      INTEGER      TIDCAT, TIDREG
      INTEGER      TINULL
C
      DOUBLE PRECISION BEGIN(3),STEP(3)
      DOUBLE PRECISION TDNULL,TDTRUE,TDFALS
      REAL         PLX, PLY
      REAL         PPP, SAV
      REAL         RVV(MXWND,MXWND)
      REAL         RXV(MXWND)
      REAL         RMV(MXDIM)
      REAL         RDUM4, RDUM5, RDUM6, RDUM7, RDUM8
      REAL         RDUM9, RDUM10, RDUM11, RDUM12, RDUM13
      REAL         TAB(12)
      REAL         TRNULL
      REAL         TBLSEL
      REAL         VIV(MXWND,MXWND)
      REAL         XXX, YYY
      REAL         XQ(MXCMP)
      REAL         YQ(MXCMP)
      REAL         RAN3
C
      CHARACTER    IDENT*72,CUNIT*80
      CHARACTER    FRAMO*60,REGFIL*60,CATFIL*60
C     CHARACTER    VET*8
      CHARACTER    STRING*80
      LOGICAL      SFLAG
      LOGICAL      NUL(12)
C
      INCLUDE     'MID_INCLUDE:ST_DEF.INC'
      COMMON      /VMR/MADRID
      INCLUDE     'MID_INCLUDE:ST_DAT.INC'
      DATA         ICOL/2,3,4,5,6,7,8,9,10,11,12,13/
C
C *** Hello, is MIDAS out there?
      CALL STSPRO('ADDST')
C
C *** get the frame and read all relevant descriptors
      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMO,KUN,KNUL,ISTAT)
      CALL STIGET(FRAMO,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
     2            3,NAXIS,NPIX,BEGIN,STEP,IDENT,CUNIT,IPOU,IMFO,ISTAT)
      NPL = NPIX(1)
      NL  = NPIX(2)
C
C *** Input table
      CALL STKRDC('IN_B',1,1,60,IAV,REGFIL,KUN,KNUL,ISTAT)
      CALL TBTOPN(REGFIL,0,TIDREG,ISTAT)
C
C *** Output cat_file
      CALL STKRDC('OUT_B',1,1,60,IAV,CATFIL,KUN,KNUL,ISTAT)
      CALL CATINI(CATFIL,TIDCAT)
C
C *** number of position and the subarray dimensions
      CALL STKRDI('INPUTI',1,2,IAV,KVI,KUN,KNUL,ISTAT)
      IDX   = KVI(1)
      IDY   = KVI(2)
      PLX   = NPL-IDX-2
      PLY   = NL -IDY-2
      CALL STKRDI('INPUTI',3,1,IAV,NCOMP,KUN,KNUL,ISTAT)
C
C *** get the registration table opened
      CALL TBMNUL(TINULL,TRNULL,TDNULL)
      CALL TBMCON(TBLSEL,TDTRUE,TDFALS)
      CALL TBIGET(TIDREG,NCO,NRO,NSC,KW,NSA,ISTAT)
      IF (ISTAT.NE.0) THEN
         STRING = '*** FATAL: Problems with opening '//
     2            'the registration file'
         CALL STTPUT(STRING,ISTAT)
         CALL STSEPI
      ENDIF
C
      ICROW = 0
      NUMST = 0
      IRN   = 0
C
      DO 100 IROW = 1,NRO
         CALL TBSGET(TIDREG,IROW,SFLAG,ISTAT)
         IF (SFLAG) THEN
            NUMST = NUMST + 1
            CALL TBRRDR(TIDREG,IROW,12,ICOL,TAB,NUL,ISTAT)
            IX0 = TAB(1)-IDX/2-BEGIN(1)+1
            IY0 = TAB(2)-IDY/2-BEGIN(2)+1
            LX  = IDX
            LY  = IDY
            JX  = LX/2
            JY  = LY/2
            IF (IX0.LT.1) THEN
               LX  = LX+IX0-1
               IX0 = 1
               JX  = TAB(1)-BEGIN(1)
            END IF
C
            IF (IX0+LX-1.GT.NPL) THEN
               LX  = NPL-IX0+1
            ENDIF
            IF (IY0.LT.1) THEN
               LY  = LY+IY0-1
               IY0 = 1
               JY  = TAB(2)-BEGIN(2)
            END IF
            IF (IY0+LY-1.GT.NL) THEN
               LY  = NL-IY0+1
            ENDIF
C
            DO 110 IBO = 1,LY
               IY = IY0+IBO-1 
               CALL REALIN(NPL,NL,IY,IX0,LX,MADRID(IPOU),RXV)
               DO 111 JBO = 1,LX
                  RVV(JBO,IBO) = RXV(JBO)-TAB(4)           ! subtract background
  111          CONTINUE
  110       CONTINUE
C
            LX0 = LX
            LY0 = LY
C
            DO 120 IBA=1,NCOMP 
               IRN      = IRN + IX0 + IY0
               XXX      = RAN3(IRN)
               YYY      = RAN3(IRN)
               XQ(IBA)  = XXX*PLX+IDX/2+1
               YQ(IBA)  = YYY*PLY+IDY/2+1
               IXA(IBA) = XQ(IBA)                          ! get randow position
               IYA(IBA) = YQ(IBA)
C
               ICROW = ICROW + 1
               IGRP  = ICROW  
C
               IX1   = IXA(IBA)-JX
               IY1   = IYA(IBA)-JY
               IF (IX1.LT.1) THEN
                  LX  = LX+IX1-1
                  IX1 = 1
               END IF
C
               IF (IX1+LX-1.GT.NPL) THEN
                  LX  = NPL-IX1+1
               ENDIF
               IF (IY1.LT.1) THEN
                  LY  = LY+IY1-1
                  IY1 = 1
               END IF
               IF (IY1+LY-1.GT.NL) THEN
                  LY  = NL-IY1+1
               ENDIF
               XQ(IBA) = XQ(IBA)+BEGIN(1)-1                      ! position in x
               YQ(IBA) = YQ(IBA)+BEGIN(2)-1                      ! position in y
               CALL CATTWR(TIDCAT,ICROW,IGRP,XQ(IBA),YQ(IBA),TAB(3),
     2                            TAB(8),RDUM5,RDUM6,RDUM7,RDUM8,RDUM9,
     3                            RDUM10,RDUM11,RDUM12,RDUM13)
C
C *** determine the value of the array VIV(i,j) where i=1 corresponds
C     with IX1 and j=1 corresponds with IY1
               IDIFX    = LX-LX0
               IDIFY    = LY-LY0
               DO 121 J = 1,LY
                  DO 122 K = 1,LX
                     VIV(K,J) = RVV(K-IDIFX,J-IDIFY)         ! input window data
  122             CONTINUE
  121          CONTINUE
C
               DO 200 J = IY1,IY1+LY-1
                  INDY = J-IY1+1
                  CALL REALIN(NPL,NL,J,1,NPL,MADRID(IPOU),RMV)    !output window
                  DO 210 I = IX1, IX1+LX-1
                     INDX = I-IX1+1
                     PPP = RMV(I) + VIV(INDX,INDY)               ! add component
C                    PPP = VIV(INDX,INDY)                
                     IF (PPP.LT.0) THEN
                        PPP = 0.
                     ENDIF
                     RMV(I) = PPP
                     IF (RMV(I).GT.0.) SAV    = RMV(I)
                     IF (RMV(I).LE.0.) RMV(I) = SAV
  210             CONTINUE
                  CALL WRILIN(NPL,NL,J,1,NPL,MADRID(IPOU),RMV)
  200          CONTINUE
  120       CONTINUE
         ENDIF
  100 CONTINUE
C
C *** write the number of objects in the catalogue descriptor
      NSTOT = NUMST*NCOMP
      CALL CATDWR(TIDCAT,NSTOT,IDUM2,IDUM3,RDUM4,RDUM5,IDUM6,IDUM7,
     2                   IDUM8,IDUM9,RDUM10,RDUM11,RDUM12,RDUM13,
     3                   IDUM14,IDUM15)
      CALL TBSINI(TIDCAT,ISTAT)
      CALL TBTCLO(TIDCAT,ISTAT)
      CALL TBTCLO(TIDREG,ISTAT)
      CALL STSEPI
      END
