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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.COPYRIGHT   (C) 1992 European Southern Observatory
C.IDENT       .for
C.AUTHOR      E. Oliva,  Firenze-Arcetri
C.KEYWORDS    Spectroscopy, IRSPEC
C
C.PURPOSE     Execute command  ...
C
C.ALGORITHM
C
C
C.INPUT/OUTPUT
C
C.VERSION     1.0    Creation     02.09.1992   E. Oliva
C 030714	last modif
C-------------------------------------------------------
C
      PROGRAM DEFBAD
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      CHARACTER*60 FRAMEI,TABLEO
      CHARACTER*1 SMODE
      CHARACTER*64 CUNIT
      CHARACTER*72 IDENT
      REAL RTH(2)
      INTEGER*8 INPNTR,OUPNTR
      INTEGER NPIX(2)
      DOUBLE PRECISION START(2),STEP(2)
      COMMON /VMR/ MADRID(1)

      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'

      DATA MAXDIM/2/

      IRET=1
      CALL STSPRO('DEFBAD')
C
C GET NAME OF INPUT FRAME AND MAP IT
C CHECK ALSO WHETHER IT IS A 2D FRAME.....
C
      CALL STKRDC('framei',1,1,60,IRET,FRAMEI,KUNIT,KNUL,ISTAT)
      CALL STIGET(FRAMEI,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,MAXDIM,
     ,            NAXIS,NPIX,START,STEP,
     ,            IDENT,CUNIT,INPNTR,NIN,ISTAT)
      IF(NAXIS.NE.2) 
     +     CALL STETER(1,'Input frame must be 2-D')
C
C Create dummy image wich will contain badpixels.
C
      CALL STIPUT('middumma',D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     ,            NAXIS,NPIX,START,STEP,
     ,            IDENT,CUNIT,OUPNTR,NOU,ISTAT)
C
C GET FLAG TO DEFINE the selection mode
C
      CALL STKRDC('smode',1,1,1,IRET,SMODE,KUNIT,KNUL,ISTAT)
C
C Get threshold value(s)
C
      CALL STKRDR('thres',1,2,IRET,RTH,KUNIT,KNUL,ISTAT)
C
C GET KEY TO DEFINE the number of subframes
C
      CALL STKRDI('nsub',1,1,IRET,NSUB,KUNIT,KNUL,ISTAT)
C
C GET THE NAME OF THE OUTPUT TABLE WHICH WILL CONTAIN THE COORDINATES OF THE
C BAD PIXELS 
C
      CALL STKRDC('tableo',1,1,60,IRET,TABLEO,KUNIT,KNUL,ISTAT)
C
C OPEN TABLE FOR UPDATE (IOPEN=2)
C
      IOPEN=2
      CALL TBTOPN(TABLEO,IOPEN,IDETAB,ISTAT)
C
C GET NUMBERS CORRESPONDING TO COLUMNS
C
      CALL TBLSER(IDETAB,'X_COORD',NXCOL,ISTAT)
      CALL TBLSER(IDETAB,'Y_COORD',NYCOL,ISTAT)
      CALL TBLSER(IDETAB,'AVALUE',NAVAL,ISTAT)
      CALL TBLSER(IDETAB,'RVALUE',NRVAL,ISTAT)
C
C CALL ROUTINE WHICH DEFINE BAD_PIXELS AND WRITE THEM IN THE TABLE
C
      CALL FIND(MADRID(INPNTR),MADRID(OUPNTR),
     ,          NPIX(1),NPIX(2),NSUB,
     ,          IDETAB,NXCOL,NYCOL,NAVAL,NRVAL,SMODE,RTH)
C
C Close table
C
      CALL TBTCLO(IDETAB,ISTAT)
C
C RELEASE FILES, UPDATE KEYWORDS AND EXIT
C
      CALL STSEPI
      END
C
C
C
      SUBROUTINE FIND (A,B,NX,NY,NSUB,
     ,                 IDETAB,NXCOL,NYCOL,NAVAL,NRVAL,SMODE,RTH)
C
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)
C
      CHARACTER*1 SMODE,DMODE(2),CMODE(2)
      DIMENSION A(NX,NY),B(NX,NY),RTH(2)
      CHARACTER*200 OUTPUT
C
      DATA DMODE/'r','a'/
      DATA CMODE/'R','A'/
C
C Check INPUT IMAGE (A) 
C
      SUM=0.0
      SUM2=0.0
      DO J=1,NY
        DO I=1,NX
          SUM=SUM+A(I,J)
          SUM2=SUM2+A(I,J)*A(I,J)
        ENDDO
      ENDDO
      STD=SUM2-SUM*SUM/(NX*NY)
      IF(STD.LE.0) THEN
        WRITE(OUTPUT,192)
192     FORMAT('Input image is a constant, check it please')
        CALL STTPUT(OUTPUT,ISTAT)
        RETURN
      ENDIF
C
      DO I=1,2
        IMODE=I
        IF(SMODE.EQ.DMODE(I) .OR. SMODE.EQ.CMODE(I)) GO TO 10
      ENDDO
      CALL STETER
     +(1,'Invalid selection mode, only -r-, -a- are recognized')
10    CONTINUE

C
C Set badpixel_image (B) to zero
C
      DO J=1,NY
        DO I=1,NX
          B(I,J)=0.
        ENDDO
      ENDDO
      NBAD=0
C
C Divide image in 16 quadrants and define standard deviation in
C each of them.
C
      IX1=1
      IX2=NX/NSUB
      DO IXQ=1,NSUB
        IY1=1
        IY2=NY/NSUB
        DO IYQ=1,NSUB
          CALL STD_ITER(A,NX,NY,IX1,IX2,IY1,IY2,AVE,STD)
          DO IY=IY1,IY2
            DO IX=IX1,IX2
              VAL=A(IX,IY)
              RVAL=(VAL-AVE)/STD
              IF(IMODE.EQ.1) THEN
                IF(ABS(RVAL).GT.RTH(1)) THEN
                  NBAD=NBAD+1
                  R4=IX
                  CALL TBRWRR(IDETAB,NBAD,1,NXCOL,R4,ISTAT)
                  R4=IY
                  CALL TBRWRR(IDETAB,NBAD,1,NYCOL,R4,ISTAT)
                  CALL TBRWRR(IDETAB,NBAD,1,NAVAL,VAL,ISTAT)
                  CALL TBRWRR(IDETAB,NBAD,1,NRVAL,RVAL,ISTAT)
                  B(IX,IY)=A(IX,IY)
                ENDIF
              ENDIF
              IF(IMODE.EQ.2) THEN
                IF(VAL.LT.RTH(1).OR.VAL.GT.RTH(2)) THEN
                  NBAD=NBAD+1
                  R4=IX
                  CALL TBRWRR(IDETAB,NBAD,1,NXCOL,R4,ISTAT)
                  R4=IY
                  CALL TBRWRR(IDETAB,NBAD,1,NYCOL,R4,ISTAT)
                  CALL TBRWRR(IDETAB,NBAD,1,NAVAL,VAL,ISTAT)
                  CALL TBRWRR(IDETAB,NBAD,1,NRVAL,RVAL,ISTAT)
                  B(IX,IY)=A(IX,IY)
                ENDIF
              ENDIF
            ENDDO
          ENDDO
          IY1=IY2+1
          IY2=IY1+NY/NSUB
          IY2=MIN(IY2,NY)
        ENDDO
        IX1=IX2+1
        IX2=IX1+NX/NSUB
        IX2=MIN(IX2,NX)
      ENDDO
      WRITE(OUTPUT,195) NBAD
195   FORMAT(I5,' bad pixels found and stored in table.')
      CALL STTPUT(OUTPUT,ISTAT)
      RETURN
      END



      SUBROUTINE STD_ITER(A,NX,NY,IX1,IX2,IY1,IY2,AVE,STD)
C
C Define the std of a given region of the image (IX1,IX2:IY1,IY2)
C by iteratively eliminating the points which are more than 3 
C sigma from the average
C
      IMPLICIT REAL(A-H,O-Z)
      IMPLICIT INTEGER(I-N)

      DIMENSION A(NX,NY)
      DIMENSION IXOUT(5000),IYOUT(5000)

      NOUT=0
      IXOUT(1)=0
      IYOUT(1)=0

10    CONTINUE
      SUM=0.0
      SUM2=0.0
      NAVE=0
      DO IY=IY1,IY2
        DO IX=IX1,IX2
          DO J=1,NOUT
            IF(IX.EQ.IXOUT(J) .AND. IY.EQ.IYOUT(J)) GO TO 11
          ENDDO
          VAL=A(IX,IY)
          SUM=SUM+VAL
          SUM2=SUM2+VAL*VAL
          NAVE=NAVE+1
11        CONTINUE
        ENDDO
      ENDDO
      IF(NAVE.LE.0) THEN
        CALL STTPUT('---------------------------------------',ISTAT)
        CALL STTPUT('  NAVE=0, problems in STD_ITER',ISTAT)
        CALL STTPUT(' You probably chose too many sub_frames.',ISTAT)
        CALL STTPUT('---------------------------------------',ISTAT)
        CALL STETER(1,'Aborted')
      ENDIF
      AVE=SUM/NAVE 
      STD=SQRT((SUM2-NAVE*AVE*AVE)/(NAVE-1))
      IF(STD.LE.0) CALL STETER(1,'STD<=0, problems in STD_ITER')
      N=0
      DO IY=IY1,IY2
        DO IX=IX1,IX2
          DIST=ABS((A(IX,IY)-AVE)/STD)
          IF(DIST.GT.3) THEN
            N=N+1
            IXOUT(N)=IX
            IYOUT(N)=IY
          ENDIF
        ENDDO
      ENDDO
      IF(N.GT.NOUT) THEN
        NOUT=N
        GO TO 10
      ENDIF
      RETURN
      END
