C @(#)tstest.for	17.1.1.1 (ESO-DMD) 01/25/02 17:12:48
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 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 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  5 JUL 1989
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: M.PERON
C
C.IDENTIFICATION
C
C  Program TSTEST
C
C.KEYWORDS
C
C means,Student
C
C.PURPOSE
C
C   Student-T Test for different means 
C
C.REFERENCES  Numerical Recipes p465
C 
C  010621		last modif
C 
C-----------------------------------------------------------
C
       PROGRAM TSTEST
C 
       REAL*4 T,PROB
C 
       INTEGER*8  ADDRIN(2)
C 
       INTEGER MADRID
       INTEGER PARVAL,NPAR,ISTAT,TID,NCOL,NROW,NACOL,NAROW,COL1,COL2
       INTEGER NBYTES,NSC,KUN,NREAL
C 
       CHARACTER*80 INTAB,OUTPUT1,OUTPUT2,OUT,OUTPUT3,OUTPUT0
       CHARACTER*17 COLIND,COLDEP
C 
       INCLUDE 'MID_INCLUDE:TABLES.INC'
       COMMON/VMR/MADRID(1)
C 
       INCLUDE 'MID_INCLUDE:TABLED.INC'
       DATA PARVAL/8/
       DATA OUTPUT0
     1   /'Null Hypothesis:The two samples have consistent means'/
       DATA OUTPUT1
     1    /'Student-t Statistic '/
       DATA OUTPUT2
     2    /'---------------------------------------------------------'/
       DATA OUTPUT3
     3  /' Probability of exceeding this value under Null Hypothesis'/
C 
C 
       CALL STSPRO('TSTEST')
       CALL TDPGET(PARVAL,NPAR,ISTAT)
       IF  (ISTAT.NE.0) GOTO 80
       INTAB = TPARBF(1)
       COLIND = TPARBF(2)
       COLDEP = TPARBF(3)
       CALL TBTOPN(INTAB,F_U_MODE,TID,ISTAT)       
       IF(ISTAT.NE.0) GOTO 80
       CALL TBIGET(TID,NCOL,NROW,NSC,NACOL,NAROW,ISTAT)
        CALL TBCSER(TID,COLIND,COL1,ISTAT)
         IF (COL1.LT.0) THEN
         CALL STTPUT('Column not found...',ISTAT)
          GOTO 80
        ENDIF
        CALL TBCSER(TID,COLDEP,COL2,ISTAT)
        IF (COL2.LT.0) THEN
          CALL STTPUT('Column not found...',ISTAT)
          GOTO 80
        ENDIF
       NBYTES = 4*NROW
       CALL TDMGET(NBYTES,ADDRIN(1),ISTAT)
       CALL TDMGET(NBYTES,ADDRIN(2),ISTAT)
       CALL TMAP2(TID,NROW,COL1,COL2,MADRID(ADDRIN(1)),
     1                MADRID(ADDRIN(2)),NREAL)
       CALL TTEST(MADRID(ADDRIN(1)),MADRID(ADDRIN(2)),NREAL,T,PROB)
       WRITE(OUT,100)T
       CALL STTPUT(OUTPUT0,ISTAT)
       CALL STTPUT(OUTPUT2,ISTAT)
       CALL STTPUT(OUTPUT1,ISTAT)
       CALL STTPUT(OUT,ISTAT)
       CALL STTPUT(OUTPUT2,ISTAT)
       CALL STTPUT(OUTPUT3,ISTAT)
       WRITE(OUT,100)PROB
       CALL STTPUT(OUT,ISTAT)
       CALL STTPUT(OUTPUT2,ISTAT)
       CALL STKWRR('OUTPUTR',T,1,1,KUN,ISTAT)
       CALL STKWRR('OUTPUTR',PROB,2,1,KUN,ISTAT)
80        CONTINUE
       CALL STSEPI
C 
100     FORMAT(1X,G15.6,15X,G15.6)
       END       



       SUBROUTINE TTEST(DATA1,DATA2,NROW,T,PROB)
       INTEGER NROW
       REAL*4 DATA1(NROW),DATA2(NROW),T,PROB
       REAL*4 AVE1,AVE2,VAR,VAR1,VAR2,DF,BETAI
       CALL AVEVAR(DATA1,NROW,AVE1,VAR1)
       CALL AVEVAR(DATA2,NROW,AVE2,VAR2)
       DF = 2*NROW-2
       VAR = ((NROW-1)*VAR1+(NROW-1)*VAR2)/DF
       T = (AVE1-AVE2)/SQRT(VAR*(1./NROW+1./NROW))
       PROB = BETAI(0.5*DF,0.5,DF/(DF+T**2))
       RETURN
       END