* +-======-+ 
*  Copyright (c) 2003-2007 United States Government as represented by 
*  the Admistrator of the National Aeronautics and Space Administration.  
*  All Rights Reserved.
*  
*  THIS OPEN  SOURCE  AGREEMENT  ("AGREEMENT") DEFINES  THE  RIGHTS  OF USE,
*  REPRODUCTION,  DISTRIBUTION,  MODIFICATION AND REDISTRIBUTION OF CERTAIN 
*  COMPUTER SOFTWARE ORIGINALLY RELEASED BY THE UNITED STATES GOVERNMENT AS 
*  REPRESENTED BY THE GOVERNMENT AGENCY LISTED BELOW ("GOVERNMENT AGENCY").  
*  THE UNITED STATES GOVERNMENT, AS REPRESENTED BY GOVERNMENT AGENCY, IS AN 
*  INTENDED  THIRD-PARTY  BENEFICIARY  OF  ALL  SUBSEQUENT DISTRIBUTIONS OR 
*  REDISTRIBUTIONS  OF THE  SUBJECT  SOFTWARE.  ANYONE WHO USES, REPRODUCES, 
*  DISTRIBUTES, MODIFIES  OR REDISTRIBUTES THE SUBJECT SOFTWARE, AS DEFINED 
*  HEREIN, OR ANY PART THEREOF,  IS,  BY THAT ACTION, ACCEPTING IN FULL THE 
*  RESPONSIBILITIES AND OBLIGATIONS CONTAINED IN THIS AGREEMENT.
*  
*  Government Agency: National Aeronautics and Space Administration
*  Government Agency Original Software Designation: GSC-15354-1
*  Government Agency Original Software Title:  GEOS-5 GCM Modeling Software
*  User Registration Requested.  Please Visit http://opensource.gsfc.nasa.gov
*  Government Agency Point of Contact for Original Software:  
*  			Dale Hithon, SRA Assistant, (301) 286-2691
*  
* +-======-+ 
***************************************************************************************
*	$Id: studt.gs,v 1.1 2010/08/17 18:24:18 ltakacs Exp $
*	Copyright (C) 2004 Bin Guan.
*	Distributed under GNU/GPL.
***************************************************************************************
function studt(arg)
*
* Calculate p for given dof and t in a two-tailed Student t-test. 
*
dof=subwrd(arg,1)
t=subwrd(arg,2)
pname=subwrd(arg,3)
if(t='')
usage()
return
endif
if(pname='')
pname='studtout'
endif

pi=3.141593
pid2=pi/2.0

t=math_abs(t)
w=t/math_sqrt(dof)
th=math_atan(w)

if(dof=1)
p=1-th/pid2
pname'='p
return
endif

sth=math_sin(th)
cth=math_cos(th)
if(math_fmod(dof,2)=1)
p=1-(th+sth*cth*statcom(cth*cth,2,dof-3,-1))/pid2
else
p=1-sth*statcom(cth*cth,1,dof-3,-1)
endif
pname'='p
return
***************************************************************************************
*
* Calculate some statistic. 
*
function statcom(q,i,j,b)
zz=1
z=zz

k=i
while(k<=j)
zz=zz*q*k/(k-b)
z=z+zz
k=k+2
endwhile

return z
***************************************************************************************
function usage()
*
* Print usage information.
*
say '  Two-tailed Student t-test.'
say ''
say '  Usage: studt <dof> <t> [<p>]'
say '         <dof>: degree of freedom.'
say '         <t>: t-statistic.'
say '         <p>: probability. Defaults to studtout.'
say ''
say '  Copyright (C) 2004 Bin Guan.'
say '  Distributed under GNU/GPL.'
return