diff --git a/doc/tutorial.qbk b/doc/tutorial.qbk
index d14c340..84a0e4f 100644
--- a/doc/tutorial.qbk
+++ b/doc/tutorial.qbk
@@ -26,3 +26,10 @@
 [password]
 
 [endsect]
+
+[section Generating quasi-random line-sphere intersections]
+
+[import ../example/intersections.cpp]
+[intersections]
+
+[endsect]
diff --git a/example/Jamfile.v2 b/example/Jamfile.v2
index a57fe3f..c598e78 100644
--- a/example/Jamfile.v2
+++ b/example/Jamfile.v2
@@ -10,3 +10,4 @@
 run die.cpp ;
 run weighted_die.cpp ;
 run password.cpp /boost//random ;
+run intersections.cpp /boost//random ;
diff --git a/example/intersections.cpp b/example/intersections.cpp
new file mode 100644
index 0000000..f4c1435
--- /dev/null
+++ b/example/intersections.cpp
@@ -0,0 +1,74 @@
+// intersections.cpp
+//
+// Copyright (c) 2018
+// Justinas V. Daugmaudis
+//
+// 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)
+
+//[intersections
+/*`
+    For the source of this example see
+    [@boost://libs/random/example/intersections.cpp intersections.cpp].
+
+    This example demonstrates generating quasi-randomly distributed chord
+    entry and exit points on an '''<code>S<superscript>2</superscript></code>''' sphere.
+
+    First we include the headers we need for __niederreiter_base2
+    and __uniform_01 distribution.
+ */
+
+#include <boost/random/niederreiter_base2.hpp>
+#include <boost/random/uniform_01.hpp>
+
+#include <boost/math/constants/constants.hpp>
+
+#include <boost/tuple/tuple.hpp>
+
+/*`
+  We use 4-dimensional __niederreiter_base2 as a source of randomness.
+ */
+boost::random::niederreiter_base2 gen(4);
+
+
+int main()
+{
+  typedef boost::tuple<double, double, double> point_t;
+
+  const std::size_t n_points = 100; // we will generate 100 points
+
+  std::vector<point_t> points;
+  points.reserve(n_points);
+
+  boost::random::uniform_01<double> dist;
+
+  for (std::size_t i = 0; i != n_points; ++i)
+  {
+    /*`
+      Using formula from J. Rovira et al., "Point sampling with uniformly distributed lines", 2005
+      to compute uniformly distributed chord entry and exit points on the surface of a sphere.
+    */
+    double cos_theta = 1 - 2 * dist(gen);
+    double sin_theta = std::sqrt(1 - cos_theta * cos_theta);
+    double phi = boost::math::constants::two_pi<double>() * dist(gen);
+    double sin_phi = std::sin(phi), cos_phi = std::cos(phi);
+
+    point_t point_on_sphere(sin_theta*sin_phi, cos_theta, sin_theta*cos_phi);
+
+    /*`
+      Here we assume that our sphere is a unit sphere at (0,0,0). If your sphere was
+      different then now would be the time to scale and translate the `point_on_sphere`.
+    */
+
+    points.push_back(point_on_sphere);
+  }
+
+  /*`
+    Vector `points` now holds generated 3D points on a unit sphere.
+  */
+
+  return 0;
+}
+
+//]