C @(#)tdstat.for	17.1.1.1 (ES0-DMD) 01/25/02 17:47:17
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
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  19:23 - 11 DEC 1987
C.VERSION: 1.1  NO MAPPING VERSION - 11 DEC 1989
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: J.D.PONZ
C
C.IDENTIFICATION        TDSTAT.FOR
C.KEYWORDS           TABLE, APPLICATIONS
C.ENVIRONMENT  MIDAS
C
C.COMMENTS 
C    SOURCE HAS TO BE MODIFIED TO USE THE TABLE TURBO VERSION
C MODIF : M Peron 30JAN91 sigma = SQRT(sum(Xi-MEAN)/(N-1))
C------------------------------------------------------------------

      SUBROUTINE TDISTA
     . (TID,ICOL,N,K2,ICLASS,CINT,IFREQ,NACT,XMEAN,XSTD)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C COMPUTE FREQUENCY DISTRIBUTION FOR VALUES IN THE ARRAY.
C OTHER STATISTICAL PARAMETERS RETURNED BY THE ROUTINE ARE
C MEAN, STD.DEV., MINIMUM AND MAXIMUM.
C       INTEGER VERSION
C
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER TID        ! IN  : TABLE IDENTIFIER
      INTEGER ICOL       ! IN  : COLUMN INDEX
      INTEGER N          ! IN  : NUMBER OF INPUT DATA
      INTEGER K2         ! IN  : NUMBER OF CLASES
      INTEGER ICLASS     ! IN  :
      REAL CINT(K2)      ! IN  :
      INTEGER IFREQ(K2)  ! IN  :
      INTEGER NACT       ! OUT : ACTUAL NUMBER OF POINTS
      REAL XMEAN         ! OUT : MEAN VALUE
      REAL XSTD          ! OUT : STANDARD DEVIATION
C
      LOGICAL LK0, INULL, ISEL
      INTEGER I,NUP,NDN,K0
      INTEGER JJ,IERR,NOB,LDN
      INTEGER LUP,J, IVAL
      INTEGER STATUS
      REAL    XMIN,XMAX,XTEMP,CR,ALPHA
      REAL    STEP
      DOUBLE PRECISION YSTD,YMEAN,YWORK
    
C
      NACT   = 0
      YSTD   = 0.00D+0
      YMEAN  = 0.00D+0
      IF (K2-2) 290,10,10
   10 IF (N-1) 300,20,20
C
C ... ZERO  FREQUENCIES
C
   20 CONTINUE
      DO 30 I = 1,K2
          IFREQ(I) = 0
   30 CONTINUE
      NUP    = 0
      NDN    = 0
      K0     = K2 - 2
      LK0    = K0 .EQ. 0
C
C ... DETERMINE  MIN  AND  MAX, MEAN AND STD DEV.
C
      I      = 1
      CALL TBSGET(TID,I,ISEL,STATUS)
      CALL TBERDI(TID,I,ICOL,IVAL,INULL,STATUS)
      XTEMP = IVAL
   40     CONTINUE
          IF (.NOT.INULL .AND. ISEL) GO TO 50
          I      = I + 1
          IF (I.GT.N) GOTO 65
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDI(TID,I,ICOL,IVAL,INULL,STATUS)
          XTEMP = IVAL
          GO TO 40 
   50 CONTINUE
      XMIN   = XTEMP
      XMAX   = XTEMP
      DO 60 J = I,N
          CALL TBSGET(TID,J,ISEL,STATUS)
          CALL TBERDI(TID,J,ICOL,IVAL,INULL,STATUS)
          XTEMP = IVAL
          IF (.NOT.INULL .AND. ISEL) THEN
              NACT   = NACT + 1
              YMEAN  = YMEAN + XTEMP
              IF (XTEMP.LT.XMIN) XMIN   = XTEMP
              IF (XTEMP.GT.XMAX) XMAX   = XTEMP
          END IF
   60 CONTINUE
   65 CONTINUE
      IF (NACT.EQ.0) THEN
          XMEAN  = 0.
          XSTD   = 0.
          RETURN
      ELSE
          YMEAN  = YMEAN/DBLE(NACT)
      END IF
      IF (ICLASS.EQ.1) GO TO 120
      IF (LK0) GO TO 260
C
C ... CALCULATE CLASS INTERVALS WITH EQUAL SPACING
C
      CR     = XMAX - XMIN
      ALPHA  = 0.001
      IF (CR.GT.0.0) GO TO 100
      IF (XMIN) 80,90,80
   80 CR     = XMIN*ALPHA
      GO TO 100
   90 CR     = ALPHA
  100 STEP   = CR* (1.00+ALPHA)/FLOAT(K0)
      CR     = XMIN - 0.50*ALPHA*CR
      CINT(1) = CR
      DO 110 I = 1,K0
          XTEMP  = CR
          CR     = XTEMP + STEP
          CINT(I+1) = CR
  110 CONTINUE
      GO TO 150
C
C ... CLASS INTERVALS SUPPLIED--CHECK IN ASCENDING ORDER
C
  120 CR     = CINT(1)
      IF (LK0) GO TO 270
      JJ     = K0 + 1
      DO 140 I = 2,JJ
          XTEMP  = CINT(I)
          IF (CR.LT.XTEMP) GO TO 130
          IERR   = 3
          GO TO 310
  130     CR     = XTEMP
  140 CONTINUE
C
C ... DETERMINE CLASS TO WHICH EACH CASE BELONGS
  150 NOB    = K2/2
      CR     = CINT(NOB)
      LDN    = NOB + 1
      LUP    = K2 - 1
      K0     = NOB - 1
      LK0    = K0 .EQ. 0
      DO 240 I = 1,N
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDI(TID,I,ICOL,IVAL,INULL,STATUS)
          XTEMP = IVAL
          IF (INULL .OR. (.NOT.ISEL)) GO TO 230
          YWORK  = XTEMP
          YWORK  = YWORK - YMEAN
          YSTD   = YSTD + YWORK*YWORK
          IF (XTEMP-CR) 190,160,160
C
C ... VARIATE VALUE ABOVE MIDDLE  INTERVAL
C
  160     DO 180 J = LDN,LUP
              IF (XTEMP.GE.CINT(J)) GO TO 170
              IFREQ(J) = IFREQ(J) + 1
              GO TO 230
  170         CONTINUE
  180     CONTINUE
          NUP    = NUP + 1
          GO TO 230
C
C ... VARIATE VALUE BELOW MIDDLE  INTERVAL
C
  190     IF (LK0) GO TO 220
          DO 210 J = 1,K0
              JJ     = NOB - J
              IF (XTEMP.LT.CINT(JJ)) GO TO 200
              JJ     = JJ + 1
              IFREQ(JJ) = IFREQ(JJ) + 1
              GO TO 230
  200         CONTINUE
  210     CONTINUE
  220     NDN    = NDN + 1
  230     CONTINUE
  240 CONTINUE
  250 IFREQ(1) = NDN
      IFREQ(K2) = NUP
      XMEAN  = YMEAN
      IF (NACT.EQ.1) THEN
         YSTD = 0
      ELSE
         YSTD   = DSQRT(YSTD/(NACT-1))
      ENDIF
      XSTD   = YSTD
      RETURN
  260 CR     = (XMAX+XMIN)*0.50
      CINT(1) = CR
  270 CONTINUE
      DO 280 I = 1,N
          CALL TBERDI(TID,I,ICOL,IVAL,INULL,STATUS)
          IF (IVAL.LT.CR) NDN    = NDN + 1
  280 CONTINUE
      NUP    = N - NDN
      GO TO 250
  290 IERR   = 1
      GO TO 310
  300 IERR   = 2
  310 NACT   = 0
      RETURN
      END

      SUBROUTINE TDRSTA
     . (TID,ICOL,N,K2,ICLASS,CINT,IFREQ,NACT,XMEAN,XSTD)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C COMPUTE FREQUENCY DISTRIBUTION FOR VALUES IN THE ARRAY.
C OTHER STATISTICAL PARAMETERS RETURNED BY THE ROUTINE ARE
C MEAN, STD.DEV., MINIMUM AND MAXIMUM.
C       SINGLE PRECISION VERSION
C
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER TID        ! IN  : TABLE IDENTIFIER
      INTEGER ICOL       ! IN  : COLUMN INDEX
      INTEGER N          ! IN  : NUMBER OF INPUT DATA
      INTEGER K2         ! IN  : NUMBER OF CLASES
      INTEGER ICLASS     ! IN  :
      REAL CINT(K2)      ! IN  :
      INTEGER IFREQ(K2)  ! IN  :
      INTEGER NACT       ! OUT : ACTUAL NUMBER OF POINTS
      REAL XMEAN         ! OUT : MEAN VALUE
      REAL XSTD          ! OUT : STANDARD DEVIATION
C
      LOGICAL LK0
      LOGICAL INULL
      LOGICAL ISEL
C      
      INTEGER I,NUP,NDN,K0
      INTEGER JJ,IERR,NOB,LDN
      INTEGER LUP,J
      INTEGER STATUS
C
      REAL    XMIN,XMAX,XTEMP
      REAL    CR         ! ???
      REAL    ALPHA
      REAL    STEP
      DOUBLE PRECISION YSTD,YMEAN,YWORK
      INTEGER ISTAT    
C
C     set to zero
      NACT   = 0
      YSTD   = 0.00D+0
      YMEAN  = 0.00D+0
      IF (K2-2) 290,10,10
   10 IF (N-1) 300,20,20
C
C ... ZERO  FREQUENCIES
C
   20 CONTINUE
      DO 30 I = 1,K2
          IFREQ(I) = 0
   30 CONTINUE
      NUP    = 0
      NDN    = 0
      K0     = K2 - 2
      LK0    = K0 .EQ. 0
C
C ... DETERMINE  MIN  AND  MAX, MEAN AND STD DEV.
C
C-     CALL STTPUT('START WITH THE CALCULATION',ISTAT)
      I      = 1
      CALL TBSGET(TID,I,ISEL,STATUS)
      CALL TBERDR(TID,I,ICOL,XTEMP,INULL,STATUS)
   40     CONTINUE
          IF (.NOT.INULL .AND. ISEL) GO TO 50
          I      = I + 1
          IF (I.GT.N) GOTO 65
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDR(TID,I,ICOL,XTEMP,INULL,STATUS)
          GO TO 40 
   50 CONTINUE
      XMIN   = XTEMP
      XMAX   = XTEMP
      DO 60 J = I,N
          CALL TBSGET(TID,J,ISEL,STATUS)
          CALL TBERDR(TID,J,ICOL,XTEMP,INULL,STATUS)
          IF (.NOT.INULL .AND. ISEL) THEN
              NACT   = NACT + 1
C              CALL STTPUT('Increase NACT by one',ISTAT)
              YMEAN  = YMEAN + XTEMP
              IF (XTEMP.LT.XMIN) XMIN   = XTEMP
              IF (XTEMP.GT.XMAX) XMAX   = XTEMP
          END IF
   60 CONTINUE
   65 CONTINUE
      IF (NACT.EQ.0) THEN
          CALL STTPUT
     +    ('No selected rows, set mean and standard dev. to 0.0',ISTAT)
          XMEAN  = 0.
          XSTD   = 0.
          RETURN
      ELSE
          YMEAN  = YMEAN/DBLE(NACT)
      END IF
      IF (ICLASS.EQ.1) GO TO 120
      IF (LK0) GO TO 260
C
C ... CALCULATE CLASS INTERVALS WITH EQUAL SPACING
C
C      write(6,*) XMAX,XMIN
      CR     = XMAX - XMIN
      ALPHA  = 0.001
      IF (CR.GT.0.0) GO TO 100
      IF (XMIN) 80,90,80
   80 CR     = XMIN*ALPHA
      GO TO 100
   90 CR     = ALPHA
  100 STEP   = CR* (1.00+ALPHA)/FLOAT(K0)
      CR     = XMIN - 0.50*ALPHA*CR
      CINT(1) = CR
      DO 110 I = 1,K0
          XTEMP  = CR
          CR     = XTEMP + STEP
          CINT(I+1) = CR
  110 CONTINUE
      GO TO 150
C
C ... CLASS INTERVALS SUPPLIED--CHECK IN ASCENDING ORDER
C
  120 CR     = CINT(1)
      IF (LK0) GO TO 270
      JJ     = K0 + 1
C      write(6,*) JJ,CINT(1),CINT(2)
      DO 140 I = 2,JJ
          XTEMP  = CINT(I)
C
C Typo origin of error for similar column entries? PN 12/98
C          IF (CR.LT.XTEMP) GO TO 130
C  LE. instead of LT.?
          IF (CR.LE.XTEMP) GO TO 130
          CALL STTPUT('Warning: CR larger XTEMP!',ISTAT)
          IERR   = 3
          GO TO 310
  130     CR     = XTEMP
  140 CONTINUE
C
C ... DETERMINE CLASS TO WHICH EACH CASE BELONGS
  150 NOB    = K2/2
      CR     = CINT(NOB)
      LDN    = NOB + 1
      LUP    = K2 - 1
      K0     = NOB - 1
      LK0    = K0 .EQ. 0
      DO 240 I = 1,N
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDR(TID,I,ICOL,XTEMP,INULL,STATUS)
          IF (INULL .OR. (.NOT.ISEL)) GO TO 230
          YWORK  = XTEMP
          YWORK  = YWORK - YMEAN
          YSTD   = YSTD + YWORK*YWORK
          IF (XTEMP-CR) 190,160,160
C
C ... VARIATE VALUE ABOVE MIDDLE  INTERVAL
C
  160     DO 180 J = LDN,LUP
              IF (XTEMP.GE.CINT(J)) GO TO 170
              IFREQ(J) = IFREQ(J) + 1
              GO TO 230
  170         CONTINUE
  180     CONTINUE
          NUP    = NUP + 1
          GO TO 230
C
C ... VARIATE VALUE BELOW MIDDLE  INTERVAL
C
  190     IF (LK0) GO TO 220
          DO 210 J = 1,K0
              JJ     = NOB - J
              IF (XTEMP.LT.CINT(JJ)) GO TO 200
              JJ     = JJ + 1
              IFREQ(JJ) = IFREQ(JJ) + 1
              GO TO 230
  200         CONTINUE
  210     CONTINUE
  220     NDN    = NDN + 1
  230     CONTINUE
  240 CONTINUE
  250 IFREQ(1) = NDN
      IFREQ(K2) = NUP
      XMEAN  = YMEAN
      IF (NACT.EQ.1) THEN
         YSTD = 0
      ELSE
         YSTD   = DSQRT(YSTD/(NACT-1))
      ENDIF
      XSTD   = YSTD
      
C      CALL STTPUT('Usually routine ends here',ISTAT)
C      WRITE(6,*) NACT
C      WRITE(6,*) XMEAN, XSTD
     
      
      RETURN
  260 CR     = (XMAX+XMIN)*0.50
      CINT(1) = CR
  270 CONTINUE      
      DO 280 I = 1,N
          CALL TBERDR(TID,I,ICOL,XTEMP,INULL,STATUS)
          IF (XTEMP.LT.CR) NDN    = NDN + 1
  280 CONTINUE
      NUP    = N - NDN
      GO TO 250
  290 IERR   = 1 
      GO TO 310
  300 IERR   = 2 
  310 NACT   = 0
C      CALL STTPUT('NACT set to 0',ISTAT)      
C      CALL STTPUT('END OF THIS ROUTINE REACHED',ISTAT)
C      WRITE(6,*) NACT
C      WRITE(6,*) XMEAN, XSTD

      RETURN
      END

      SUBROUTINE TDDSTA
     .  (TID,ICOL,N,K2,ICLASS,CINT,IFREQ,NACT,XMEAN,XSTD)
C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
C.PURPOSE
C COMPUTE FREQUENCY DISTRIBUTION FOR VALUES IN THE ARRAY.
C OTHER STATISTICAL PARAMETERS RETURNED BY THE ROUTINE ARE
C MEAN, STD.DEV., MINIMUM AND MAXIMUM.
C       DOUBLE PRECISION VERSION
C
C
C------------------------------------------------------------------
      IMPLICIT NONE
      INTEGER TID        ! IN : TABLE IDENT
      INTEGER ICOL       ! IN : COLUMN INDEX
      INTEGER N          ! IN : number of input data
      INTEGER K2         ! IN : number of classes
      INTEGER ICLASS     !  IN :
      REAL CINT(K2)      !  IN :
      INTEGER IFREQ(K2)  ! IN  :
      INTEGER NACT       ! OUT : actual number of points
      REAL XMEAN         ! OUT : mean value
      REAL XSTD          ! OUT : standard deviation
C
      LOGICAL LK0
      INTEGER I,NUP,NDN,K0,STATUS
      INTEGER JJ,IERR,NOB,LDN
      INTEGER LUP,J
      LOGICAL ISEL, INULL
      REAL XMIN          ! Why real in double version?
      REAL XMAX          ! Why real in double version?
      REAL CR            ! Why real in double version?
      REAL ALPHA         ! Why real in double version?
      REAL STEP          ! Why real in double version?
      DOUBLE PRECISION YSTD
      DOUBLE PRECISION YMEAN
      DOUBLE PRECISION YWORK
      DOUBLE PRECISION XTEMP
C
      NACT   = 0
      YSTD   = 0.00D+0
      YMEAN  = 0.00D+0
      IF (K2-2) 290,10,10
   10 IF (N-1) 300,20,20
C
C ... ZERO  FREQUENCIES
C
   20 CONTINUE
      DO 30 I = 1,K2
          IFREQ(I) = 0
   30 CONTINUE
      NUP    = 0
      NDN    = 0
      K0     = K2 - 2
      LK0    = K0 .EQ. 0
C
C ... DETERMINE  MIN  AND  MAX, MEAN AND STD DEV.
C
      I      = 1
      CALL TBSGET(TID,I,ISEL,STATUS)
      CALL TBERDD(TID,I,ICOL,XTEMP,INULL,STATUS)
   40     CONTINUE
          IF (.NOT.INULL .AND. ISEL) GO TO 50
          I      = I + 1
          IF (I.GT.N) GOTO 65
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDD(TID,I,ICOL,XTEMP,INULL,STATUS)
          GO TO 40 
   50 CONTINUE
      XMIN   = XTEMP
      XMAX   = XTEMP
      DO 60 J = I,N
          CALL TBSGET(TID,J,ISEL,STATUS)
          CALL TBERDD(TID,J,ICOL,XTEMP,INULL,STATUS)
          IF (.NOT.INULL .AND. ISEL) THEN
              NACT   = NACT + 1
              YMEAN  = YMEAN + XTEMP
              IF (XTEMP.LT.XMIN) XMIN   = XTEMP
              IF (XTEMP.GT.XMAX) XMAX   = XTEMP
          END IF
   60 CONTINUE
   65 CONTINUE
      IF (NACT.EQ.0) THEN
          XMEAN  = 0.
          XSTD   = 0.
          RETURN
      ELSE
          YMEAN  = YMEAN/DBLE(NACT)
      END IF
      IF (ICLASS.EQ.1) GO TO 120
      IF (LK0) GO TO 260
C
C ... CALCULATE CLASS INTERVALS WITH EQUAL SPACING
C
      CR     = XMAX - XMIN
      ALPHA  = 0.001
      IF (CR.GT.0.0) GO TO 100
      IF (XMIN) 80,90,80
   80 CR     = XMIN*ALPHA
      GO TO 100
   90 CR     = ALPHA
  100 STEP   = CR* (1.00+ALPHA)/FLOAT(K0)
      CR     = XMIN - 0.50*ALPHA*CR
      CINT(1) = CR
      DO 110 I = 1,K0
          XTEMP  = CR
          CR     = XTEMP + STEP
          CINT(I+1) = CR
  110 CONTINUE
      GO TO 150
C
C ... CLASS INTERVALS SUPPLIED--CHECK IN ASCENDING ORDER
C
  120 CR     = CINT(1)
C      write(6,*) XMAX,XMIN
      IF (LK0) GO TO 270
      JJ     = K0 + 1
      DO 140 I = 2,JJ
          XTEMP  = CINT(I)
C 
C Typo origin of error? PN 12/98
C          IF (CR.LT.XTEMP) GO TO 130
C  LE. instead of LT.?
          IF (CR.LE.XTEMP) GO TO 130
          IERR   = 3
          GO TO 310

  130     CR     = XTEMP
  140 CONTINUE
C
C ... DETERMINE CLASS TO WHICH EACH CASE BELONGS
  150 NOB    = K2/2
      CR     = CINT(NOB)
      LDN    = NOB + 1
      LUP    = K2 - 1
      K0     = NOB - 1
      LK0    = K0 .EQ. 0
      DO 240 I = 1,N
          CALL TBSGET(TID,I,ISEL,STATUS)
          CALL TBERDD(TID,I,ICOL,XTEMP,INULL,STATUS)
          IF (INULL .OR. (.NOT.ISEL)) GOTO 230
          YWORK  = XTEMP
          YWORK  = YWORK - YMEAN
          YSTD   = YSTD + YWORK*YWORK
          IF (XTEMP-CR) 190,160,160
C
C ... VARIATE VALUE ABOVE MIDDLE  INTERVAL
C
  160     DO 180 J = LDN,LUP
              IF (XTEMP.GE.CINT(J)) GO TO 170
              IFREQ(J) = IFREQ(J) + 1
              GO TO 230
  170         CONTINUE
  180     CONTINUE
          NUP    = NUP + 1
          GO TO 230
C
C ... VARIATE VALUE BELOW MIDDLE  INTERVAL
C
  190     IF (LK0) GO TO 220
          DO 210 J = 1,K0
              JJ     = NOB - J
              IF (XTEMP.LT.CINT(JJ)) GO TO 200
              JJ     = JJ + 1
              IFREQ(JJ) = IFREQ(JJ) + 1
              GO TO 230

  200         CONTINUE
  210     CONTINUE
  220     NDN    = NDN + 1
  230     CONTINUE
  240 CONTINUE
  250 IFREQ(1) = NDN
      IFREQ(K2) = NUP
      XMEAN  = YMEAN
      IF (NACT.EQ.1) THEN
         YSTD = 0
      ELSE
         YSTD   = DSQRT(YSTD/(NACT-1))
      ENDIF
      XSTD   = YSTD
      RETURN
  260 CR     = (XMAX+XMIN)*0.50
      CINT(1) = CR
  270 CONTINUE
      DO 280 I = 1,N
          CALL TBERDD(TID,I,ICOL,XTEMP,INULL,STATUS)
          IF (XTEMP.LT.CR) NDN    = NDN + 1
  280 CONTINUE
      NUP    = N - NDN
      GO TO 250
  290 IERR   = 1
      GO TO 310
  300 IERR   = 2
  310 NACT   = 0
      RETURN
      END