function theResult = inv_student_t(p, dof, tol)

% inv_student_t -- Inverse Student_t probability.
%  inv_student_t(p, dof, tol) returns the Student-t value
%   whose integral [-t, t] is probability p, for dof degrees
%   of freedom, computed to tolerance tol (default = 1e-8).
 
% Copyright (C) 1999 Dr. Charles R. Denham, ZYDECO.
%  All Rights Reserved.
%   Disclosure without explicit written consent from the
%    copyright owner does not constitute publication.
 
% Version of 25-Jan-1999 14:06:53.

MAXIMUM_ITERATIONS = 100;
SUB_SAMPLES = 100;
DEFAULT_TOLERANCE = 1.0e-8;

if nargin < 2, help(mfilename), return, end
if nargin < 3, tol = DEFAULT_TOLERANCE; end

if ischar(p), p = eval(p); end
if ischar(dof), dof = eval(dof); end
if ischar(tol), tol = eval(tol); end

if length(p) > 1 | length(dof) > 1
	if length(p) > 1
		result = zeros(size(p));
		if length(dof) == 1, dof = dof .* ones(size(p)); end
	elseif length(dof) > 1
		result = zeros(size(dof));
		if length(p) == 1, p = p .* ones(size(dof)); end
	end
	for i = 1:prod(size(p))
		result(i) = inv_student_t(p(i), dof(i), tol);
	end
	if nargout > 0
		theResult = result;
	else
		disp(result)
	end
	return
end

a = 0;
b = dof;

w = warning;
warning off;

while student_t(b, dof) < p
	a = b;
	b = 2 * b;
end

if a == 0, a = eps; end

result = [];

for i = 1:MAXIMUM_ITERATIONS
	delta = (b - a) ./ SUB_SAMPLES;
	a = a:delta:b;
	if a(1) == 0, a(1) = eps; end
	a(end) = b;
	c = student_t(a, dof);
	k = 1:length(c)-1;
	f = find(c(k) <= p & c(k+1) >= p);
	f = f(1);
	b = a(f+1);
	a = a(f);
	c = (c(f) + c(f+1))/2;
	err = c - p;
	if abs(err) < tol
		result = (a+b)/2;
		break;
	end
end

warning(w)

if nargout > 0
	theResult = result;
else
	disp(result)
end
