C @(#)rfotsky.for	19.1 (ES0-DMD) 02/25/03 13:30:15
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 ROMSKY
C++++
C.PURPOSE: Register the sky background counts in rectangular or
C          annular regions. The sky is defined as the mode of the
C          pixel intensity distribution in the region.
C.AUTHORS: R. Buonanno, G. Buscema, C. Corsi, I. Ferraro, G. Iannicola
C          Osservatorio Astronomico di Roma
C.VERSION: 16.01.87 First running version outside MIDAS         R. Buonanno
C          17.12.87 Installation in MIDAS                       R.H. Warmels
C          --.11.88 new version                                 R. Buonanno
C          16.12.88 rewritten for portability and readability   R.H. Warmels
C          02.01.89 sky background written in MIDAS tables      R.H. Warmels
C.VERSION  910123 RHW IMPLICIT NONE included; all variables defined
C---
       IMPLICIT     NONE
C
       INTEGER      NT
       INTEGER      NZ
       INTEGER      NDATA
       INTEGER      NBIN
       INTEGER      NTZ
       PARAMETER    (NT=256)
       PARAMETER    (NZ=1024)
       PARAMETER    (NBIN=1024)
       PARAMETER    (NDATA=2048)
C      PARAMETER    (NDATA=100000)
       PARAMETER    (NTZ=NT*NZ)
C
       INTEGER      MADRID(1)
       INTEGER      EC, ED, EL
       INTEGER      I, J, IY, K
       INTEGER      MIK, KMOD, LMOD, IOLD, II
       INTEGER      IAC, KUN, KNUL, ISTAT
       INTEGER      NAXIS, NPL, NL, NC
       INTEGER      IXW, IYW
       INTEGER      IPX, IPY
       INTEGER      IX0, IY0, IX1, IY1
       INTEGER      LX, LY, LXY
       INTEGER      NRX, NRY
       INTEGER      NDUM
       INTEGER      LV
       INTEGER      MI(NT,NZ)
       INTEGER      IVI(NBIN)
       INTEGER      NPIX(2)
       INTEGER      AREA(6)
       INTEGER      ICOL(256)
       INTEGER      MAK
       INTEGER      TIDSKY
       INTEGER*8    IPNTR
       INTEGER      IMF
C
       REAL         AMI, AMA
       REAL         BACKGR(NT)
       REAL         BINS, BINL, BINH, BINSS
       REAL         DATA(NDATA)
       REAL         DENS(4)
       REAL         FO(NTZ)
       REAL         PSX, PSY
       REAL         RBIN(2)
       DOUBLE PRECISION BEGIN(2),STEP(2)
C
       CHARACTER    IDENT*72
       CHARACTER    CUNIT*72
       CHARACTER    FRAME*60
       CHARACTER    STRING*80
       CHARACTER    SKYFIL*60
       CHARACTER*16 TUNIT,TLABEL,TTYPE,TFORM
C      
       INCLUDE      'MID_INCLUDE:ST_DEF.INC/NOLIST'
       COMMON       /VMR/MADRID
       INCLUDE      'MID_INCLUDE:ST_DAT.INC/NOLIST'
C
       DATA         TUNIT/' '/
       DATA         TLABEL/'BACKGR_LEVEL'/
       DATA         TTYPE/'R*4'/
       DATA         TFORM/'G13.6'/
C
  901  FORMAT(10X,'Minimum and maximum used:', G13.6,', ', G13.6)
  902  FORMAT(10X,'Area in pixel units: [', 
     2            I5, ',', I5, ':', I5, ',', I5, ']')
C 903  FORMAT('*** FATAL: Too many histogram bins; maximum is ',I5)
  904  FORMAT(10X,'Number of cells in x: ',I3, '; in y:', I3)
C
       CALL STSPRO('RFOTSKY')
C
       CALL STKRDC('IN_A',1,1,60,IAC,FRAME,KUN,KNUL,ISTAT)
       CALL STIGET(FRAME,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     2             2,NAXIS,NPIX,BEGIN,STEP,IDENT,CUNIT,
     3             IPNTR,IMF,ISTAT)
       NPL = NPIX(1)
       NL  = NPIX(2)
       IXW = BEGIN(1)
       IYW = BEGIN(2)
C
       CALL STKRDC('OUT_A',1,1,60,IAC,SKYFIL,KUN,KNUL,ISTAT)
C
       CALL STKRDC('INPUTC',1,1,80,IAC,STRING,KUN,KNUL,ISTAT)
       CALL EXTCOO(IMF,STRING,NAXIS,NDUM,AREA(1),AREA(3),ISTAT)
       IX0 = AREA(1)
       IY0 = AREA(2)
       IX1 = AREA(3)
       IY1 = AREA(4)
       LX  = AREA(3) - AREA(1) + 1
       LY  = AREA(4) - AREA(2) + 1
       LXY = LX*LY
C
       CALL STKRDI('INPUTI',1,2,IAC,AREA(5),KUN,KNUL,ISTAT)   !# cells in #x,#y
       NRX = AREA(5)
       NRY = AREA(6)  
       NC  = NRX*NRY                               !total number of rectanguals
C
C *** INFO to the use
       STRING = '*** INFO: Sky values will be calculated in '//
     2          'rectangular grid'
       CALL STTPUT(STRING,ISTAT)
C
C *** get binning parameters for the histogram determination
       CALL STKRDR('INPUTR',1,2,IAC,RBIN,KUN,KNUL,ISTAT)                !# bins
       BINL = RBIN(1)
       BINH = RBIN(2)
       IF (ABS(BINH-BINL).LT.10.0E-10) THEN
C
C *** set min and max pixel values
          CALL STDRDR(IMF,'LHCUTS',1,4,IAC,DENS,KUN,KNUL,ISTAT)
          IF (DENS(1).LT.DENS(2)) THEN
             BINL = DENS(1)
             BINH = DENS(2)
          ELSE
             BINL = DENS(3)
             BINH = DENS(4)
          ENDIF
C
          IF ((BINH-BINL).LT.10.0E-10) THEN
             CALL MINMAX(MADRID(IPNTR),LXY,BINL,BINH)
             IF ((BINH-BINL).LT.10.0E-10) THEN
                STRING = '*** FATAL: frame minimum = frame maximum ...'
                CALL STTPUT(STRING,ISTAT)
                CALL STSEPI
             ENDIF
          ENDIF
       ENDIF
       BINS = (BINH-BINL)/NBIN
C
       WRITE (STRING,901) BINL, BINH
       CALL STTPUT(STRING,ISTAT)
C
C ***  Create output table
       CALL STECNT('GET',EC,ED,EL)
       CALL STECNT('PUT',1,0,0)
       CALL TBTINI(SKYFIL, 0, F_O_MODE, NRX, NRY, TIDSKY, ISTAT)
       IF (ISTAT.NE.0) THEN
          STRING = '*** FATAL: Problems with sky file; try again ...'
          CALL STTPUT(STRING,ISTAT)
          CALL STSEPI
       ENDIF
       CALL STECNT('PUT',EC,ED,EL)
       DO 10 I = 1, NRX
          CALL TBCINI(TIDSKY, D_R4_FORMAT, 1, TFORM, TUNIT,
     &                TLABEL, ICOL(I), ISTAT)
   10  CONTINUE
C
       WRITE (STRING,902) IX0,IY0,IX1,IY1
       CALL STTPUT(STRING,ISTAT)
       WRITE(STRING,904) NRX,NRY
       CALL STTPUT(STRING,ISTAT)
C
       DO 20 I=1,NT
          DO 21 J=1,NZ
             MI(I,J) = 0
   21     CONTINUE
   20  CONTINUE
        DO 30 I=1,NBIN
          IVI(I) = 0
   30  CONTINUE
C
C *** determine the intensity histogram
       DO 40 IY=IY0,IY1
          CALL REALIN(NPL,NL,IY,IX0,LX,MADRID(IPNTR),DATA)
          MAK = -1000000
          MIK = -MAK
          DO 41 J = 1,LX
             IF (DATA(J).GT.BINL .AND. DATA(J).LT.BINH) THEN
                K      = (DATA(J)-BINL)/BINS+1
                IVI(K) = IVI(K)+1
                MAK    = MAX0(MAK,K)
                MIK    = MIN0(MIK,K)
             END IF
   41     CONTINUE
   40  CONTINUE
C
C ** get the maximum in the histogram
       KMOD = 0
       LMOD = 0
       DO 50 I = MIK,MAK
          IF (IVI(I).GT.LMOD) THEN
             LMOD = IVI(I)
             KMOD = I
          END IF
   50  CONTINUE
C
       BINSS = BINS/10.0
       AMI = BINL+KMOD*BINS-NZ*BINSS/2
       IF (AMI.LT.BINL) THEN
          AMI=BINL
       ENDIF
       AMA = AMI+NZ*BINSS
C
       PSX = FLOAT(LX)/NRX
       PSY = FLOAT(LY)/NRY
       DO 60 IY=IY0,IY1
          IOLD = IY-IY0+1
          CALL REALIN(NPL,NL,IY,IX0,LX,MADRID(IPNTR),DATA)
          DO 61 J = 1,LX
             IF (DATA(J).GE.AMI .AND. DATA(J).LT.AMA) THEN
                IF (PSX-INT(PSX).EQ.0.) THEN
                   IPX=(J-1)/INT(PSX)+1
                ELSE
                   IPX=J/PSX+1
                END IF
                IF (PSY-INT(PSY).EQ.0.) THEN
                   IPY=(IY-IY0)/INT(PSY)+1
                ELSE
                   IPY = IOLD/PSY+1
                END IF
                K        = IPX+(IPY-1)*NRX
                II       = (DATA(J)-AMI)/BINS + 1
                MI(K,II) = MI(K,II)+1
             END IF
   61     CONTINUE
   60  CONTINUE
C
       DO 70 I = 1,NC                                          !  to trova fondi
          LV = 0
          DO 72 J = 1,NZ
             IF (MI(I,J).GT.LV) THEN
                LV = MI(I,J)
                FO(I) = (J-1)*BINS+BINS/2.+AMI
             END IF
   72     CONTINUE
   70  CONTINUE
C
       DO 80 J = 1,NRY
          DO 81 I = 1,NRX
             BACKGR(I) = FO((J-1)*NRX+I)
   81     CONTINUE
          CALL TBRWRR(TIDSKY,J,NRX,ICOL,BACKGR,ISTAT)
   80  CONTINUE
C
       CALL STDWRI(IMF,'SKY_AREA',AREA,1,6,KUN,ISTAT)
       CALL TBSINI(TIDSKY,ISTAT)
       CALL TBTCLO(TIDSKY,ISTAT)
       CALL STSEPI
       END
