diff --git a/src/basics/Vector3.cpp b/src/basics/Vector3.cpp
index bcd81d7..92db091 100644
--- a/src/basics/Vector3.cpp
+++ b/src/basics/Vector3.cpp
@@ -2,6 +2,7 @@
 
 #include <cmath>
 #include "PackStream.h"
+#include "RandomGenerator.h"
 
 const Vector3 paysages::basics::VECTOR_ZERO(0.0, 0.0, 0.0);
 const Vector3 paysages::basics::VECTOR_DOWN(0.0, -1.0, 0.0);
@@ -81,3 +82,15 @@ Vector3 Vector3::midPointTo(const Vector3 &other) const
 {
     return Vector3((other.x + x) * 0.5, (other.y + y) * 0.5, (other.z + z) * 0.5);
 }
+
+Vector3 Vector3::randomInSphere(double radius, bool only_surface)
+{
+    // TODO More uniform spatial repartition
+    //  The current randomization clusters result near the center and at the poles
+    VectorSpherical vec = {
+        only_surface ? radius : RandomGenerator::random() * radius,
+        (RandomGenerator::random() - 0.5) * M_PI,
+        RandomGenerator::random() * M_2PI
+    };
+    return Vector3(vec);
+}
diff --git a/src/basics/Vector3.h b/src/basics/Vector3.h
index 2c8a40c..c380241 100644
--- a/src/basics/Vector3.h
+++ b/src/basics/Vector3.h
@@ -68,6 +68,13 @@ public:
 
     VectorSpherical toSpherical() const;
 
+    /**
+     * Produce a random vector in a sphere domain.
+     *
+     * If *only_surface* is true, produce a vector with *radius* as length.
+     */
+    static Vector3 randomInSphere(double radius=1.0, bool only_surface=false);
+
 public:
     // TODO Make private
     double x;
diff --git a/src/basics/basics_global.h b/src/basics/basics_global.h
index 591f672..3711655 100644
--- a/src/basics/basics_global.h
+++ b/src/basics/basics_global.h
@@ -38,5 +38,8 @@ using namespace paysages::basics;
 #ifndef M_PI_4
 #define M_PI_4 0.78539816339744830962
 #endif
+#ifndef M_2PI
+#define M_2PI 6.28318530717958647692
+#endif
 
 #endif // BASICS_GLOBAL_H
diff --git a/src/tests/Vector3_Test.cpp b/src/tests/Vector3_Test.cpp
index a03d157..c6be931 100644
--- a/src/tests/Vector3_Test.cpp
+++ b/src/tests/Vector3_Test.cpp
@@ -8,3 +8,17 @@ TEST(Vector3, midPointTo)
     Vector3 vm = v1.midPointTo(v2);
     EXPECT_VECTOR3_COORDS(vm, 2.5, 2.25, 3.5);
 }
+
+TEST(Vector3, randomInSphere)
+{
+    Vector3 v;
+
+    v = Vector3::randomInSphere();
+    EXPECT_LE(v.getNorm(), 1.0);
+
+    v = Vector3::randomInSphere(0.01);
+    EXPECT_LE(v.getNorm(), 0.01);
+
+    v = Vector3::randomInSphere(0.5, true);
+    EXPECT_DOUBLE_EQ(v.getNorm(), 0.5);
+}