C @(#)bgval.for	19.1 (ESO-DMD) 02/25/03 14:01:03
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
