C===========================================================================
C Copyright (C) 1995,2004 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 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence 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,2004 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  16:48 - 19 NOV 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION
C
C  program TDTSTAT.FOR
C
C.PURPOSE
C
C  Execute the command
C  STAT/TABLE table column-ref [step [min-val  [max-val]]] LOG10
C
C.KEYWORDS
C
C  histogram, table
C
C.ALGORITHM
C
C  use table interface routines
C
C
C.VERSION
C 04021	last modif
C
C-----------------------------------------------------------
C
      SUBROUTINE TDSTAT
      IMPLICIT NONE
C
      LOGICAL TRIV
C
      INTEGER IOUT(2), ISTAT
      INTEGER PARVAL, NPAR, TID, NROW, NCOL, NSC, NAC, NAR
      INTEGER IFREQ(512)
      INTEGER I, II, II1, II2, DUMMY, NVALS, ICOL
      INTEGER IAC, NACT, INS, IX, NNN, NBINS, IDUMMY, LL
      INTEGER DUNIT
      INTEGER INDEX
      INTEGER TINULL
      INTEGER KUN, KNUL, DTYPE
C
      REAL TRNULL
      REAL TSTAT(8),XMEAN,XSTD
      REAL CR,STEP,X0,X1,XT,AXMAX,AXMIN
      REAL CINT(512),XMIN,XMAX
      REAL ROUT(4),RPAR(4)
C
      DOUBLE PRECISION TDNULL, DXMIN, DXMAX
C
      CHARACTER*16 LABEL
      CHARACTER*16 MSG
      CHARACTER*80 TABLE, LINE
      CHARACTER*17 COLREF
      CHARACTER*80 WORK
      CHARACTER*10 DESNA1, DESNA2, DESNA3
      CHARACTER*4  WS,COMLIN
      CHARACTER*8  TYPE,FORM
C
      INCLUDE 'MID_INCLUDE:TABLES.INC'
      INCLUDE 'MID_INCLUDE:TABLED.INC'
      DATA PARVAL/7/
      DATA DESNA1/'TSTATxxx  '/
      DATA DESNA2/'TCLASxxx  '/
      DATA DESNA3/'TFREQxxx  '/
      DATA MSG/'ERR:TSTATTBLxxxx'/
C
C ... GET MACHINE CONSTANTS
C
      CALL TBMNUL(TINULL, TRNULL, TDNULL)
C
C ... get input parameters + default
C
      CALL TDPGET(PARVAL,NPAR,ISTAT)
      IF (ISTAT.NE.0) GO TO 3000
      TABLE   = TPARBF(1)
      COLREF  = TPARBF(2)
      RPAR(2) = TRNULL
      RPAR(3) = RPAR(2)
      RPAR(1) = 0.0
C 
      IF (NPAR .GE. 3) THEN
         II1 = 1
         DO 10, II2 = 3,NPAR
            II     = INDEX(TPARBF(II2),' ') - 1
            WORK   = TPARBF(II2) (1:II)
            IF (WORK(1:1).NE.'L' .AND. WORK(1:1).NE.'?') THEN
               CALL GENCNV(WORK,2,1,IDUMMY,RPAR(II1),DUMMY,NVALS)
               II1    = II1 + 1
               IF (NVALS.NE.1) THEN
                  CALL STTPUT('Invalid input param. ...',ISTAT)
                  GO TO 3000
               ENDIF
            ENDIF
10       CONTINUE
      ENDIF
C 
      STEP   = RPAR(1)
      AXMIN  = RPAR(2)
      AXMAX  = RPAR(3)
C
C ... init input table
C
      TID = -1
      CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT)
      CALL TBIGET(TID,NCOL,NROW,NSC,NAC,NAR,ISTAT)
      IF (NROW.EQ.0) THEN
          CALL STTPUT(' Table empty ...',ISTAT)
          GO TO 3000
      ENDIF
C
C ... find input column
C
      CALL TBCSER(TID,COLREF,ICOL,ISTAT)
      IF (ICOL.EQ.-1) THEN
          CALL STTPUT('Column not found...',ISTAT)
          GO TO 3000
      ENDIF
      CALL TBFGET(TID,ICOL,FORM,LL,DTYPE,ISTAT)
      IF (DTYPE.EQ.D_C_FORMAT) THEN
          CALL STTPUT('Error, Column of type Character',ISTAT)
          GO TO 3000
      ENDIF
C
C ... build the grid
C
      IF (DTYPE.EQ.D_R8_FORMAT) THEN
          CALL TDMXDS(TID,ICOL,NROW,0,DXMIN,DXMAX)
          TYPE = 'R*8'
          XMIN = DXMIN
          XMAX = DXMAX
      ELSE IF (DTYPE.EQ.D_R4_FORMAT) THEN
          CALL TDMXRS(TID,ICOL,NROW,0,XMIN,XMAX)
          TYPE = 'R*4'
      ELSE
          CALL TDMXIS(TID,ICOL,NROW,0,XMIN,XMAX)
          IF (DTYPE.EQ.D_I1_FORMAT) THEN
              TYPE = 'I*1'
          ELSE IF (DTYPE.EQ.D_I2_FORMAT) THEN
              TYPE = 'I*2'
          ELSE
              TYPE = 'I*4'
          END IF
      ENDIF
C 
      IF (XMAX.LE.(XMIN+STEP)) THEN
          XMAX = XMIN + STEP
          IF (XMIN.GE.XMAX) THEN
             CALL STTPUT
     2       ('WARNING: XMAX = XMIN; increased XMAX by 1.0',ISTAT)
             XMAX = XMIN + 1
          ENDIF
          TRIV = .TRUE.
      ELSE
          TRIV = .FALSE.
      ENDIF
C 
      CALL TDSCAL(XMIN,XMAX,1.,X0,X1,IX,XT,NNN)
      IF (NNN .LE. 0) THEN                      !no dynamic range...
         XMEAN = XMIN
         XSTD = 0
         NBINS = 1
         NACT = 1
         GOTO 1000
      ENDIF
C 
      IF (AXMIN.EQ.TRNULL) THEN
          AXMIN  = X0*10.**IX
      END IF
      IF (AXMAX.EQ.TRNULL) THEN
          AXMAX  = X1*10.**IX
      END IF
      IF (STEP.LE.0.0) STEP = XT
C 
      NBINS  = (AXMAX-AXMIN)/STEP + 2
      NBINS  = MIN(NBINS,512)
      CR     = AXMIN
      CINT(1) = CR
      DO 20, I = 1,NBINS - 2
         CR = CR + STEP
         CINT(I+1) = CR
20    CONTINUE
C
C ... frequency analysis
C
      IF (DTYPE.EQ.D_R8_FORMAT) THEN
          CALL TDDSTA(TID,ICOL,NROW,NBINS,1,CINT,
     +                IFREQ,NACT,XMEAN,XSTD)
      ELSE IF (DTYPE.EQ.D_R4_FORMAT) THEN
          CALL TDRSTA(TID,ICOL,NROW,NBINS,1,CINT,
     +                IFREQ,NACT,XMEAN,XSTD)
      ELSE
          CALL TDISTA(TID,ICOL,NROW,NBINS,1,CINT,
     +                IFREQ,NACT,XMEAN,XSTD)
      ENDIF
      IF (TRIV) XMAX = XMIN
C 
      IF (NACT.EQ.0) 
     +   CALL STTPUT(' No selected entries or empty column',ISTAT)
C
C ... write results as descriptors
C
1000  INS = 1000 + ICOL
      WRITE (WS,9000) INS
      DESNA1(6:8) = WS(2:4)
      DESNA2(6:8) = WS(2:4)
      DESNA3(6:8) = WS(2:4)
      CALL STKRDC('MID$CMND',1,1,4,IAC,COMLIN,KUN,KNUL,ISTAT)
C
C ... STATISTIC/TABLE
C
      IF (COMLIN(1:4).EQ.'STAT') THEN
          LINE =  'Table : '//TABLE
          CALL STTPUT(LINE,ISTAT)
          WRITE (WORK,9000) ICOL
          CALL TBLGET(TID,ICOL,LABEL,ISTAT)
          LINE = 'Column #'//WORK(1:4)//' Label :'//LABEL//
     +           ' Type :'//TYPE
          CALL STTPUT(LINE, ISTAT)
          WRITE (WORK,9020) NROW,NACT
          CALL STTPUT(WORK,ISTAT)
          WRITE (WORK,9030) XMIN,XMAX
          CALL STTPUT(WORK,ISTAT)
          WRITE (WORK,9040) XMEAN,XSTD
          CALL STTPUT(WORK,ISTAT)
      ELSE
          TSTAT(1) = 1.
          TSTAT(2) = NACT
          TSTAT(3) = XMEAN
          TSTAT(4) = XSTD
          TSTAT(5) = XMIN
          TSTAT(6) = XMAX
          TSTAT(7) = NBINS
          TSTAT(8) = 0.
          CALL STDWRR(TID,DESNA1,TSTAT,1,8,DUNIT,ISTAT)
          IF (ISTAT.NE.0) GO TO 3000
          CALL STDWRR(TID,DESNA2,CINT,1,NBINS,DUNIT,ISTAT)
          IF (ISTAT.NE.0) GO TO 3000
          CALL STDWRI(TID,DESNA3,IFREQ,1,NBINS,DUNIT,ISTAT)
          IF (ISTAT.NE.0) GO TO 3000
      ENDIF
      IOUT(1) = ICOL
      IOUT(2) = NACT
      ROUT(1) = XMIN
      ROUT(2) = XMAX
      ROUT(3) = XMEAN
      ROUT(4) = XSTD
      CALL STKWRI('OUTPUTI',IOUT,1,2,KUN,ISTAT)
      CALL STKWRR('OUTPUTR',ROUT,1,4,KUN,ISTAT)
C
C ... end
C
      CALL TBTCLO(TID,ISTAT)
3000  IF (ISTAT.NE.0) THEN
         WRITE (MSG(13:16),9000) ISTAT
         CALL TDERRR(ISTAT,MSG,I)
      ENDIF
      RETURN
 9000 FORMAT (I4)
 9020 FORMAT ('Total no. of entries :',I8,', selected no. of entries :',
     +       I8)
C 9030 FORMAT ('Minimum value :',E14.5,', Maximum value:',E14.5)
C 9040 FORMAT ('Mean value    :',E14.5,', Standard dev.:',E14.5)
 9030 FORMAT ('Minimum value :',g16.7,', Maximum value:',g16.7)
 9040 FORMAT ('Mean value    :',g16.8,', Standard dev.:',g16.8)
      END
