C @(#)necavmed.for	19.1 (ES0-DMD) 02/25/03 14:20:21
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
C @(#)necavmed.for	19.1 (ESO)   02/25/03    14:20:21
C
C+
C.COPYRIGHT: Copyright (c) 1991 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0          23-JULY-1991
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: P.BALLESTER
C
C.IDENTIFICATION
C
C.KEYWORDS
C
C  ECHELLE, CASPEC, BLAZE FUNCTION
C
C.PURPOSE
C
C  COMPARE SUCCESSIVE ORDERS TO GET THE MEDIAN AND AVERAGE BLAZE FUNCTION
C  
C.ALGORITHM
C
C  IN DEVELOPMENT ...
C
C.INPUT/OUTPUT
C
C-
      PROGRAM AVEMED

      IMPLICIT  NONE

      INTEGER   NAXISA,NPIXA(2),IAV,STAT,ACTVAL,MAXORD
      INTEGER*8   PNTRA,PNTRB
      INTEGER   IMNOA,IMNOB
      INTEGER   KNULL,KUNIT(1)
      INTEGER   MADRID(1)

      PARAMETER (MAXORD=500)

      CHARACTER FRAMEA*60,FRAMEB*60
      CHARACTER CUNIT*64,IDENTA*72

      INTEGER          ORDSTA(MAXORD),ORDEND(MAXORD),DELTAY
      REAL             THRES
      DOUBLE PRECISION STEPA(2),STARTA(2)

      INCLUDE 'MID_INCLUDE:st_def.inc'

      COMMON /VMR/ MADRID

      INCLUDE 'MID_INCLUDE:st_dat.inc'

      CALL STSPRO('AVEMED')


      CALL STKRDC('IN_A',1,1,60,IAV,FRAMEA,KUNIT,KNULL,STAT)

      CALL STIGET(FRAMEA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     +            2,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNIT,
     +            PNTRA,IMNOA,STAT)

      IF (NPIXA(2).LT.MAXORD) THEN

          CALL STDRDI(IMNOA,'ORDSTA',1,NPIXA(2),ACTVAL,ORDSTA,
     +                                       KUNIT,KNULL,STAT)
          CALL STDRDI(IMNOA,'ORDEND',1,NPIXA(2),ACTVAL,ORDEND,
     +                                       KUNIT,KNULL,STAT)

      ELSE

          CALL STETER(10,'Buffer overflow in AVEMED.')

      ENDIF

      CALL STKRDC('OUT_A',1,1,60,IAV,FRAMEB,KUNIT,KNULL,STAT)

      CALL STIPUT(FRAMEB,D_R4_FORMAT,F_IO_MODE,F_IMA_TYPE,
     +            NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNIT,
     +            PNTRB,IMNOB,STAT)

      CALL STKRDI('INPUTI',1,1,IAV,DELTAY,KUNIT,KNULL,STAT)


      IF (DELTAY.GT.0) THEN

         CALL STKRDR('INPUTR',1,1,IAV,THRES,KUNIT,KNULL,STAT)

         CALL MEDIAN(MADRID(PNTRA),NPIXA(1),NPIXA(2),MADRID(PNTRB),
     +            ORDSTA,ORDEND,DELTAY,THRES)

      ELSE 

         CALL AVERAGE(MADRID(PNTRA),NPIXA(1),NPIXA(2),MADRID(PNTRB),
     +            ORDSTA,ORDEND)

      ENDIF

      CALL STSEPI

      END

C ======================Median Filtering====================

      SUBROUTINE MEDIAN(INPFRAM,NX,NY,OUTFRAM,ORDSTA,ORDEND,
     +                  DELY,THRES)

      IMPLICIT   NONE

      INTEGER    NX,NY,DELY,ROW,COL
      INTEGER    ORDSTA(NY),ORDEND(NY)

      REAL       INPFRAM(NX,NY),OUTFRAM(NX,NY),THRES,RATIO

      INTEGER    STAMIN,STAMAX,ENDMIN,ENDMAX
      INTEGER    ROWINF,ROWSUP,FLAG,ROWCUR,Y
      REAL       TRANS,BUFFER(-50:50)

      IF (DELY.GT.50) THEN

         CALL STETER(9,'Width too big in AVEMED. Limited to 50.')

      ENDIF

C --- Initialisations

      STAMIN = ORDSTA(1)
      STAMAX = ORDSTA(1)
      ENDMIN = ORDEND(1)
      ENDMAX = ORDEND(1)
      ROWINF = 0
      ROWSUP = 0

      DO 5  ROW = 1,NY

             IF (ORDSTA(ROW).LT.STAMIN) STAMIN = ORDSTA(ROW)
             IF (ORDSTA(ROW).GT.STAMAX) STAMAX = ORDSTA(ROW)
             IF (ORDEND(ROW).LT.ENDMIN) ENDMIN = ORDEND(ROW)
             IF (ORDEND(ROW).GT.ENDMAX) ENDMAX = ORDEND(ROW)

5     CONTINUE

      DO 10  COL = 1,NX

C --- Find out what is ROWINF and ROWSUP for column COL
C --- This part is easier to program with a DO WHILE. Consider C language.

      IF (COL.GE.STAMAX.AND.COL.LE.ENDMIN) THEN

          ROWINF = 1
          ROWSUP = NY

      ELSE IF (COL.LT.STAMAX) THEN

          ROWINF = 0
          ROWSUP = 0
          DO 20  ROW = 1,NY
             IF (ORDSTA(ROW).LT.COL.AND.ROWINF.EQ.0)  ROWINF = ROW
20        CONTINUE
          DO 25  ROW = NY,1,-1
             IF (ORDSTA(ROW).LT.COL.AND.ROWSUP.EQ.0)  ROWSUP = ROW
25        CONTINUE

      ELSE IF (COL.GT.ENDMIN) THEN

          ROWINF = 0
          ROWSUP = 0
          DO 30  ROW = 1,NY
             IF (ORDEND(ROW).GT.COL.AND.ROWINF.EQ.0)  ROWINF = ROW
30        CONTINUE
          DO 35  ROW = NY,1,-1
             IF (ORDEND(ROW).GT.COL.AND.ROWSUP.EQ.0)  ROWSUP = ROW
35        CONTINUE

      ENDIF

C --- Fill in the buffer with defined values, prolonging the frame by a
C --- mirror effect.

      DO 40 ROW = ROWINF,ROWSUP

         DO 50 Y = (-1)*DELY , DELY

             ROWCUR = ROW+Y

             IF (ROWCUR.LT.ROWINF)
     +               BUFFER(Y) = INPFRAM(COL,2*ROWINF-ROWCUR)

             IF (ROWCUR.GT.ROWSUP)
     +           BUFFER(Y) = INPFRAM(COL,2*ROWSUP-ROWCUR)

             IF (ROWCUR.GE.ROWINF.AND.ROWCUR.LE.ROWSUP)
     +           BUFFER(Y) = INPFRAM(COL,ROWCUR)

50       CONTINUE

C ---    Sort the buffer values and get the median value.
C ---    Consider using HeapSort algorithm (Numerical Recipes)

60       CONTINUE

         FLAG = 0

         DO 70 Y = (-1)*DELY , DELY-1

            IF (BUFFER(Y).GT.BUFFER(Y+1)) THEN
                TRANS = BUFFER(Y)
                BUFFER(Y) = BUFFER(Y+1)
                BUFFER(Y+1) = TRANS
                FLAG = 1
            ENDIF

70       CONTINUE

         IF (FLAG.NE.0) GOTO 60

         RATIO = (INPFRAM(COL,ROW)-BUFFER(0))/(INPFRAM(COL,ROW)
     +                                                   +BUFFER(0))

         IF (ABS(RATIO).GT.THRES)   THEN
                 OUTFRAM(COL,ROW) = BUFFER(0)
         ELSE
                 OUTFRAM(COL,ROW) = INPFRAM(COL,ROW)
         ENDIF

40    CONTINUE
10    CONTINUE

      RETURN
      END

C ====================== Average ====================

      SUBROUTINE AVERAGE(INPFRAM,NX,NY,OUTFRAM,ORDSTA,ORDEND)

      IMPLICIT   NONE

      INTEGER    NX,NY,ROW,COL
      INTEGER    ORDSTA(NY),ORDEND(NY)

      REAL       INPFRAM(NX,NY),OUTFRAM(NX,NY)

      INTEGER    STAMIN,STAMAX,ENDMIN,ENDMAX
      INTEGER    ROWINF,ROWSUP,CNT
      REAL       SUM

      DO 5  ROW = 1,NY

             IF (ORDSTA(ROW).LT.STAMIN) STAMIN = ORDSTA(ROW)
             IF (ORDSTA(ROW).GT.STAMAX) STAMAX = ORDSTA(ROW)
             IF (ORDEND(ROW).LT.ENDMIN) ENDMIN = ORDEND(ROW)
             IF (ORDEND(ROW).GT.ENDMAX) ENDMAX = ORDEND(ROW)

5     CONTINUE

      DO 10  COL = 1,NX

C --- Find out what is ROWINF and ROWSUP for column COL
C --- This part is easier to program with a DO WHILE. Consider C language.

      IF (COL.GE.STAMAX.AND.COL.LE.ENDMIN) THEN

          ROWINF = 1
          ROWSUP = NY

      ELSE IF (COL.LT.STAMAX) THEN

          ROWINF = 0
          ROWSUP = 0
          DO 20  ROW = 1,NY
             IF (ORDSTA(ROW).LT.COL.AND.ROWINF.EQ.0)  ROWINF = ROW
20        CONTINUE
          DO 25  ROW = NY,1,-1
             IF (ORDSTA(ROW).LT.COL.AND.ROWSUP.EQ.0)  ROWSUP = ROW
25        CONTINUE

      ELSE IF (COL.GT.ENDMIN) THEN

          ROWINF = 0
          ROWSUP = 0
          DO 30  ROW = 1,NY
             IF (ORDEND(ROW).GT.COL.AND.ROWINF.EQ.0)  ROWINF = ROW
30        CONTINUE
          DO 35  ROW = NY,1,-1
             IF (ORDEND(ROW).GT.COL.AND.ROWSUP.EQ.0)  ROWSUP = ROW
35        CONTINUE

      ENDIF

C --- Fill in the buffer with defined values, prolonging the frame by a
C --- mirror effect.

      SUM = 0.
      DO 40 ROW = ROWINF,ROWSUP

          SUM = SUM + INPFRAM(COL,ROW)

40    CONTINUE

      CNT = ROWSUP - ROWINF + 1
      SUM = SUM/CNT

      DO 50 ROW = 1,NY

         OUTFRAM(COL,ROW) = SUM

50    CONTINUE
10    CONTINUE

      RETURN
      END
