C @(#)ustat.for	17.1.1.1 (ES0-DMD) 01/25/02 17:12:45
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 @(#)ustat.for	5.1 (ESO-IPG) 4/5/93 15:15:11
       SUBROUTINE AVEVAR(DATA,N,AVE,VAR)
       INTEGER J,N
       REAL DATA(N),AVE,VAR,S
       AVE = 0.
       VAR = 0.
       DO 10 J = 1,N
           AVE = AVE+DATA(J)
10       CONTINUE
       AVE=AVE/N
       DO 20 J = 1,N
           S = DATA(J) - AVE
           VAR =  VAR+S*S
20       CONTINUE
       VAR = VAR/(N-1)
       RETURN
       END



       SUBROUTINE CRANK(N,W,INDX,IRANK,S)
       INTEGER  N,INDX(N),J,JI,JT
       REAL IRANK(N),W(N),RANK,S,T
       S = 0
       J = 1
1       IF (J.LT.N) THEN
          IF(W(INDX(J)).NE.W(INDX(J+1))) THEN
              IRANK(INDX(J)) = J
              J = J+1
          ELSE
              DO 10 JT = J+1,N
                    IF(W(INDX(JT)).NE.W(INDX(J))) GOTO 20

10            CONTINUE
             JT = N+1
20             RANK = 0.5*(J+JT-1)
             DO 30 JI = J,JT-1
                 IRANK(INDX(JI)) = RANK
30             CONTINUE
             T = JT-J
             S = S + T**3-T
             J = JT
          ENDIF
       GOTO 1
       ENDIF
       IF(J.EQ.N) IRANK(INDX(N)) = N
       RETURN
       END




       SUBROUTINE INDEXX(N,ARRIN,INDX)
       INTEGER J,N,I
       REAL ARRIN(N),Q,INDXT
       INTEGER L,IR,INDX(N)
       DO J=1,N
          INDX(J)=J
       ENDDO
       L = N/2 + 1       
       IR = N
10      CONTINUE
       IF (L.GT.1)THEN
           L = L-1
           INDXT = INDX(L)
           Q = ARRIN(INDX(L))
       ELSE
           INDXT = INDX(IR)
           Q = ARRIN(INDX(IR))
           INDX(IR) = INDX(1)
           IR = IR-1
            IF (IR.EQ.1) THEN
              INDX(1) = INDXT
              RETURN
           ENDIF
       ENDIF       
       I = L
       J = L+L
20      IF (J.LE.IR) THEN
            IF(J.LT.IR)THEN
               IF(ARRIN(INDX(J)).LT.ARRIN(INDX(J+1)))J=J+1
            ENDIF
            IF (Q.LE.ARRIN(INDX(J)))THEN
                INDX(I) = INDX(J)
                I = J
                J = J+J
            ELSE
                J = IR+1
            ENDIF
       GOTO 20
       ENDIF
       INDX(I) = INDXT
       GOTO 10
       END




        SUBROUTINE SORT(LISTE,NBINT)
       INTEGER NBINT,I1,J,I2,JNEW
       REAL LISTE(NBINT)
       J=NBINT-1
2       CONTINUE
           JNEW=0
           DO 3 I1=1,J
              I2=I1+1
              IF(LISTE(I1).LE.LISTE(I2)) GOTO 3
               CALL SWAP(LISTE(I1),LISTE(I2))
               JNEW=I1-1
3           CONTINUE
            J=JNEW
       IF(J.GT.0) GOTO 2
       END

       SUBROUTINE SWAP(C1,C2)
       REAL*4 C1,C2,KAUX
       KAUX=C1
       C1=C2
       C2=KAUX
        END
       
       SUBROUTINE TMAP(TID,NROW,COL,X,NREAL)
       INTEGER TID,NROW,ISTAT,I,COL,NREAL
       LOGICAL NULL,ISEL
       REAL X(NROW)
       NREAL = 0
       DO I=1,NROW
          CALL TBSGET(TID,I,ISEL,ISTAT)
          IF (ISEL)  THEN 
            NREAL = NREAL+1
            CALL TBERDR(TID,I,COL,X(NREAL),NULL,ISTAT)
            IF (NULL) NREAL = NREAL-1
          ENDIF
       ENDDO
       RETURN
       END


       SUBROUTINE TMAP2(TID,NROW,COL1,COL2,X1,X2,NREAL)
       INTEGER TID,NROW,COL1,COL2,ISTAT,I,NREAL
       LOGICAL NULL1,NULL2,ISEL
       REAL*4 X1(NROW),X2(NROW)
       NREAL = 0 
       DO I=1,NROW
          CALL TBSGET(TID,I,ISEL,ISTAT)
          IF (ISEL)  THEN 
            NREAL = NREAL+1
            CALL TBERDR(TID,I,COL1,X1(NREAL),NULL1,ISTAT)
            CALL TBERDR(TID,I,COL2,X2(NREAL),NULL2,ISTAT)
            IF (NULL1 .OR. NULL2) NREAL = NREAL-1 
          ENDIF 
       ENDDO
       RETURN
       END