C @(#)mode.for	19.1 (ES0-DMD) 02/25/03 13:30:47
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+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.IDENTIFICATION
C  subroutine MODE         version 3         820518
C  A. Kruszewski           ESO - Garching
C.PURPOSE
C  calculates mode of data contained in array "A"
C.ALGORITHM
C  mode is calculated with use of relation between mean, median
C  and mode, see:  Kendall, M., Stuart, A. 1977, The Advanced Theory
C  of Statistics, 4-th edition, p. 40
C  as this relation is valid only for mildly assymetric distributions
C  the data is repeatedly clipped around previous value of mode
C  first truncation is made around starting value of mean
C  if the dimension "N" of array "A" is greater than 5000 then only
C  every II-th element is used where II=N/5000+1
C
C  control array "CRMD" holds three parameters:
C  "CRMD(1)" introduces possibility of choice between mode and mean
C  subroutine returns mode when "CRMD(1)=3.0"
C  subroutine returns clipped mean when "CRMD(1)=0.0"
C  other values of "CRMD(1)" can be used when necessary
C  "CRMD(3)" is a clipping factor - data is clipped by "SIGMA*CRMD(3)"
C  from both sides of recent value of "AMODE"
C  iterations are terminated when the difference between two last
C  values of mode divided by last value of sigma is smaller than
C  "CRMD(2)"
C  maximum number of iterations is 10
C.INPUT/OUTPUT
C  input arguments:
C  A           real array        data array
C  N           integer           number of values in "A"
C  CRMD        real array        values of control parameters
C  output arguments:
C  AMODE       real              mode of array "A"
C  SIGMA       real              sigma of "A" after last clipping
C-----------------------------------------------------------------------
      SUBROUTINE MODE(A,N,CRMD,AMODE,SIGM)
C
      IMPLICIT NONE
      INTEGER  N
      REAL     A(N)
      REAL     CRMD(3)
      REAL     SIGM
C
      INTEGER  I, IC, II, IN, INN, L
      REAL     ALC, AHC, AMEAN, AMED, AMODE
      REAL     B(5000), C(5000), CFTR
      REAL     DFTR, SCR, SIGMA, SFTR
C
C
      DFTR = CRMD(1)
      CFTR = CRMD(2)
      SFTR = CRMD(3)
C
C *** Define array B of dimension IN.
      IF (N.GT.5000) THEN
         II = (N-1) / 5000 + 1
         IN = N / II
         INN = II * IN
         DO 10 I = II , INN , II
            B(I/II) = A(I)
   10    CONTINUE
      ELSE
         DO 11 I = 1 , N
            B(I) = A(I)
   11    CONTINUE  
         IN = N
      ENDIF
C
C *** Calculate mean of non truncated data and
C *** set it as a first approximation for mode.
      CALL MEAN(B, IN, AMODE, SIGMA)
      SIGM = SIGMA
C
C *** Set iteration control.
      IC = 9
C
C *** Perform iterations.
   12 CONTINUE 
C
C *** Save recent value of mode.
      SCR = AMODE
C
C *** Set limits of truncation.
      ALC = AMODE - SFTR * SIGMA
      AHC = AMODE + SFTR * SIGMA
C
C *** Continue iterations.
      L=0
      DO 13 I = 1 , IN
         IF ( B(I) .GT. ALC .AND. B(I) .LT. AHC ) THEN
            L = L + 1
            C(L) = B(I)
         ENDIF
   13 CONTINUE 
C
C *** Stops iterations if number of data points
C *** after recent truncation is less than 2.
      IF ( L .LT. 2 ) THEN
         IC = 0
      ELSE
         CALL MEAN(C, L, AMEAN, SIGMA)
C
C *** Replaces mode with mean when DFTR.LE.0.1.
         IF ( DFTR .GT. 0.1 ) THEN
            CALL MEDIAN(C, L, ALC, AHC, AMED)
            IF ( L .GE. 10 ) THEN
               AMODE = DFTR * AMED - (DFTR-1.0) * AMEAN
            ELSE
               AMODE = AMED
            ENDIF
         ELSE
            AMODE=AMEAN
         ENDIF
C   
C *** Compare two most recent values of mode and
C *** stop iterations if they differ by less than
C *** expected mean error of the mean of truncated
C *** distribution multiplied by factor 'SGFACTOR'.
         SCR = MAX( ABS(SCR-AMODE) , ABS(AMODE-AMEAN) )
C
C *** Stops iterations if sigma is not positive.
         IF (SIGMA.LE.0.0) THEN
            IC=0
         ELSE
            SCR = SCR / ( SQRT( (SIGMA*SIGMA)/L ) )
            IF ( SCR .LT. CFTR ) THEN
               IC = 0
            ELSE
               IC = IC - 1
            ENDIF
         ENDIF
      ENDIF
      IF (IC.GT.0) GOTO 12
C
      RETURN
      END



