C @(#)rfotsky.for 17.1.1.1 (ES0-DMD) 01/25/02 17:18:17 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