C @(#)bintabl.for	17.1.1.1 (ES0-DMD) 01/25/02 17:12:50
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.COPYRIGHT:  Copyright (c) 1987 European Southern Observatory,
C.                               all rights reserved
C.IDENTIFICATION: BINTBL.FOR
C.PURPOSE:    Produces table with binned averages of input table columns
C             together with some statistics. For more info see MIDAS help
C             file.
C.LANGUAGE:   F77+ESOext
C.AUTHOR:     J. Edwin Huizinga
C.ALGORITHM:  Straightforward
C.VERSION:    910109 JEH private version converted to MIDAS standards
C.VERSION:    910624 RHW include in the MIDAS applic area
C.VERSION:    910701 RHW standard fortran code
C             920908 MP correct the type of CT1 and LF1
C---------------------------------------------------------------------
      PROGRAM BINTBL
C
      IMPLICIT NONE
C
      INTEGER       MAXBIN
      INTEGER       MAXROW
      PARAMETER     (MAXBIN=1024)  ! maximum number of bins
      PARAMETER     (MAXROW=16000) ! maximum number of selected table rows
C
      INTEGER       TID
      INTEGER       J
      INTEGER       STATUS,IAV,NC,MADRID(1)
      INTEGER       COLS(5),ROW,I,IND,NVAL,IPOINT(MAXROW)
      INTEGER       NCOL,NROW,NSC,ACOL,AROW,IBIN
      INTEGER       UNIT(1),ICOL(2),TINULL,NBIN
      REAL          SUM(MAXbiN),COUNT(MAXBIN),SD(MAXBIN),sig(MAXBIN)
      REAL          VAL(2),XX(5),BINCEN(MAXBIN)
      REAL          X(MAXROW),Y(MAXROW),MEAN(MAXBIN)
      REAL          TRNULL,USRMIN,USRMAX,USRBSI
      REAL          RMIN,RMAX,DUMMY,SIGMA
      DOUBLE PRECISION TDNULL
      LOGICAL       NULL(2),FLAG
      CHARACTER*1   METHOD,STABIN
      CHARACTER*4   DEFAUL
      CHARACTER*60  TABLE
      CHARACTER*40  COLUM1,COLUM2
      CHARACTER*80  LABEL1,LABEL2,TEXT
      CHARACTER*16  UNIT1,UNIT2,FORM1
      INTEGER       CT1,LF1
C
C
C set up MIDAS:
C
      COMMON        /VMR/MADRID
      INCLUDE       'MID_INCLUDE:TABLES.INC/NOLIST'
      INCLUDE       'MID_INCLUDE:TABLED.INC/NOLIST'
C
   10 FORMAT('*** FATAL: To many selected entries, maximum: ',I6)
C
c *** start midas communications
      CALL STSPRO('BINTBL')
C
C init some values:
C
      DO I=1, MAXBIN
         SD(I)     = 0.0
         SUM(I)    = 0.0
         SIG(I)    = 0.0
         COUNT(I)  = 0
         BINCEN(I) = 0.0
      ENDDO
      CALL TBMNUL(TINULL,TRNULL,TDNULL)
      METHOD = 'A'
      STABIN = 'L'
C
C get table name, and column:
C
      CALL STKRDC('P1',1,1,60,IAV,TABLE,UNIT,NULL,STATUS)
      IF (STATUS.NE.0) THEN
         TEXT = '*** FATAL: Problems with table parameter'//TABLE
         GOTO 9000
      ENDIF
C
      CALL STKRDC('P2',1,1,40,IAV,COLUM1,UNIT,NULL,STATUS)
      CALL STKRDC('P3',1,1,40,IAV,COLUM2,UNIT,NULL,STATUS)
C
C get user binsize,min,max parameters, null if omitted
C
      USRBSI  = TRNULL
      USRMIN  = TRNULL
      USRMAX  = TRNULL
      SIGMA   = TRNULL
      CALL STKRDC('DEFAULT',1,1,4,IAV,DEFAUL,UNIT,NULL,STATUS)
      IF (DEFAUL(1:1).NE.'Y')
     +   CALL STKRDR('INPUTR',1,1,IAV,USRBSI,UNIT,NULL,STATUS)
      IF (DEFAUL(2:2).NE.'Y')
     +   CALL STKRDR('INPUTR',2,1,IAV,USRMIN,UNIT,NULL,STATUS)
      IF (DEFAUL(3:3).NE.'Y')
     +   CALL STKRDR('INPUTR',3,1,IAV,USRMAX,UNIT,NULL,STATUS)
      IF (DEFAUL(4:4).NE.'Y')
     +   CALL STKRDR('INPUTR',4,1,IAV,SIGMA,UNIT,NULL,STATUS)
C
C check user options
C
      IF ((USRMIN.NE.TRNULL).AND.(USRMAX.NE.TRNULL)) THEN
         IF (USRMAX.LT.USRMIN) THEN
            STABIN = 'U'
            DUMMY  = USRMAX
            USRMAX = USRMIN
            USRMIN = DUMMY
         ENDIF
      ENDIF
      IF ((USRBSI.NE.TRNULL).AND.(USRBSI.LT.0.0)) THEN
         METHOD = 'B'
         IBIN = INT(-1.0*USRBSI)
         IF ((USRBSI+IBIN).NE.0.0) THEN
            TEXT = '*** FATAL: Only integer bin_frequencies, dumbo'
            GOTO 9000
         ENDIF
         IF (IBIN.EQ.1) THEN
            TEXT = '*** FATAL: Binsize = 1, does not make sence, dumbo'
            GOTO 9000
         ENDIF
         USRBSI = 1.0*IBIN
      ENDIF
C
C open table, read info and labels, units, formats etc:
C
      CALL TBTOPN(TABLE,F_I_MODE,TID,STATUS)
      IF (STATUS.NE.0) THEN
         TEXT = '*** FATAL: Problems opening table:'//TABLE
         GOTO 9000
      ENDIF
C
      CALL TBIGET(TID,NCOL,NROW,NSC,ACOL,AROW,STATUS)
      IF (STATUS.NE.0) THEN
         TEXT = '*** FATAL: Problems with getting table info'//TABLE
         GOTO 9000
      ENDIF
      IF (NROW.EQ.0) THEN
         TEXT = '*** FATAL: No selected entries, or empty table'
         CALL STTPUT(TEXT,STATUS)
      ENDIF
C
      CALL TBCSER(TID,COLUM1,ICOL(1),STATUS)
      IF (STATUS.NE.0) THEN
         TEXT = '*** FATAL: Problems finding column:'//COLUM1
         GOTO 9000
      ENDIF
C
      CALL TBCSER(TID,COLUM2,ICOL(2),STATUS)
      IF (STATUS.NE.0) THEN
         TEXT = '*** FATAL: Problems finding column:'//COLUM2
         GOTO 9000
      ENDIF
C
      FORM1 = '                '
      CALL TBLGET(TID,ICOL(1),LABEL1,STATUS)
      CALL TBLGET(TID,ICOL(2),LABEL2,STATUS)
      CALL TBUGET(TID,ICOL(1),UNIT1,STATUS)
      CALL TBUGET(TID,ICOL(2),UNIT2,STATUS)
      CALL TBFGET(TID,ICOL(1),FORM1,LF1,CT1,STATUS)
C
C get min,max from table column
C
      CALL TDMXRS(TID,ICOL(1),NROW,0,RMIN,RMAX)
      IF (RMAX.EQ.RMIN) THEN
         TEXT = '*** FATAL: This is useless, column max = min'
         GOTO 9000
      ENDIF
C
C make bins, if no user defaults are given
C
      IF (USRMIN.EQ.TRNULL)  USRMIN = RMIN
      IF (USRMAX.EQ.TRNULL)  USRMAX = RMAX
      IF (USRBSI.EQ.TRNULL)  USRBSI = (USRMAX-USRMIN)/9.99
      NBIN = INT((USRMAX-USRMIN)/USRBSI)+1
      IF ((NBIN.GT.MAXBIN).AND.(METHOD.EQ.'A')) THEN
         TEXT = ' *** FATAL: To many bins, max bins= '
         WRITE(TEXT(37:42),'(I5)') MAXBIN
         GOTO 9000
      ENDIF
C
C doit:
C
      NVAL = 0
      DO ROW=1, NROW
         CALL TBSGET(TID,ROW,FLAG,STATUS)
         IF (FLAG) THEN
            CALL TBRRDR(TID,ROW,2,ICOL,VAL,NULL,STATUS)
            IF (.NOT.(NULL(1).OR.NULL(2))) THEN
               IF ((VAL(1).GE.USRMIN).AND.(VAL(1).LE.USRMAX)) THEN
                  NVAL = NVAL+1
                  IF (NVAL.GT.MAXROW) THEN
                     WRITE(TEXT,10) MAXROW
                     GOTO 9000
                  ENDIF
                  X(NVAL) = VAL(1)
                  Y(NVAL) = VAL(2)
               ENDIF
            ENDIF
         ENDIF
      ENDDO
      CALL TBTCLO(TID,STATUS)
C
C bin the data one way or the other:
C
      NBIN = 0
      IF (METHOD.EQ.'A') THEN
         DO I=1, NVAl
            IND = INT((X(I)-USRMIN)/USRBSI)+1
            IF ((IND.GT.0).AND.(IND.LE.MAXBIN)) THEN
               COUNT(IND) = COUNT(IND) + 1
               SUM(IND)   = SUM(IND) + Y(I)
               IF (IND.GT.NBIN) NBIN = IND
            ELSE
               TEXT= ' *** WARNING: This should never happen'
               CALL STTPUT(TEXT,STATUS)
            ENDIF
         ENDDO
      ELSE                       ! contant number of points per bin
C
C sort x and get rearranged indexing in ipoint
C ipoint can than be used to index y as well
C user indicates starting point by order of min,max
C
         IF (STABIN.EQ.'U') THEN
            CALL TBSRRD(X,IPOINT,NVAL,STATUS)
         ELSE
            CALL TBSRRA(X,IPOINT,NVAL,STATUS)
         ENDIF
         DO I=1, NVAL
            NBIN = INT(1.0*(I-1)/IBIN)+1
            IF (NBIN.GT.MAXBIN) THEN
               TEXT = '*** FATAL: to many bins, maximum = '
               WRITE(TEXT(36:41),'(I6)') MAXBIN
               GOTO 9000
            ENDIF
            SUM(NBIN) = SUM(NBIN)+Y(IPOINT(I))
            BINCEN(NBIN) = BINCEN(NBIN)+X(I)
         ENDDO
         DO I=1, NBIN
            COUNT(I) = IBIN
         ENDDO
         COUNT(NBIN) = MOD(NVAL,IBIN)
         IF (COUNT(NBIN).EQ.0) COUNT(NBIN)=IBIN
      ENDIF
C
C compute mean:
C
      DO I=1, NBIN
         IF (COUNT(I).GT.0) MEAN(I) = SUM(I)/COUNT(I)
      ENDDO
      IF (METHOD.NE.'A') THEN
         DO I=1, NBIN
            IF (COUNT(I).GT.0) BINCEN(I) = BINCEN(I)/COUNT(I)
         ENDDO
      ENDIF
C
C compute standard deviation:
C
      IF (METHOD.EQ.'A') THEN
         DO I=1, NVAL
            IND = INT((X(I)-USRMIN)/USRBSI+1.0)
            IF ((IND.GT.0).AND.(IND.LE.NBIN))
     +         SD(IND) = SD(IND) + (MEAN(IND)-Y(I))**2.0
         ENDDO
      ELSE
         DO I=1, NVAL
            J = INT(1.0*I/IBIN)+1
            SD(J) = SD(J) + (MEAN(J)-Y(IPOINT(I)))**2.0
         ENDDO
      ENDIF
C
C compute sigma of distribution
C
      DO I=1, NBIN
         IF (COUNT(I).GT.1) THEN
            SD(I)=SQRT(SD(I)/(COUNT(I)*(COUNT(I)-1.0)))
            SIG(I)=SD(I)*SQRT(COUNT(I))
         ENDIF
      ENDDO
C
C check for x-sigma deviations from the mean
C
      IF (SIGMA.NE.TRNULL) THEN
         IF (METHOD.EQ.'A') THEN
            DO I=1, NVAL
               IND = INT((X(I)-USRMIN)/USRBSI+1.0)
               IF (ABS(MEAN(IND)-Y(I)).GT.SIGMA*SIG(IND)) THEN
                  IF ((IND.GT.0).AND.(IND.LE.NBIN)) THEN
                     COUNT(IND) = COUNT(IND) - 1
                     SUM(IND)   = SUM(IND) - Y(I)
                  ELSE
                     TEXT= ' *** WARNING: This should never happen'
                     CALL STTPUT(TEXT,STATUS)
                  ENDIF
               ENDIF
            ENDDO
         ELSE                       ! constant number of points per bin
            DO I=1, NBIN
               IF (COUNT(I).GT.1) BINCEN(I) = BINCEN(I)*COUNT(I)
            ENDDO
            DO I=1, NVAL
               NBIN = INT(1.0*(I-1)/IBIN)+1
               IF (ABS(MEAN(NBIN)-Y(IPOINT(I))).GT.SIGMA*SIG(NBIN)) THEN
                  SUM(NBIN) = SUM(NBIN)-Y(IPOINT(I))
                  BINCEN(NBIN) = BINCEN(NBIN)-X(I)
                  COUNT(NBIN) = COUNT(NBIN)-1
               ENDIF
            ENDDO
         ENDIF
         DO i=1, NBIN
            IF (COUNT(I).GT.0) MEAN(I) = SUM(I)/COUNT(I)
         ENDDO
         IF (METHOD.NE.'A') THEN
            DO I=1, NBIN
               IF (COUNT(I).GT.0) BINCEN(I) = BINCEN(I)/COUNT(I)
            ENDDO
         ENDIF
C
C recompute standard deviation:
C
         IF (METHOD.EQ.'A') THEN
            DO I=1, NVAL
               IND = INT((X(I)-USRMIN)/USRBSI+1.0)
               IF ((IND.GT.0).AND.(IND.LE.NBIN))
     +            SD(IND) = SD(IND) + (MEAN(IND)-Y(I))**2.0
            ENDDO
         ELSE
            DO I=1, NVAL
               J = INT(1.0*I/IBIN)+1
               SD(J) = SD(J) + (MEAN(J)-Y(IPOINT(I)))**2.0
            ENDDO
         ENDIF
C
C recompute sigma of distribution
C
         DO I=1, NBIN
            IF (COUNT(I).GT.1) THEN
               SD(I)=SQRT(SD(I)/(COUNT(I)*(COUNT(I)-1.0)))
               SIG(I)=SD(I)*SQRT(COUNT(I))
            ENDIF
         ENDDO
      ENDIF
C
C create output tabel:
C
      CALL TBTINI('bin',F_TRANS,F_O_MODE,5,NBIN,TID,STATUS)
      IF (STATUS.NE.0) THEN
         TEXT = '*** FATAL: Problems creating output table'
         GOTO 9000
      ENDIF
      CALL TBCINI(TID,D_R4_FORMAT,1,FORM1,UNIT1,LABEL1,NC,STATUS)
      CALL TBCINI(TID,D_R4_FORMAT,1,'E15.7',UNIT2,LABEL2,NC,STATUS)
      CALL TBCINI(TID,D_R4_FORMAT,1,'E15.7',UNIT2,'ERR_AV',NC,STATUS)
      CALL TBCINI(TID,D_R4_FORMAT,1,'E15.7',UNIT2,'SIGMA',NC,STATUS)
      CALL TBCINI(TID,D_R4_FORMAT,1,'I8',' ','FREQUENCY',NC,STATUS)
C
C write results
C
      DO I=1,5
         COLS(I) = I
      ENDDO
      DO I=1, NBIN
         IF (METHOD.EQ.'A') THEN
            XX(1) = USRMIN + (I-0.5)*USRBSI
         ELSE
            XX(1) = BINCEN(I)
         ENDIF
         XX(2) = MEAN(I)
         XX(3) = SD(I)
         XX(4) = SIG(I)
         XX(5) = COUNT(I)
         CALL TBRWRR(TID,I,5,COLS,XX,STATUS)
      ENDDO
      CALL TBTCLO(TID,STATUS)
C
c finito:
C
      GOTO 9001
C
C here for errors:
C
 9000 CALL STTPUT(TEXT,STATUS)
 9001 CALL STSEPI
      END
 
 
      SUBROUTINE TBSRRD(A, IPOINT, N, IFAIL)
C
C     SORT REAL ARRAY A(N) IN DESCENDING ORDER
C     POINT(N) CONTAINS THE PREMUTATION SORTING THE ARRAY
C
      IMPLICIT NONE
      INTEGER   N
      INTEGER   IPOINT(N)
      REAL      A(N)
      INTEGER IFAIL
C
      INTEGER   I, II, J, JJ, K, M, IJ, IT, ITT, L
      INTEGER   IU(22),IL(22)
      REAL      T, TT
C
      DO 5 I = 1, N
        IPOINT(I) = I
 5    CONTINUE
C
      II = 1
      JJ = N
      K  = 4
      J  = II - JJ
      M  = 1
      I  = II
      J  = JJ
120   IF (I-J) 140, 400, 400
140   K  = I
      IJ = (J+I)/2
      T  = A(IJ)
      IT = IPOINT(IJ)
      IF (A(I)-T) 160, 180, 180
160   A(IJ) = A(I)
      A(I)  = T
      T     = A(IJ)
      IPOINT(IJ) = IPOINT(I)
      IPOINT(I)  = IT
      IT         = IPOINT(IJ)
180   L = J
      IF (A(J)-T) 260, 260, 200
200   A(IJ) = A(J)
      A(J)  = T
      T     = A(IJ)
      IPOINT(IJ) = IPOINT(J)
      IPOINT(J)  = IT
      IT         = IPOINT(IJ)
      IF (A(I)-T) 220, 260, 260
220   A(IJ) = A(I)
      A(I)  = T
      T     = A(IJ)
      IPOINT(IJ) = IPOINT(I)
      IPOINT(I)  = IT
      IT         = IPOINT(IJ)
      GO TO 260
240   A(L) = A(K)
      A(K) = TT
      IPOINT(L) = IPOINT(K)
      IPOINT(K) = ITT
260   L = L - 1
      IF (A(L)-T) 260, 280, 280
280   TT  = A(L)
      ITT = IPOINT(L)
300   K   = K + 1
      IF (A(K)-T) 320, 320, 300
320   IF (K-L) 240, 240, 340
340   IF (L-I-J+K) 380, 380, 360
360   IL(M) = I
      IU(M) = L
      I = K
      M = M + 1
      GO TO 440
380   IL(M) = K
      IU(M) = J
      J     = L
      M     = M + 1
      GO TO 440
400   M = M - 1
      IF (M) 420, 640, 420
420   I = IL(M)
      J = IU(M)
440   IF (J-I-11) 460, 140, 140
460   IF (I-II) 500, 120, 500
480   I = I + 1
500   IF (I-J) 520, 400, 520
520   T  = A(I+1)
      IT = IPOINT(I+1)
      IF (A(I)-T) 540, 480, 480
540   K           = I
560   A(K+1)      = A(K)
      IPOINT(K+1) = IPOINT(K)
      K           = K - 1
      IF (T-A(K)) 580, 580, 560
580   A(K+1)      = T
      IPOINT(K+1) = IT
      GO TO 480
640   RETURN
      END
 
 
      SUBROUTINE TBSRRA(A, IPOINT, N, IFAIL)
C
C     SORT REAL ARRAY A(N) IN ASCENDING ORDER
C     POINT(N) CONTAINS THE PREMUTATION SORTING THE ARRAY
C
      IMPLICIT NONE
      INTEGER   N
      INTEGER   IPOINT(N)
      REAL      A(N)
      INTEGER IFAIL
C
      INTEGER   I, II, J, JJ, K, M, IJ, IT, ITT, L
      INTEGER   IU(22),IL(22)
      REAL      T, TT

      DO 5 I = 1, N
        IPOINT(I) = I
5     CONTINUE
C
      II = 1
      JJ = N
      K  = 4
      J  = II - JJ
      M  = 1
      I  = II
      J  = JJ
120   IF (I-J) 140, 400, 400
140   K  = I
      IJ = (J+I)/2
      T  = A(IJ)
      IT = IPOINT(IJ)
      IF (T-A(I)) 160, 180, 180
160   A(IJ) = A(I)
      A(I)  = T
      T     = A(IJ)
      IPOINT(IJ) = IPOINT(I)
      IPOINT(I)  = IT
      IT         = IPOINT(IJ)
180   L = J
      IF (T-A(J)) 260, 260, 200
200   A(IJ) = A(J)
      A(J)  = T
      T     = A(IJ)
      IPOINT(IJ) = IPOINT(J)
      IPOINT(J)  = IT
      IT         = IPOINT(IJ)
      IF (T-A(I)) 220, 260, 260
220   A(IJ) = A(I)
      A(I)  = T
      T     = A(IJ)
      IPOINT(IJ) = IPOINT(I)
      IPOINT(I)  = IT
      IT         = IPOINT(IJ)
      GO TO 260
240   A(L) = A(K)
      A(K) = TT
      IPOINT(L) = IPOINT(K)
      IPOINT(K) = ITT
260   L = L - 1
      IF (T-A(L)) 260, 280, 280
280   TT  = A(L)
      ITT = IPOINT(L)
300   K   = K + 1
      IF (T-A(K)) 320, 320, 300
320   IF (K-L) 240, 240, 340
340   IF (L-I-J+K) 380, 380, 360
360   IL(M) = I
      IU(M) = L
      I = K
      M = M + 1
      GO TO 440
380   IL(M) = K
      IU(M) = J
      J     = L
      M     = M + 1
      GO TO 440
400   M = M - 1
      IF (M) 420, 640, 420
420   I = IL(M)
      J = IU(M)
440   IF (J-I-11) 460, 140, 140
460   IF (I-II) 500, 120, 500
480   I = I + 1
500   IF (I-J) 520, 400, 520
520   T  = A(I+1)
      IT = IPOINT(I+1)
      IF (T-A(I)) 540, 480, 480
540   K           = I
560   A(K+1)      = A(K)
      IPOINT(K+1) = IPOINT(K)
      K           = K - 1
      IF (A(K)-T) 580, 580, 560
580   A(K+1)      = T
      IPOINT(K+1) = IT
      GO TO 480
640   RETURN
      END