diff --git a/doc/sf_and_dist/owens_t.qbk b/doc/sf_and_dist/owens_t.qbk
new file mode 100644
index 000000000..06a73d015
--- /dev/null
+++ b/doc/sf_and_dist/owens_t.qbk
@@ -0,0 +1,112 @@
+[section:owens_t Owen's T function]
+
+[h4 Synopsis]
+
+``
+#include <boost/math/special_functions/owens_t.hpp>
+``
+
+  namespace boost{ namespace math{
+  
+  template <class T>
+  ``__sf_result`` owens_t(T h, T a);
+  
+  template <class T, class ``__Policy``>
+  ``__sf_result`` owens_t(T h, T a, const ``__Policy``&);
+  
+  }} // namespaces
+  
+[h4 Description]
+
+Returns the
+[@http://en.wikipedia.org/wiki/Owen%27s_T_function Owens_t function]
+of ['h] and ['a].
+
+[optional_policy]
+
+[sixemspace][sixemspace][equation owens_t]
+
+[$plot_owens_t.png]
+
+The function `owens_t(h, a)` gives the probability
+of the event ['(X > h and 0 < Y < a * X)],
+where ['X] and ['Y] are independent standard normal random variables.
+
+For h and a > 0, T(h,a),
+gives the volume of an uncorrelated bivariate normal distribution
+with zero means and unit variances over the area between
+['y = ax] and ['y = 0] and to the right of ['x = h].
+
+That is the area shaded in the figure below (Owens 1956).
+
+[graph owens_integration_area]
+
+and is also illustrated by a 3D plot.
+
+[$plot_owens_3d_xyp.png]
+
+This function is used in the computation of the __skew_normal_distrib.
+It is also used in the computation of bivariate and
+multivariate normal distribution probabilities.
+The return type of this function is computed using the __arg_pomotion_rules:
+the result is of type `double` when T is an integer type, and type T otherwise.
+
+[h4 Accuracy]
+
+The Patefield-Tandy algorithm provides six methods of evalualution;
+the best method is selected according to the values of ['a] and ['h].
+See the original paper and the source in
+[@../../../../../boost/math/special_functions/owens_t.hpp owens_t.hpp] for details.
+
+Over the built-in types and range tested,
+errors are less than 10 * std::numeric_limits<RealType>::epsilon().
+
+There is no fully generic version of this function:  all the implementations
+are tuned to specific accuracy levels, the most precise of which delivers
+about 30-digits of precision.
+
+For arbitrary precision calculation we just provide a
+T1 method as the default case and warn that only absolute errors can be guaranteed.
+We could maybe include T1b as well for h < 1,
+but the real problems come when h gets large.
+
+[h4 Testing]
+
+Test data was generated by T1 and T4.
+
+* T1 was rejected if the result was too small compared to `atan(a)` (ie cancellation),
+* T4 was rejected if there was no convergence,
+* Both were rejected if they didn't agree.
+
+Over the built-in types and range tested,
+errors are less than 10 std::numeric_limits<RealType>::epsilon().
+
+However, that there was a whole domain (large ['h], small ['a])
+where it was not possible to generate any reliable test values
+(all the methods got rejected for one reason or another).
+
+There are also two sets of sanity tests: spot values are computed using __Mathematica and ???? TODO.
+
+[h4 Implementation]
+
+The function was proposed and evaluated by
+D. B. Owen, Tables for computing bivariate normal probabilities, 
+Ann. Math. Statist., 27, 1075-1090 (1956).
+
+The algorithms of Patefield, M. and Tandy, D.
+"Fast and accurate Calculation of Owen's T-Function", Journal of Statistical Software, 5 (5), 1 - 25 (2000)
+are adapted for C++ with arbitrary RealType.
+
+However, at present, some coefficients are limited to long double precision of 20 decimal digits,
+so this places a constraint on the accuracy of results,
+even if RealType provides high precision, for example 50 decimal digits.   
+
+[endsect] [/section:owens_t The owens_t Function]
+
+[/ 
+  Copyright 2012 Bejamin Sobotta, John Maddock and Paul A. Bristow.
+  Distributed under the Boost Software License, Version 1.0.
+  (See accompanying file LICENSE_1_0.txt or copy at
+  http://www.boost.org/LICENSE_1_0.txt).
+]
+