C @(#)mode.for	17.1.1.1 (ES0-DMD) 01/25/02 17:15:43
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
C  A           real*4 array        data array
C  N           integer*4           number of values in "A"
C  CRMD        real*4 array        values of control parameters
C  output arguments
C  AMODE       real*4              mode of array "A"
C  SIGMA       real*4              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.
C
      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.
C
      CALL MEAN( B , IN , AMODE , SIGMA )
      SIGM = SIGMA
C
C ******      Set iteration control.
C
      IC = 9
C
C ******      Perform iterations.
C
   12      CONTINUE 
C
C ******          Save recent value of mode.
C
            SCR = AMODE
C
C ******          Set limits of truncation.
C
          ALC = AMODE - SFTR * SIGMA
          AHC = AMODE + SFTR * SIGMA
C
C ******          Continue iterations.
C
          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.
C
          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.
C
            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'.
C   
            SCR = MAX( ABS(SCR-AMODE) , ABS(AMODE-AMEAN) )
C
C ******            Stops iterations if sigma is not positive.
C
            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
C
      END