C @(#)bgval.for 17.1.1.1 (ESO-DMD) 01/25/02 17:40:31 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 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence 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 BGVAL(MODE,AIMG,NP,NL,IAREA,ITER,IMIN,LEAP,BGVLS,RMNS, + NPON,ISTOP,ILAST) C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine BGVAL version 2.00 880630 C K. Banse (cleaned version of H. Waldthausen) C C.KEYWORDS C kap*sig clipping C C.PURPOSE C C this routine computes a constant background value (bgval) of C rows or subimages C C.ALGORITHM C iterarate on bgval using a kap*sig clipping C C.INPUT/OUTPUT parameters C C Using mode 0 or -1 the routine can quickly select the most probable C bgvals of an image if LEAP is set. C C -> MODE operation mode C :0 background of rows C :1 background of an area C -> AIMG input image C -> NP x-pixels and C -> NL y-pixels of AIMG C -> IAREA pixel indices of selected image area : C x-start,x-end,y-start,y-end C -> ITER number of iterations C -> IMIN minimum number of points to derive bgval C -> ISTOP resticts execution of BGVAL to ISTOP C consistent and consecutive background scanlines. C if neg. or 0 bgval is calculated for all rows. C -> LEAP number of rows to jump when the present C row-bgval exeeds the last min-bgval+3*sig threshold. C this will happen if a row has to many/extended sources. C if neg. or 0 no jumping;if ISTOP neg no jumping. C ouput MODE = 0 C <- BGVLS background values of rows C : -99. if outside IAREA or C iteration deleted all points (idap). C <- RMNS standard deviation of bgval points C : -1. if outside IAREA or idap. C <- NPON number of bgval points found C : -1 if outside IAREA or idap. C <- ILAST column index of last computed bgval. C C output MODE = 1 C <-BGVLS(1) background value of area : -99. if idap C <-RMNS(1) standard deviation of bgval points : -1 if idap C <-NPON(1) number of bgval points found : -1 if idap C <- ILAST number of iterations used C C------------------------------------------------------------------------- C IMPLICIT NONE C INTEGER MODE,NP,NL,IAREA(4),ITER,IMIN,LEAP,NPON(*),ISTOP,ILAST INTEGER IXA,IXE,IYA,IYE INTEGER K,KL,JUMP,ISTOR,NPA,J,I INTEGER IX,IY,NPP,NPC,IXP,IYP C REAL BGVLS(*),RMNS(*),AIMG(NP,NL) REAL SM,DUMMY,SMIN,HCUT,AKAP REAL SMEAN,S,SQ,V,AP,DEL,DELS,DELA,DELV,ARG,RM,BG REAL SMSQ,RMS,CLIP C C.. initialize C IXA = IAREA(1) IXE = IAREA(2) IYA = IAREA(3) IYE = IAREA(4) ILAST = 0 JUMP = 0 ISTOR = 0 SM = 0. DUMMY = 1.0E37 SMIN = DUMMY HCUT = DUMMY IF (ISTOP.EQ.0) ISTOP = -1 IF (LEAP.LT.0) LEAP = 0 KL = 1 IF (MODE.EQ.0) KL = NL IF (MODE.EQ.-1) KL = NP C DO 20, K=1,KL BGVLS(K) = -99. RMNS(K) = -1. NPON(K) = -1 20 CONTINUE C IF (MODE.EQ.0) GOTO 2000 IF (MODE.EQ.1) GOTO 4000 C C*** MODE = -1 compute background levels and rms of columns C DO 1000, IX=IXA,IXE IF (ISTOR.EQ.ISTOP) RETURN C ILAST = ILAST + 1 IF (SM.GT.HCUT) THEN IF ( (JUMP.EQ.LEAP) .OR. (ISTOP.LT.0) ) GOTO 50 JUMP = JUMP + 1 GOTO 1000 ENDIF C 50 ISTOR = ISTOR + 1 JUMP = 0 SMEAN = 0. AKAP = 1.E37 NPA = 0 C C.. iterate DO 500, J=1,ITER S = 0. SQ = 0. NPP = 0 DO 200, IY=IYA,IYE V = AIMG(IX,IY) IF (ABS(V-SMEAN) .LE. AKAP) THEN NPP = NPP + 1 S = S + V SQ = SQ + V*V ENDIF 200 CONTINUE C IF (NPP.GT.IMIN) THEN AP = FLOAT(NPP) SMEAN = S/AP DEL = (SQ/AP) - (SMEAN*SMEAN) ARG = MAX(1.E-37,DEL) RM = SQRT(ARG) BG = SMEAN AKAP = RM + RM IF (NPA.EQ.NPP) GOTO 600 NPA = NPP ELSE NPON(IX) = NPP AP = 0. ISTOR = 0 GOTO 1000 ENDIF 500 CONTINUE C 600 SM = SMEAN BGVLS(IX) = BG RMNS(IX) = RM NPON(IX) = NPP IF (SM.LT.SMIN) THEN ISTOR = 0 SMIN = SM HCUT = SMIN + (3.*RM) ENDIF 1000 CONTINUE RETURN C C*** MODE = 0 compute background levels and rms of rows C 2000 DO 3000, IY=IYA,IYE IF (ISTOR.EQ.ISTOP) RETURN ILAST = ILAST + 1 IF (SM.GT.HCUT) THEN IF ( (JUMP.EQ.LEAP) .OR. (ISTOP.LT.0) ) GOTO 2100 JUMP = JUMP + 1 GOTO 3000 ENDIF C 2100 ISTOR = ISTOR + 1 JUMP = 0 SMEAN = 0. AKAP = 1.E37 NPA = 0 C C.. iterate C DO 2600, J=1,ITER S = 0. SQ = 0. NPP = 0 DO 2300, IX=IXA,IXE V = AIMG(IX,IY) IF (ABS(V-SMEAN) .LE. AKAP) THEN NPP = NPP + 1 S = S + V SQ = SQ + V*V ENDIF 2300 CONTINUE C IF (NPP.GT.IMIN) THEN AP = FLOAT(NPP) SMEAN = S/AP DEL = (SQ/AP) - (SMEAN*SMEAN) ARG = MAX(1.E-37,DEL) RM = SQRT(ARG) BG = SMEAN AKAP = RM + RM IF (NPA.EQ.NPP) GOTO 2900 NPA = NPP ELSE NPON(IY) = NPP AP = 0. ISTOR = 0 GOTO 3000 ENDIF 2600 CONTINUE C 2900 SM = SMEAN BGVLS(IY) = BG RMNS(IY) = RM NPON(IY) = NPP IF (SM.LT.SMIN) THEN ISTOR = 0 SMIN = SM HCUT = SMIN + (3.*RM) ENDIF 3000 CONTINUE RETURN C C*** MODE = 1 background value of an area C C.. get initial smean C 4000 SMEAN = 0. DO 4100, IY=IYA,IYE DO 4040 IX=IXA,IXE SMEAN = SMEAN + AIMG(IX,IY) 4040 CONTINUE 4100 CONTINUE C IXP = IXE - IXA + 1 IYP = IYE - IYA + 1 NPC = IXP*IYP SMEAN = SMEAN/FLOAT(NPC) C C.. get initial rms DELS = 0. SMSQ = SMEAN*SMEAN C DO 4300, IY=IYA,IYE DO 4200, IX=IXA,IXE V = AIMG(IX,IY) DELV = V*V + SMSQ - 2.*V*SMEAN IF (DELV.LT.0.) DELV = 0. DELS = DELS + DELV 4200 CONTINUE 4300 CONTINUE C IF (NPC.NE.1) THEN DELA = DELS/FLOAT(NPC-1) RMS = SQRT(DELA) AKAP = 2.0 CLIP = AKAP*RMS ELSE BGVLS(1) = SMEAN RMNS(1) = -1.0 NPON(1) = NPC RETURN ENDIF C C.. exit if no iter BGVLS(1) = SMEAN RMNS(1) = RMS NPON(1) = NPC IF (ITER.LE.0) RETURN C C.. iterate C NPA = NPC DO 5000, I=1,ITER S = 0. NPP = 0 DELS = 0. SMSQ = SMEAN*SMEAN DO 4600, IY=IYA,IYE DO 4500, IX=IXA,IXE V = AIMG(IX,IY) IF ( ABS(V-SMEAN) .LT. CLIP ) THEN NPP = NPP + 1 S = S + V DELV = (V*V) + SMSQ - (2.*SMEAN*V) IF (DELV.LT.0.) DELV = 0. DELS = DELS + DELV ENDIF 4500 CONTINUE 4600 CONTINUE C C.. calculate new CLIP C IF (NPP.GT.IMIN) THEN IF (NPA.EQ.NPP) GOTO 6500 NPA = NPP DELA = DELS/FLOAT(NPP-1) RMS = SQRT(DELA) SMEAN = S/FLOAT(NPP) CLIP = AKAP*RMS ELSE NPON(1) = NPA RMNS(1) = -1. BGVLS(1) = SMEAN RETURN ENDIF 5000 CONTINUE C C.. store results 6500 ILAST = I IF (ILAST.GT.ITER) ILAST = ITER BGVLS(1) = SMEAN RMNS(1) = RMS NPON(1) = NPA C RETURN END SUBROUTINE LINFOL(A,NDIM,START,STEP,TAB) C C++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION C subroutine LINFOL version 3.30 850729 C K. Banse ESO - Garching C C.KEYWORDS C C.PURPOSE C C.ALGORITHM C C.INPUT/OUTPUT C call as LINFOL(A,NDIM,START,STEP,TAB) C C.VERSIONS C C-------------------------------------------------- C IMPLICIT NONE C INTEGER NDIM INTEGER N,NEXT,NROW,NCOL,NMAL INTEGER STAT INTEGER TID,TBCLNM(2) INTEGER TABNUL(2) C CHARACTER*(*) TAB CHARACTER CBUF*80 CHARACTER TLABL(2)*16 C DOUBLE PRECISION STEP,START DOUBLE PRECISION RM,XA,XB,YA,YB,XD,XN C REAL A(NDIM) REAL RBUF(2) C INCLUDE 'MID_INCLUDE:ST_DEF.INC' C DATA TLABL /'X ','Y '/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C init XN = START NEXT = 1 C C table CALL TBTOPN(TAB,'I ',TID,STAT) CALL TBIGET(TID,NCOL,NROW,N,N,N,STAT) !get no. of rows + columns IF (NROW.LE.1) THEN CALL STTPUT('Not enough rows in table...',STAT) GOTO 8000 ENDIF CALL TBLSER(TID,TLABL(1),TBCLNM(1),STAT) IF (TBCLNM(1).LE.0) THEN CBUF(1:) = 'column :'//TLABL(1)//'is missing ... ' CALL STTPUT(CBUF,STAT) GOTO 8000 ENDIF CALL TBLSER(TID,TLABL(2),TBCLNM(2),STAT) IF (TBCLNM(2).LE.0) THEN CBUF(1:) = 'column :'//TLABL(2)//'is missing ... ' CALL STTPUT(CBUF,STAT) GOTO 8000 ENDIF C C get first row of values NMAL = 1 CALL TBRGET(TID,'R*4 ',NMAL,2,TBCLNM,RBUF,TABNUL,STAT) XA = RBUF(1) YA = RBUF(2) C 500 A(NEXT) = YA NEXT = NEXT + 1 IF (NEXT.GT.NDIM) GOTO 7000 !already all points filled up C XN = XN + STEP IF (XN.LE.XA) GOTO 500 !set all pixels before XA to YA C C get next row of values 1000 NMAL = NMAL + 1 IF (NMAL.GT.NROW) GOTO 5000 C CALL TBRGET(TID,'R*4 ',NMAL,2,TBCLNM,RBUF,TABNUL,STAT) XB = RBUF(1) YB = RBUF(2) XD = XB - XA IF (XB.LT.XN) THEN !we need XA < XN < XB IF (XD.LT.1.D-30) THEN WRITE(*,10001) NMAL !should be XB > XA ELSE XA = XB YA = YB ENDIF GOTO 1000 !omit that pair X,Y ENDIF RM = (YB - YA) / XD C C here we loop and fill the data 1500 A(NEXT) = SNGL(RM*(XN - XA) + YA) NEXT = NEXT + 1 IF (NEXT.GT.NDIM) GOTO 7000 !already all points filled up C XN = XN + STEP IF (XN.LE.XB) GOTO 1500 C C move to next pair XA = XB YA = YB GOTO 1000 C C here, if table processed 5000 IF (NEXT.LT.NDIM) THEN !missing values set to last one DO 5500, N=NEXT,NDIM A(N) = A(NEXT-1) 5500 CONTINUE ENDIF C C that's it folks 7000 CALL TBTCLO(TID,STAT) RETURN C C Problems with table 8000 DO 8080, N=1,NDIM A(N) = 0.0 8080 CONTINUE RETURN C 10001 FORMAT('invalid x,y value in row ',I4) END