C @(#)avarea.for	17.1.1.1 (ES0-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 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
      REAL FUNCTION AVAREA(METHOD,A,NPIX,SUBPIX,NDIM,O)
C
C ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C	
C.IDENTIFICATION
C  real function AVAREA		version 2.80	880630
C  K. Banse			ESO - Garching
C 
C.KEYWORDS
C  median filter, average
C 
C.PURPOSE
C  find median value or average
C 
C.ALGORITHM
C  as in [MIDAS.PRIM.SUB1]AVERAGE.FOR
C	
C.INPUT/OUTPUT
C  use as VALUE = AV_AREA(METHOD,A,NPIX,SUBPIX,NDIM,O)
C  
C input par:
C  METHOD:	char. exp.	A for average, M for median
C  A:		R*4 array	input array
C  NPIX:	I*4		no. of pixels per line
C  SUBPIX:	I*4 array	start pixels of area
C  NDIM:	I*4 array	no. of pixels in subarea
C  O:		R*4 array	working buffer
C	
C  AVAREA:	R*4		median value to be returned
C	
C --------------------------------------------------------------------
C
      IMPLICIT NONE
C 
      INTEGER      NPIX,NDIM(*),SUBPIX(*)
      INTEGER      INDX,NH,NN,NNN,NTOT,NX,NY,OFF,TOP
C 
      REAL         A(*),O(*)
C 
      CHARACTER*(*)  METHOD
C 
      DOUBLE PRECISION   SUM
C      
C  point to first pixel (in lower left corner)
      OFF = (SUBPIX(2)-1)*NPIX + SUBPIX(1)
      NTOT = NDIM(1) * NDIM(2)
C      
C  branch according to method
      IF (METHOD(1:1).EQ.'A') GOTO 5000
C      
C  here for median finding
C      
      O(1) = A(OFF)
      TOP = 1
      NH = (NTOT + 1)/2             !index of median
C
C  move through rest of first line
      DO 1200 NX=1,NDIM(1)-1
         INDX = OFF + NX                  !point to "first" pixel in row
C      
C  fill array O only from 1 to NH
         DO 800 NN=1,TOP                  !loop through already ordered set
            IF (A(INDX).LT.O(NN)) THEN
               DO 600 NNN=TOP,NN,-1
                  O(NNN+1) = O(NNN)       !shift up ...
600            CONTINUE
             O(NN) = A(INDX)              !merge value into array O
             GOTO 1000
            ENDIF
800      CONTINUE
C 
         O(TOP+1) = A(INDX)               !add value on top
1000     IF (TOP.LT.NH) TOP = TOP + 1
C      
1200  CONTINUE
C      
C  test, if we are only handling single line
      IF (NDIM(2).LE.1) THEN
         AVAREA = O(NH)
         RETURN
      ELSE
         OFF = OFF + NPIX - 1             !point to next line (start value - 1)
      ENDIF
C      
C  handle remaining lines
      DO 2300 NY=1,NDIM(2)-1
         DO 2200 NX=1,NDIM(1)
            INDX = OFF + NX               !point to "first" pixel in row
C      
C  fill array O only from 1 to NH
            DO 1800 NN=1,TOP              !loop through already ordered set
               IF (A(INDX).LT.O(NN)) THEN
                  DO 1600 NNN=TOP,NN,-1
                     O(NNN+1) = O(NNN)    !shift up ...
1600              CONTINUE
                  O(NN) = A(INDX)         !merge value into array O
                  GOTO 2000
               ENDIF
1800        CONTINUE
C 
            O(TOP+1) = A(INDX)            !add value on top
2000        IF (TOP.LT.NH) TOP = TOP + 1
2200     CONTINUE
         OFF = OFF + NPIX
2300  CONTINUE
C      
C  return median value
      AVAREA = O(NH)
      RETURN
C      
C  here for simple averaging
C      
5000  OFF = OFF - 1
      SUM = 0.D0
      DO 6000 NY=1,NDIM(2)
         DO 5500 NX=1,NDIM(1)
            SUM = SUM + A(OFF+NX)
5500     CONTINUE
         OFF = OFF + NPIX
6000  CONTINUE
C      
C  return average value
      AVAREA = SNGL(SUM/DBLE(NTOT))
      RETURN
C      
      END