C @(#)necfilt.for	17.1.1.1 (ES0-DMD) 01/25/02 17:51:29
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 ECHFILT
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  1 JUL 1989
C
C.LANGUAGE: F77+ESOext
C
C.IDENTIFICATION
C
C  Program ECHFILT
C
C  Otmar Stahl                     Landessternwarte Heidelberg
C  M Peron
C.KEYWORDS
C
C  ECHELLE, CASPEC, COSMIC RAY FILTER
C
C.PURPOSE
C
C  FILTER ECHELLE IMAGES FOR COSMIC RAY HITS AND OTHER DEFECTS. ORDER 
C  POSITIONS ARE DEFINED  BY COEFFICIENTS OF A REGRESSION. THE OFFSET 
C  IS DEFINED BY THE PARAMETER OFFSET.
C            
C
C.ALGORITHM
C   
C  TWO FILTERING METHODS ARE COMBINED: OUTSIDE THIS PROGRAM, THE RAW
C  DATA FRAME IS FILTERED WITH A STANDARD MEDIAN FILTER. THE BACKGROUND
C  IS DETERMINED ON THE FILTERED FRAME AND SUBTRACTED FROM BOTH FRAMES. 
C  BOTH THE UNFILTERED FRAME AND THE MEDIAN FILTERED FRAME ARE USED AS INPUT
C  FOR THIS PROGRAM. SINCE MEDIAN  FILTERING IS NOT GOOD CLOSE TO THE 
C  SPECTRAL ORDERS, THESE REGIONS ARE REPLACED (IN THE MEDIAN FILTERED
C  FRAME) WITH RAW DATA FILTERED THIS PROGRAM. THE METHOD MAKES USE OF
C  THE FACT, THAT IN ECHELLE FRAMES THE PROFILE PERPENDICULAR TO THE 
C  MAIN DISPERSION VARIES USALLY ONLY SLOWLY WITH WAVELENGTH. THIS 
C  PROFILE IS FOUND BY MEDIAN FILTERING A NUMBER OF  SPATIAL PROFILES
C  IN THE NEIGHBORHOOD. EACH INDIVIDUAL PROFILE IS THEN COMPARED PIXEL
C  BY PIXEL WITH THE MEDIAN PROFILE. OUTLIERS ABOVE A GIVEN THRESHOLD 
C  ARE DETECTED AND REMOVED BY REPLACING THEM WITH THE APPROPRIATE VALUE
C  FROM THE MEDIAN PROFILE. 
C
C
C  C.INPUT/OUTPUT
C
C  FILTER/ECHELLE output input widx,widy,no_iter ron,g,ethresh 
C                                             deg1,deg2 radx,rady,mthresh
C
C  IN_A               input frame
C  OUT_A               output frame 
C  INPUTI              wx,wy,iter
C  INPUTR              ron,g,ethresh
CC
C
C-----------------------------------------------------------


        IMPLICIT NONE
        INTEGER                I,IORD,IORD1,IORD2,NCOL,NROW,NOB
        INTEGER                INDX,NA,NPAR,ISTAT,ICOSM,NSORT,NOA
        INTEGER*8   PNTRA,PNTRB
        INTEGER                NDIM,NAXISA,STATUS,IDUMMY
        INTEGER                NPIXA(3), IPAR(10),JPAR(10),KUN,KNUL,IACT
        INTEGER                TID,ACOL,AROW,MASK,MADRID,MAXDIM
        PARAMETER              (NDIM = 3)
        CHARACTER*60              FRAMEIN, FRAMEOUT, INA1, INA2, TABLE
        CHARACTER*8               TYPE
        CHARACTER*9               NAME
        CHARACTER*60              TABLEFULL
        CHARACTER*1 DSC
        REAL                        RPAR(10)
        REAL                       PARAMS(12)
        DOUBLE PRECISION           DPAR(100)
        CHARACTER*72              INSTRUMA
        CHARACTER*64              CUNITA, DESNAM
        DOUBLE PRECISION          STARTA(3),STEPA(3)        
        INCLUDE 'MID_INCLUDE:TABLES.INC'
        COMMON /VMR/ MADRID(1)
        INCLUDE 'MID_INCLUDE:TABLED.INC'
        DATA MASK/3/
c
        CALL STSPRO('ECHFILT')
C
C  get input
C
        CALL STKRDC('IN_A', 1,1,60,IACT,INA1,KUN,KNUL,STATUS)
        CALL STKRDC('IN_B', 1,1,60,IACT,TABLE,KUN,KNUL,STATUS)
        CALL STKRDC('INPUTC', 1,1,9,IACT,NAME,KUN,KNUL,STATUS)
        CALL STKRDC('OUT_A',1,1,60,IACT,INA2,KUN,KNUL,STATUS)
        CALL STKRDI('INPUTI',1,3,IACT,JPAR,KUN,KNUL,STATUS)
        CALL STKRDR('INPUTR',1,3,IACT,RPAR,KUN,KNUL,STATUS)       
C
C CHECK WINDOW SIZE
C
        IF(JPAR(1).GT.21)      JPAR(1) = 21
        IF(JPAR(2).GT.21)      JPAR(2) = 21
        IF(JPAR(3).LE.JPAR(2)) JPAR(3) = JPAR(2)-1
C
       CALL CLNFRA(INA1,FRAMEIN,0)
       CALL CLNFRA(INA2,FRAMEOUT,0)
       MAXDIM = 3
       DSC = ' '
C
C  map input and output frame
C
      CALL STIGET(FRAMEIN,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,MAXDIM,
     1    NAXISA,NPIXA,STARTA,STEPA,INSTRUMA,CUNITA,PNTRA,NOA,STATUS)
       CALL STIGET(FRAMEOUT,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,MAXDIM,
     1   NAXISA,NPIXA,STARTA,STEPA,INSTRUMA,CUNITA,PNTRB,NOB,STATUS)

C  input parameters in table
C
       I = INDEX(TABLE,'.') 
       IF (I.EQ.0) THEN
         TABLE     = TABLE(1:59)//' '
         IDUMMY    = INDEX(TABLE,' ') - 1
         TABLEFULL = TABLE(1:IDUMMY)//'.tbl'
       ELSE
         TABLEFULL = TABLE
         TABLE     = TABLEFULL(1:I-1)
       ENDIF
       CALL TBTOPN(TABLE,F_IO_MODE,TID, STATUS)
        CALL TBIGET(TID,NCOL,NROW,NSORT,ACOL,AROW,STATUS) 
c
C  read coeffs
C
       NAME = NAME(1:8)//' '
       INDX = INDEX(NAME,' ')-1
       DESNAM = NAME(1:INDX)//'C'
        CALL STDRDC(TID,DESNAM,1,17,4,IACT,TYPE,KUN,KNUL,STATUS)
       TYPE = 'MULT'
       IF (TYPE.EQ.'MULT') THEN
         DESNAM = NAME(1:INDX)//'I'
          CALL STDRDI(TID,DESNAM,4,4,IACT,IPAR,KUN,KNUL,STATUS)
         NA = (IPAR(3)+1)*(IPAR(4)+1)
         DESNAM = NAME(1:INDX)//'R'
          CALL STDRDR(TID,DESNAM,1,4,IACT,PARAMS(4),KUN,KNUL,STATUS)
         DESNAM = NAME(1:INDX)//'D'
          CALL STDRDD(TID,DESNAM,1,NA,IACT,DPAR,KUN,KNUL,STATUS)
         CALL TBTCLO(TID, STATUS)
C
         IORD1 = PARAMS(6)
         IORD2 = PARAMS(7)
c
C  loop over the orders
C
         NPAR      = 10
         PARAMS(8) = IPAR(3)
         PARAMS(9) = IPAR(4)
          ICOSM     = 0
         DO 20 I = IORD1, IORD2
           PARAMS(10) = I 
           IORD = I - IORD1 + 1
           CALL EXTR(MADRID(PNTRA), NPIXA(1), NPIXA(2), 
     .      MADRID(PNTRB),IORD,NPAR,
     .      PARAMS,NA,DPAR,JPAR,RPAR,ICOSM) 
20     CONTINUE
       ELSE
         CALL STTPUT(' Error in dispersion coefficients ', ISTAT)
         GOTO 1000
       ENDIF
c
C
C  descriptor handling
C
       CALL STDCOP(NOA,NOB,MASK,DSC,STATUS)
C
C  free data
C
1000       CALL STSEPI
       END


       SUBROUTINE EXTR(IN,NPIXA1,NPIXA2,OUT,INDX,
     .  NPAR,PAR,ND,DPAR,JPAR,RPAR,ICOSM)
C
        IMPLICIT NONE
        INTEGER               JPAR(10),npixa1,npixa2,iwind,iwindh,nww
        INTEGER               iwidth,iter,iylim,k,l,ix1,ip,l1,k1,iy
        INTEGER               iy2,indx,icosm,ii,ixx,iyy,nd,npar,jj
        DOUBLE PRECISION      DPAR(ND)
        REAL                       OUT(NPIXA1,NPIXA2),PAR(NPAR)
        REAL                    SUBFR(21,21),FILT(21)
        REAL                       IN(NPIXA1,NPIXA2),RPAR(10)
        DOUBLE PRECISION           XVAL, YVAL, AX
        DOUBLE PRECISION           WIDTH,HWIDTH,OFFSET,DC
        DOUBLE PRECISION           YV,WW(50),SUM
       NWW = 50
C
C ... define extraction constants
C
        IWIND  = JPAR(1)
        IWIDTH = JPAR(2)
        ITER   = JPAR(3)
       WIDTH  = IWIDTH
       HWIDTH = 0.5*WIDTH
        IWINDH = IWIND/2
C
C       OFFSET = PAR(3) - HWIDTH + 1.
       OFFSET = PAR(3) - HWIDTH + 0.5
       IYLIM  = NPIXA2 - IWIDTH - 1
        ax     = 1.d0
C
C ... get ordinate values
C
       K     = PAR(8)
       L     = PAR(9)
       YV    = PAR(10)
C
C ... follow the order with the filter box
C
         DO IX1 = 1+IWINDH,NPIXA1-IWINDH
           IP  = 0
           DC  = 1.D0
           YVAL= 0.D0
            XVAL= ax
           DO L1 = 0, L
             IP = IP + 1
             WW(IP) = DC
             YVAL   = YVAL + WW(IP)*DPAR(IP)
             DO K1 = 1, K
               IP = IP + 1
               WW(IP) = WW(IP-1)*XVAL
               YVAL   = YVAL + WW(IP)*DPAR(IP)
             ENDDO
             DC = DC*YV
           ENDDO
            YVAL = YVAL + OFFSET
           SUM = 0.D0
C
            IY = YVAL
            IF (IY.GT.2 .AND. IY.LT.IYLIM) THEN
              IY2 = IY + IWIDTH
              IF(IY2.GT.2 .AND. IY2.LT.IYLIM) THEN
                DO II = 1,IWIND
                  IXX = IX1-IWINDH+(II-1)
                  DO JJ = 1,IWIDTH
                    IYY = IY+JJ-1
                    SUBFR(II,JJ) = IN(IXX,IYY)
                  ENDDO
                ENDDO
                CALL FILTER(SUBFR,FILT,JPAR,RPAR,ICOSM)
                DO JJ = 1,IWIDTH
                  IYY = IY+JJ-1
                  OUT(IX1,IYY) = filt(JJ)
                ENDDO
              ENDIF
            ENDIF
            AX = AX + 1.D0
c
         ENDDO
       RETURN
       END

      SUBROUTINE FILTER(SUBFR,Y,JPAR,RPAR,ICOSM)
C
      IMPLICIT NONE
      INTEGER         JPAR(10),ICOSM,IWIND,IWIDTH,ITER,J,K,IX,I,L,KK
      INTEGER*2       MASK(21)
      REAL            RPAR(10),RON,G,SIGMA,SIGMA2,V0
      REAL            SUBFR(21,21),VARI(21),PROF(21)
      REAL            SPEC(21),Y(21),X(21),XMED
      REAL            DIFF,DIFF2,QUOT,SUM,SUM1,SUM2,RMAX
C
      RON    = RPAR(1)
      G      = RPAR(2)
      SIGMA  = RPAR(3)
      IWIND  = JPAR(1)
      IWIDTH = JPAR(2)
      ITER   = JPAR(3)
      SIGMA2 = SIGMA*SIGMA
      V0     = RON**2/G**2
C
C SAVE COLUMN AT THE CENTER
C
      IX = IWIND/2 + 1
      DO I = 1,IWIDTH
        Y(I) = SUBFR(IX,I)
      ENDDO
C
C  EXTRACT INITIAL SPECTRUM (JUST THE SUM) AND SET mask TO 1
C
      DO J = 1,IWIND
        SUM = 0.
        DO K = 1,IWIDTH
          SUM = SUM+SUBFR(J,K)
        ENDDO
        SPEC(J) = SUM
        IF(SPEC(J).LT.1.) SPEC(J) = 1.
      ENDDO
C
      DO K = 1,IWIDTH
        MASK(K) = 1
      ENDDO
C
C COMPUTE SPATIAL PROFILE
C                                                 
      DO J = 1,IWIND
        DO K = 1,IWIDTH
          SUBFR(J,K) = SUBFR(J,K)/SPEC(J)
        ENDDO
      ENDDO 
c
      DO K = 1,IWIDTH
        DO J = 1,IWIND
          X(J) = SUBFR(J,K)
        ENDDO
C
C  MEDIAN FILTER OF ROWS
C
       CALL MDIAN1(X,IWIND,XMED)
       PROF(K) = XMED
       IF(PROF(K).LT.1.E-4) PROF(K)=1.E-4
      ENDDO
C
C NORMALIZE SPATIAL PROFILE AND ENFORCE POSITIVITY
C
      DO J = 1,IWIND
        SUM=0.
        DO K = 1,IWIDTH
          SUM = SUM+PROF(K)
        ENDDO
        DO K = 1,IWIDTH
          PROF(K) = PROF(K)/SUM
        ENDDO
      ENDDO 
C
c iterate (IN THE PRESENT VERSION, THE SPATIAL PROFILE IS NOT 
C ITERATED, BUT ONLY THE MASK AND THE VARIANCE)
c                                                                       
      DO L = 1,ITER                                                  
c                                                                       
c  compute variance estimate                 
C
        DO K = 1,IWIDTH 
          VARI(K) = V0 + ABS(SPEC(IX)*PROF(K))/G
        ENDDO
c                                                                       
c MASK COSMIC RAY HITS (ONLY ONE PER COLUMN AND ITERATION)
c
        RMAX = 1.
        KK = 0
        DO K = 1,IWIDTH
          DIFF = Y(K)-SPEC(IX)*PROF(K)
          DIFF2 = DIFF*DIFF
          QUOT = (DIFF2/(SIGMA2*VARI(K)))*MASK(K)
          IF(QUOT.GT.RMAX) THEN
            ICOSM = ICOSM+1
            RMAX = QUOT
            KK = K
          ENDIF
        ENDDO 
        IF(KK.NE.0) THEN
          MASK(KK) = 0
        ENDIF
C
C COMPUTE WEIGHTED MEAN WITH COSMICS MASKED
C
        SUM1 = 0.
        SUM2 = 0.            
        DO K = 1,IWIDTH
          SUM1 = SUM1 + MASK(K)*PROF(K)*Y(K)   /VARI(K)
          SUM2 = SUM2 + MASK(K)*PROF(K)*PROF(K)/VARI(K)
        ENDDO
        SPEC(IX) = SUM1/SUM2
      ENDDO
C
C FINALLY, REMOVE COSMICS
C
      DO K = 1,IWIDTH
        IF(MASK(K).EQ.0) THEN
          Y(K) = SPEC(IX)*PROF(K)
          ELSE 
          Y(K) = Y(K)
        ENDIF
      ENDDO      
c                                                                       
      RETURN
      END                                                           
C
C*************************************************************************
      SUBROUTINE MDIAN1(ARR,N,XMED)
C
      IMPLICIT NONE
      INTEGER      N,N2
      REAL         ARR(N),XMED
C
C  SORT DATA
C
      CALL SORT(N,ARR)
C
C  COMPUTE MEDIAN
C
      N2 = N/2
      IF(2*N2.EQ.N) THEN
        XMED = 0.5*(ARR(N2)+ARR(N2+1))
      ELSE
        XMED = ARR(N2+1)
      ENDIF
C
      RETURN
      END


      SUBROUTINE SORT(N,ARR)
C
      IMPLICIT NONE
      INTEGER       I,J,N
      REAL          ARR(N),A
C
      DO J = 2,N
        A = ARR(J)      
        DO I = J-1,1,-1
          IF(ARR(I).LE.A) GOTO 10
          ARR(I+1) = ARR(I)
        ENDDO
      I = 0
10    ARR(I+1) = A
      ENDDO
C
      RETURN
      END