function value = tfn ( x, fx )

%% TFN calculates the T-function of Owen.
%
%  Modified:
%
%    20 January 2008
%
%  Author:
%
%    JC Young, Christoph Minder
%    MATLAB version by John Burkardt
%
%  Reference:
%
%    MA Porter, DJ Winstanley,
%    Remark AS R30:
%    A Remark on Algorithm AS76:
%    An Integral Useful in Calculating Noncentral T and Bivariate
%    Normal Probabilities,
%    Applied Statistics,
%    Volume 28, Number 1, 1979, page 113.
%
%    JC Young, Christoph Minder,
%    Algorithm AS 76:
%    An Algorithm Useful in Calculating Non-Central T and
%    Bivariate Normal Distributions,
%    Applied Statistics,
%    Volume 23, Number 3, 1974, pages 455-457.
%
%  Parameters:
%
%    Input, real X, FX, the parameters of the function.
%
%    Output, real VALUE, the value of the T-function.
%
  ng = 5;

  r = [ ...
    0.1477621, ...
    0.1346334, ...
    0.1095432, ...
    0.0747257, ...
    0.0333357 ];
  tp = 0.159155;
  tv1 = 1.0E-35;
  tv2 = 15.0;
  tv3 = 15.0;
  tv4 = 1.0E-05;
  u = [ ...
    0.0744372, ...
    0.2166977, ...
    0.3397048, ...
    0.4325317, ...
    0.4869533 ];
%
%  Test for X near zero.
%
  if ( abs ( x ) < tv1 )
    value = tp * atan ( fx );
    return
  end
%
%  Test for large values of abs(X).
%
  if ( tv2 < abs ( x ) )
    value = 0.0;
    return
  end
%
%  Test for FX near zero.
%
  if ( abs ( fx ) < tv1 )
    value = 0.0;
    return
  end
%
%  Test whether abs ( FX ) is so large that it must be truncated.
%
  xs = - 0.5 * x * x;
  x2 = fx;
  fxs = fx * fx;
%
%  Computation of truncation point by Newton iteration.
%
  if ( tv3 <= log ( 1.0 + fxs ) - xs * fxs )

    x1 = 0.5 * fx;
    fxs = 0.25 * fxs;

    while ( 1 )

      rt = fxs + 1.0;

      x2 = x1 + ( xs * fxs + tv3 - log ( rt ) ) ...
      / ( 2.0 * x1 * ( 1.0 / rt - xs ) );

      fxs = x2 * x2;

      if ( abs ( x2 - x1 ) < tv4 )
        break
      end

      x1 = x2;

    end

  end
%
%  Gaussian quadrature.
%
  rt = 0.0;

  for i = 1 : ng

    r1 = 1.0 + fxs * ( 0.5 + u(i) )^2;
    r2 = 1.0 + fxs * ( 0.5 - u(i) )^2;

    rt = rt + r(i) * ( exp ( xs * r1 ) / r1 + exp ( xs * r2 ) / r2 );

  end

  value = rt * x2 * tp;