diff --git a/CMakeLists.txt b/CMakeLists.txt
index 5258282e4..eb425175d 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -57,8 +57,21 @@ if (SWIG_FOUND)
   else()
     message (STATUS "Searching for Ruby - found.")
   endif()
+
+  ########################################
+  # Include python
+  find_package(PythonLibs QUIET)
+  if (NOT PYTHONLIBS_FOUND)
+    message (STATUS "Searching for Python - not found.")
+  else()
+    message (STATUS "Searching for Python - found.")
+  endif()
 endif()
 
+# Location of "fake install folder" used in tests
+# Defined here at root scope so it is available for tests in src and test folders
+set(FAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/fake/install")
+
 #============================================================================
 # Configure the build
 #============================================================================
diff --git a/examples/angle_example.py b/examples/angle_example.py
new file mode 100644
index 000000000..bbe0fe0a1
--- /dev/null
+++ b/examples/angle_example.py
@@ -0,0 +1,43 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# This example will only work if the Python interface library was compiled and
+# installed.
+#
+# Modify the PYTHONPATH environment variable to include the ignition math
+# library install path. For example, if you install to /user:
+#
+# $ export PYTHONPATH=/usr/lib/python:$PYTHONPATH
+#
+
+import ignition.math
+
+print("PI in degrees = {}\n".format(ignition.math.Angle.Pi.Degree()))
+
+a1 = ignition.math.Angle(1.5707)
+a2 = ignition.math.Angle(0.7854)
+print("a1 = {} radians, {} degrees\n".format(a1.Radian(), a1.Degree()))
+print("a2 = {} radians, {} degrees\n".format(a2.Radian(), a2.Degree()))
+print("a1 * a2 = {} radians, {} degrees\n".format((a1 * a2).Radian(),
+      (a1 * a2).Degree()))
+print("a1 + a2 = {} radians, {} degrees\n".format((a1 + a2).Radian(),
+      (a1 + a2).Degree()))
+print("a1 - a2 = {} radians, {} degrees\n".format((a1 - a2).Radian(),
+      (a1 - a2).Degree()))
+
+a3 = ignition.math.Angle(15.707)
+print("a3 = {} radians, {} degrees\n".format(a3.Radian(), a3.Degree()))
+a3.Normalize()
+print("a3.Normalize = {} radians, {} degrees\n".format(a3.Radian(),
+                                                       a3.Degree()))
diff --git a/examples/vector2_example.py b/examples/vector2_example.py
new file mode 100644
index 000000000..ecdfc3120
--- /dev/null
+++ b/examples/vector2_example.py
@@ -0,0 +1,39 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# This example will only work if the Python interface library was compiled and
+# installed.
+#
+# Modify the PYTHONPATH environment variable to include the ignition math
+# library install path. For example, if you install to /user:
+#
+# $ export PYTHONPATH=/usr/lib/python:$PYTHONPATH
+#
+import ignition.math
+
+va = ignition.math.Vector2d(1, 2)
+vb = ignition.math.Vector2d(3, 4)
+vc = ignition.math.Vector2d(vb)
+
+print("va = {} {}\n".format(va.X(), va.Y()))
+print("vb = {} {}\n".format(vb.X(), vb.Y()))
+print("vc = {} {}\n".format(vc.X(), vc.Y()))
+
+vb += va
+print("vb += va: {} {}\n".format(vb.X(), vb.Y()))
+
+vb.Normalize()
+print("vb.Normalize = {} {}\n".format(vb.X(), vb.Y()))
+
+print("vb.Distance(va) = {}\n".format(vb.Distance(va)))
diff --git a/examples/vector3_example.py b/examples/vector3_example.py
new file mode 100644
index 000000000..5e0b22f9b
--- /dev/null
+++ b/examples/vector3_example.py
@@ -0,0 +1,34 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+# This example will only work if the Python interface library was compiled and
+# installed.
+#
+# Modify the PYTHONPATH environment variable to include the ignition math
+# library install path. For example, if you install to /user:
+#
+# $ export PYTHONPATH=/usr/lib/python:$PYTHONPATH
+#
+import ignition.math
+
+v1 = ignition.math.Vector3d(0, 0, 3)
+print("v =: {} {} {}\n".format(v1.X(), v1.Y(), v1.Z()))
+
+v2 = ignition.math.Vector3d(4, 0, 0)
+print("v2 = {} {} {}\n".format(v2.X(), v2.Y(), v2.Z()))
+
+v3 = v1 + v2
+print("v1 + v2 = {} {} {}\n".format(v3.X(), v3.Y(), v3.Z()))
+
+print("v1.Distance(v2) = {}\n".format(v1.Distance(v2)))
diff --git a/src/Angle_TEST.py b/src/Angle_TEST.py
new file mode 100644
index 000000000..00d1c7323
--- /dev/null
+++ b/src/Angle_TEST.py
@@ -0,0 +1,90 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#       http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+import math
+from ignition.math import Angle
+
+
+class TestAngle(unittest.TestCase):
+
+    def test_angle(self):
+
+        angle1 = Angle()
+        self.assertEqual(0.0, angle1.Radian())
+
+        angle1.SetDegree(90.0)
+        self.assertTrue(angle1 == Angle.HalfPi)
+
+        angle1.SetDegree(180.0)
+        self.assertTrue(angle1 == Angle.Pi)
+        self.assertFalse(angle1 == Angle.Pi + Angle(0.1))
+        self.assertTrue(angle1 == Angle.Pi + Angle(0.0001))
+        self.assertTrue(angle1 == Angle.Pi - Angle(0.0001))
+        self.assertTrue(Angle(0) == Angle(0))
+        self.assertTrue(Angle(0) == Angle(0.001))
+
+        angle1 = Angle(0.1) - Angle(0.3)
+        self.assertAlmostEqual(angle1.Radian(), -0.2)
+
+        angle = Angle(0.5)
+        self.assertEqual(0.5, angle.Radian())
+
+        angle.SetRadian(math.pi/2)
+        self.assertAlmostEqual(math.degrees(math.pi/2), angle.Degree())
+
+        angle.SetRadian(math.pi)
+        self.assertAlmostEqual(math.degrees(math.pi), angle.Degree())
+
+    def test_normalized_angles(self):
+
+        angle = Angle(Angle.Pi)
+        normalized = angle.Normalized()
+
+        angle.Normalized()
+        self.assertEqual(math.degrees(math.pi), angle.Degree())
+        self.assertEqual(normalized, angle)
+
+    def test_angle_operations(self):
+
+        angle = Angle(0.1) + Angle(0.2)
+        self.assertAlmostEqual(0.3, angle.Radian())
+
+        angle = Angle(0.1) * Angle(0.2)
+        self.assertAlmostEqual(0.02, angle.Radian())
+
+        angle = Angle(0.1) / Angle(0.2)
+        self.assertAlmostEqual(0.5, angle.Radian())
+
+        angle -= Angle(0.1)
+        self.assertAlmostEqual(0.4, angle.Radian())
+
+        angle += Angle(0.2)
+        self.assertAlmostEqual(0.6, angle.Radian())
+
+        angle *= Angle(0.5)
+        self.assertAlmostEqual(0.3, angle.Radian())
+
+        angle /= Angle(0.1)
+        self.assertAlmostEqual(3.0, angle.Radian())
+        self.assertTrue(angle == Angle(3))
+        self.assertTrue(angle != Angle(2))
+        self.assertTrue(angle < Angle(4))
+        self.assertTrue(angle > Angle(2))
+        self.assertTrue(angle >= Angle(3))
+        self.assertTrue(angle <= Angle(3))
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index 192fbddc7..dcc5bf3c1 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -26,9 +26,9 @@ if (SWIG_FOUND)
   set(CMAKE_SWIG_FLAGS "")
 
   include_directories(${PROJECT_SOURCE_DIR}/include)
+  include_directories(${PYTHON_INCLUDE_PATH})
 
   set(swig_files
-    ruby
     Angle
     GaussMarkovProcess
     Rand
@@ -50,27 +50,32 @@ if (RUBY_FOUND)
     # Generate the list if .i files
     list(APPEND swig_i_files ${swig_file}.i)
   endforeach()
+  list(APPEND ruby_tests ruby_TEST)
 
   # Turn on c++
-  set_source_files_properties(${swig_i_files} PROPERTIES CPLUSPLUS ON)
+  set_source_files_properties(${swig_i_files} ruby/ruby.i PROPERTIES CPLUSPLUS ON)
 
   # Create the ruby library
 
+  set(CMAKE_SWIG_OUTDIR "${CMAKE_BINARY_DIR}/lib/ruby")
   if(CMAKE_VERSION VERSION_GREATER 3.8.0)
-    SWIG_ADD_LIBRARY(math LANGUAGE ruby SOURCES ruby.i ${swig_i_files} ${sources})
+    SWIG_ADD_LIBRARY(math LANGUAGE ruby SOURCES ruby/ruby.i ${swig_i_files})
   else()
-    SWIG_ADD_MODULE(math ruby ruby.i ${swig_i_files} ${sources})
+    SWIG_ADD_MODULE(math ruby ruby/ruby.i ${swig_i_files})
   endif()
 
   # Suppress warnings on SWIG-generated files
-  target_compile_options(math PRIVATE 
+  target_compile_options(math PRIVATE
     $<$<CXX_COMPILER_ID:GNU>:-Wno-pedantic -Wno-shadow -Wno-maybe-uninitialized -Wno-unused-parameter>
     $<$<CXX_COMPILER_ID:Clang>:-Wno-shadow -Wno-maybe-uninitialized -Wno-unused-parameter>
     $<$<CXX_COMPILER_ID:AppleClang>:-Wno-shadow -Wno-maybe-uninitialized -Wno-unused-parameter>
   )
   target_include_directories(math SYSTEM PUBLIC ${RUBY_INCLUDE_DIRS})
 
-  SWIG_LINK_LIBRARIES(math ${RUBY_LIBRARY})
+  SWIG_LINK_LIBRARIES(math
+    ${RUBY_LIBRARY}
+    ignition-math${PROJECT_VERSION_MAJOR}
+  )
   target_compile_features(math PUBLIC ${IGN_CXX_${c++standard}_FEATURES})
   install(TARGETS math DESTINATION ${IGN_LIB_INSTALL_DIR}/ruby/ignition)
 
@@ -81,3 +86,73 @@ if (RUBY_FOUND)
    	  --gtest_output=xml:${CMAKE_BINARY_DIR}/test_results/${test}rb.xml)
   endforeach()
 endif()
+
+#################################
+# Create and install Python interfaces
+# Example usage
+# $ export PYTHONPATH=/ws/install/lib/python/:$PYTHONPATH
+if (PYTHONLIBS_FOUND)
+  set_source_files_properties(python/python.i PROPERTIES CPLUSPLUS ON)
+  set_source_files_properties(python/python.i PROPERTIES SWIG_FLAGS "-includeall")
+  set_source_files_properties(python/python.i PROPERTIES SWIG_MODULE_NAME "math")
+  set(SWIG_PY_LIB pymath)
+  set(SWIG_PY_LIB_OUTPUT math)
+
+  set(CMAKE_SWIG_OUTDIR "${CMAKE_BINARY_DIR}/lib/python")
+  if(CMAKE_VERSION VERSION_GREATER 3.8.0)
+    SWIG_ADD_LIBRARY(${SWIG_PY_LIB} LANGUAGE python SOURCES python/python.i)
+  else()
+    SWIG_ADD_MODULE(${SWIG_PY_LIB} python python/python.i)
+  endif()
+
+  SWIG_LINK_LIBRARIES(${SWIG_PY_LIB}
+    ${PYTHON_LIBRARIES}
+    ignition-math${PROJECT_VERSION_MAJOR}
+  )
+
+  # When using SWIG 3 to build a python interface there is an extra underscore between the name of
+  # the library given to swig_add_library macro and the generated .so
+  if(NOT ${SWIG_VERSION} VERSION_GREATER 4.0.0)
+    set(SWIG_PY_LIB "_${SWIG_PY_LIB}")
+    set(SWIG_PY_LIB_OUTPUT "_${SWIG_PY_LIB_OUTPUT}")
+  endif()
+
+  set_target_properties(${SWIG_PY_LIB}
+    PROPERTIES
+      OUTPUT_NAME ${SWIG_PY_LIB_OUTPUT}
+  )
+
+  # Suppress warnings on SWIG-generated files
+  target_compile_options(${SWIG_PY_LIB} PRIVATE
+    $<$<CXX_COMPILER_ID:GNU>:-Wno-pedantic -Wno-shadow -Wno-maybe-uninitialized -Wno-unused-parameter -Wno-cast-function-type -Wno-missing-field-initializers>
+    $<$<CXX_COMPILER_ID:Clang>:-Wno-shadow -Wno-maybe-uninitialized -Wno-unused-parameter -Wno-cast-function-type -Wno-missing-field-initializers>
+    $<$<CXX_COMPILER_ID:AppleClang>:-Wno-shadow -Wno-maybe-uninitialized -Wno-unused-parameter -Wno-cast-function-type -Wno-missing-field-initializers>
+  )
+  install(TARGETS ${SWIG_PY_LIB} DESTINATION ${IGN_LIB_INSTALL_DIR}/python/ignition)
+  install(FILES ${CMAKE_BINARY_DIR}/lib/python/math.py DESTINATION ${IGN_LIB_INSTALL_DIR}/python/ignition)
+
+  if (BUILD_TESTING)
+    # Add the Python tests
+    set(python_tests
+      Angle_TEST
+      GaussMarkovProcess_TEST
+      python_TEST
+      Rand_TEST
+      Vector2_TEST
+      Vector3_TEST
+      Vector4_TEST
+    )
+
+    foreach (test ${python_tests})
+      add_test(NAME ${test}.py COMMAND
+        python3 ${CMAKE_SOURCE_DIR}/src/${test}.py)
+
+      set(_env_vars)
+      list(APPEND _env_vars "PYTHONPATH=${FAKE_INSTALL_PREFIX}/lib/python/")
+      list(APPEND _env_vars "LD_LIBRARY_PATH=${FAKE_INSTALL_PREFIX}/lib:$ENV{LD_LIBRARY_PATH}")
+      set_tests_properties(${test}.py PROPERTIES
+        ENVIRONMENT "${_env_vars}")
+    endforeach()
+  endif()
+
+endif()
diff --git a/src/GaussMarkovProcess_TEST.py b/src/GaussMarkovProcess_TEST.py
new file mode 100644
index 000000000..965961329
--- /dev/null
+++ b/src/GaussMarkovProcess_TEST.py
@@ -0,0 +1,58 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+from ignition.math import GaussMarkovProcess
+
+
+class TestGaussMarkovProcess(unittest.TestCase):
+
+    def test_default_constructor(self):
+        gmp = GaussMarkovProcess()
+        self.assertAlmostEqual(0.0, gmp.Start())
+        self.assertAlmostEqual(0.0, gmp.Value())
+        self.assertAlmostEqual(0.0, gmp.Theta())
+        self.assertAlmostEqual(0.0, gmp.Mu())
+        self.assertAlmostEqual(0.0, gmp.Sigma())
+
+    def test_no_noise(self):
+        # Start value of -1.2
+        # Theta (rate at which the process should approach the mean) of 1.0
+        # Mu (mean value) 2.5
+        # Sigma (volatility) of 0.0
+        gmp = GaussMarkovProcess(-1.2, 1.0, 2.5, 0)
+        self.assertAlmostEqual(-1.2, gmp.Start())
+        self.assertAlmostEqual(-1.2, gmp.Value())
+        self.assertAlmostEqual(1.0, gmp.Theta())
+        self.assertAlmostEqual(2.5, gmp.Mu())
+        self.assertAlmostEqual(0.0, gmp.Sigma())
+
+        # This process should steadily increase to the mean value of 2.5 since
+        # there is no noise.
+        for i in range(200):
+            value = gmp.Update(0.1)
+            self.assertGreater(value, -1.2)
+
+        self.assertAlmostEqual(2.5, value, delta=1e-4)
+
+        gmp.Reset()
+        for i in range(200):
+            value = gmp.Update(0.1)
+            self.assertGreater(value, -1.2)
+
+        self.assertAlmostEqual(2.5, value, delta=1e-4)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/src/Rand_TEST.py b/src/Rand_TEST.py
new file mode 100644
index 000000000..2b2d5d60b
--- /dev/null
+++ b/src/Rand_TEST.py
@@ -0,0 +1,47 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+from ignition.math import Rand
+
+class TestRand(unittest.TestCase):
+
+  def test_rand(self):
+    d = Rand.DblUniform(1, 2)
+    self.assertGreaterEqual(d, 1)
+    self.assertLessEqual(d, 2)
+
+    i = Rand.IntUniform(1, 2)
+    self.assertGreaterEqual(i, 1)
+    self.assertLessEqual(i, 2)
+
+  def test_set_seed(self):
+    N = 10
+    first = []
+    second = []
+
+    for i in range(N):
+      Rand.Seed(i)
+      first.append(Rand.IntUniform(-10, 10))
+      second.append(Rand.IntUniform(-10, 10))
+
+    for i in range(N):
+      Rand.Seed(i)
+      self.assertEqual(Rand.Seed(), i)
+      self.assertEqual(first[i], Rand.IntUniform(-10, 10))
+      self.assertEqual(second[i], Rand.IntUniform(-10, 10))
+
+
+if __name__ == '__main__':
+  unittest.main()
diff --git a/src/Vector2_TEST.py b/src/Vector2_TEST.py
new file mode 100644
index 000000000..9ac6d667c
--- /dev/null
+++ b/src/Vector2_TEST.py
@@ -0,0 +1,300 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+import math
+from ignition.math import Vector2d
+
+
+class TestVector2(unittest.TestCase):
+
+    def test_construction(self):
+        v = Vector2d()
+        self.assertAlmostEqual(0.0, v.X())
+        self.assertAlmostEqual(0.0, v.Y())
+
+        vec = Vector2d(1, 0)
+        self.assertEqual(vec.X(), 1)
+        self.assertEqual(vec.Y(), 0)
+
+        vec2 = Vector2d(vec)
+        self.assertEqual(vec2, vec)
+
+        # Copy
+        vec3 = vec
+        self.assertEqual(vec3, vec)
+
+        # Inequality
+        vec4 = Vector2d()
+        self.assertNotEqual(vec, vec4)
+
+    def test_vector2(self):
+        v = Vector2d(1, 2)
+
+        # Distance
+        self.assertAlmostEqual(2.236, v.Distance(Vector2d()), delta=1e-2)
+
+        # Normalize
+        v.Normalize()
+        self.assertTrue(v.Equal(Vector2d(0.447214, 0.894427), 1e-4))
+
+        # Set
+        v.Set(4, 5)
+        self.assertTrue(v.Equal(Vector2d(4, 5), 1e-4))
+
+        # Abs
+        v.Set(-1, -2)
+        self.assertTrue(v.Abs().Equal(Vector2d(1, 2), 1e-4))
+
+        # _eq_
+        v = Vector2d(6, 7)
+        self.assertTrue(v.Equal(Vector2d(6, 7), 1e-4))
+
+        # _add_
+        v = v + Vector2d(1, 2)
+        self.assertTrue(v.Equal(Vector2d(7, 9), 1e-4))
+
+        v += Vector2d(5, 6)
+        self.assertTrue(v.Equal(Vector2d(12, 15), 1e-4))
+
+        # __sub__
+        v = v - Vector2d(2, 4)
+        self.assertTrue(v.Equal(Vector2d(10, 11), 1e-4))
+
+        v.Set(2, 4)
+        v -= Vector2d(1, 6)
+        self.assertTrue(v.Equal(Vector2d(1, -2), 1e-4))
+
+        # __truediv__
+        v.Set(10, 6)
+        v = v / Vector2d(2, 3)
+        self.assertTrue(v.Equal(Vector2d(5, 2), 1e-4))
+
+        v.Set(10, 6)
+        v /= Vector2d(2, 3)
+        self.assertTrue(v.Equal(Vector2d(5, 2), 1e-4))
+
+        # __truediv__ int
+        v.Set(10, 6)
+        v = v / 2
+        self.assertTrue(v.Equal(Vector2d(5, 3), 1e-4))
+
+        v.Set(10, 6)
+        v /= 2
+        self.assertTrue(v.Equal(Vector2d(5, 3), 1e-4))
+
+        # __mul__
+        v.Set(10, 6)
+        v = v * Vector2d(2, 4)
+        self.assertTrue(v.Equal(Vector2d(20, 24), 1e-4))
+
+        v.Set(10, 6)
+        v *= Vector2d(2, 4)
+        self.assertTrue(v.Equal(Vector2d(20, 24), 1e-4))
+
+        # __mul__ int
+        v.Set(10, 6)
+        v = v * 2
+        self.assertTrue(v.Equal(Vector2d(20, 12), 1e-4))
+
+        v.Set(10, 6)
+        v *= 2
+        self.assertTrue(v.Equal(Vector2d(20, 12), 1e-4))
+
+        # IsFinite
+        self.assertTrue(v.IsFinite())
+
+    def test_max(self):
+        vec1 = Vector2d(0.1, 0.2)
+        vec2 = Vector2d(0.3, 0.5)
+        vec3 = Vector2d(0.4, 0.2)
+
+        self.assertAlmostEqual(vec1.Max(), 0.2)
+        self.assertAlmostEqual(vec3.Max(), 0.4)
+
+        vec1.Max(vec2)
+        self.assertAlmostEqual(vec1, Vector2d(0.3, 0.5))
+
+        vec1.Max(vec3)
+        self.assertAlmostEqual(vec1, Vector2d(0.4, 0.5))
+
+    def test_min(self):
+        vec1 = Vector2d(0.3, 0.5)
+        vec2 = Vector2d(0.1, 0.2)
+        vec3 = Vector2d(0.05, 0.1)
+
+        self.assertAlmostEqual(vec1.Min(), 0.3)
+        self.assertAlmostEqual(vec3.Min(), 0.05)
+
+        vec1.Min(vec2)
+        self.assertAlmostEqual(vec1, Vector2d(0.1, 0.2))
+
+        vec1.Min(vec3)
+        self.assertAlmostEqual(vec1, Vector2d(0.05, 0.1))
+
+    def test_equal_tolerance(self):
+        # Test Equal function with specified tolerance
+        self.assertFalse(Vector2d.Zero.Equal(Vector2d.One, 1e-6))
+        self.assertFalse(Vector2d.Zero.Equal(Vector2d.One, 1e-3))
+        self.assertFalse(Vector2d.Zero.Equal(Vector2d.One, 1e-1))
+        self.assertTrue(Vector2d.Zero.Equal(Vector2d.One, 1))
+        self.assertTrue(Vector2d.Zero.Equal(Vector2d.One, 1.1))
+
+    def test_dot(self):
+        v = Vector2d(1, 2)
+        self.assertAlmostEqual(v.Dot(Vector2d(3, 4)), 11.0)
+        self.assertAlmostEqual(v.Dot(Vector2d(0, 0)), 0.0)
+        self.assertAlmostEqual(v.Dot(Vector2d(1, 0)), 1.0)
+        self.assertAlmostEqual(v.Dot(Vector2d(0, 1)), 2.0)
+
+    def test_correct(self):
+        vec1 = Vector2d(0, float("nan"))
+        vec2 = Vector2d(float("inf"), -1)
+        vec3 = Vector2d(10, -2)
+
+        vec1.Correct()
+        vec2.Correct()
+        vec3.Correct()
+
+        self.assertAlmostEqual(vec1, Vector2d(0, 0))
+        self.assertAlmostEqual(vec2, Vector2d(0, -1))
+        self.assertAlmostEqual(vec3, Vector2d(10, -2))
+
+    def test_abs_dot(self):
+        v = Vector2d(1, -2)
+
+        self.assertAlmostEqual(v.AbsDot(Vector2d(3, 4)), 11.0)
+        self.assertAlmostEqual(v.AbsDot(Vector2d(0, 0)), 0.0)
+        self.assertAlmostEqual(v.AbsDot(Vector2d(1, 0)), 1.0)
+        self.assertAlmostEqual(v.AbsDot(Vector2d(0, 1)), 2.0)
+
+    def test_add(self):
+        vec1 = Vector2d(0.1, 0.2)
+        vec2 = Vector2d(1.1, 2.2)
+
+        vec3 = vec1
+        vec3 += vec2
+
+        self.assertAlmostEqual(vec1 + vec2, Vector2d(1.2, 2.4))
+        self.assertAlmostEqual(vec3, Vector2d(1.2, 2.4))
+
+        # Add zero
+        # Scalar right
+        self.assertEqual(vec1 + 0, vec1)
+
+        # Vector left and right
+        self.assertAlmostEqual(Vector2d.Zero + vec1, vec1)
+        self.assertAlmostEqual(vec1 + Vector2d.Zero, vec1)
+
+        # Addition assigment
+        vec4 = Vector2d(vec1)
+        vec4 += 0
+        self.assertEqual(vec4, vec1)
+        vec4 += Vector2d.Zero
+        self.assertAlmostEqual(vec4, vec1)
+
+        # Add non-trivial scalar values left and right
+        self.assertEqual(vec1 + 2.5, Vector2d(2.6, 2.7))
+
+        vec1 = vec4
+        vec4 += 2.5
+        self.assertEqual(vec4, Vector2d(2.6, 2.7))
+
+    def test_sub(self):
+        vec1 = Vector2d(0.1, 0.2)
+        vec2 = Vector2d(1.1, 2.2)
+
+        vec3 = vec2
+        vec3 -= vec1
+
+        self.assertAlmostEqual(vec2 - vec1, Vector2d(1.0, 2.0))
+        self.assertAlmostEqual(vec3, Vector2d(1.0, 2.0))
+
+        # Subtraction with zeros
+        # Scalar right
+        self.assertEqual(vec1 - 0, vec1)
+
+        # Vector left and right
+        self.assertAlmostEqual(Vector2d.Zero - vec1, -vec1)
+        self.assertAlmostEqual(vec1 - Vector2d.Zero, vec1)
+
+        # Subtraction assignment
+        vec4 = Vector2d(vec1)
+        vec4 -= 0
+        self.assertEqual(vec4, vec1)
+        vec4 -= Vector2d.Zero
+        self.assertAlmostEqual(vec4, vec1)
+
+        # Subtract non-trivial scalar values left and right
+        self.assertEqual(vec1 - 2.5, -Vector2d(2.4, 2.3))
+
+        vec4 = vec1
+        vec4 -= 2.5
+        self.assertEqual(vec4, -Vector2d(2.4, 2.3))
+
+    def test_multiply(self):
+        v = Vector2d(0.1, -4.2)
+
+        vec2 = v * 2.0
+        self.assertEqual(vec2, Vector2d(0.2, -8.4))
+
+        vec2 *= 4.0
+        self.assertEqual(vec2, Vector2d(0.8, -33.6))
+
+        # Multiply by zero
+        # Scalar right
+        self.assertEqual(v * 0, Vector2d.Zero)
+
+        # Element-wise vector multiplication
+        self.assertEqual(v * Vector2d.Zero, Vector2d.Zero)
+
+        # Multiply by one
+        # Scalar right
+        self.assertEqual(v * 1, v)
+
+        # Element-wise vector multiplication
+        self.assertEqual(v * Vector2d.One, v)
+
+        # Multiply by non-trivial scalar value
+        scalar = 2.5
+        expect = Vector2d(0.25, -10.5)
+        self.assertEqual(v * scalar, expect)
+
+        # Multiply by itself element-wise
+        v.Set(0.1, 0.5)
+        self.assertAlmostEqual(v * v, Vector2d(0.01, 0.25))
+
+    def test_lenght(self):
+        # Zero vector
+        self.assertAlmostEqual(Vector2d.Zero.Length(), 0.0)
+        self.assertAlmostEqual(Vector2d.Zero.SquaredLength(), 0.0)
+
+        # One vector
+        self.assertAlmostEqual(Vector2d.One.Length(),
+                               math.sqrt(2), delta=1e-10)
+        self.assertAlmostEqual(Vector2d.One.SquaredLength(), 2.0)
+
+        # Arbitrary vector
+        v = Vector2d(0.1, -4.2)
+        self.assertAlmostEqual(v.Length(), 4.20119030752, delta=1e-10)
+        self.assertAlmostEqual(v.SquaredLength(), 17.65)
+
+        # Integer vector
+        v = Vector2d(3, 4)
+        self.assertAlmostEqual(v.Length(), 5)
+        self.assertAlmostEqual(v.SquaredLength(), 25)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/src/Vector3_TEST.py b/src/Vector3_TEST.py
new file mode 100644
index 000000000..a7f1a84a5
--- /dev/null
+++ b/src/Vector3_TEST.py
@@ -0,0 +1,355 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+import math
+from ignition.math import Vector3d
+
+
+class TestVector3(unittest.TestCase):
+
+    def test_construction(self):
+        v = Vector3d()
+        self.assertAlmostEqual(0.0, v.X())
+        self.assertAlmostEqual(0.0, v.Y())
+        self.assertAlmostEqual(0.0, v.Z())
+
+        vec = Vector3d(1, 0, 0)
+        self.assertEqual(vec.X(), 1)
+        self.assertEqual(vec.Y(), 0)
+        self.assertEqual(vec.Z(), 0)
+
+        vec2 = Vector3d(vec)
+        self.assertEqual(vec2, vec)
+
+        # Copy
+        vec3 = vec
+        self.assertEqual(vec3, vec)
+
+        # Inequality
+        vec4 = Vector3d()
+        self.assertNotEqual(vec, vec4)
+
+    def test_vector3(self):
+        # Distance, Length()
+        v = Vector3d(1, 2, 3)
+        self.assertEqual(v.Length(), v.Distance(Vector3d.Zero))
+
+        # Rounded
+        v.Set(1.23, 2.34, 3.55)
+        self.assertEqual(v.Rounded(), Vector3d(1, 2, 4))
+
+        # Round
+        v.Round()
+        self.assertEqual(v.Round(), Vector3d(1, 2, 4))
+
+        # DotProd
+        self.assertEqual(v.Dot(Vector3d(1, 2, 3)), 17.0)
+
+        # DistToLine
+        v.Set(0, 0, 0)
+        self.assertEqual(1.0, v.DistToLine(Vector3d(1, -1, 0),
+                         Vector3d(1, 1, 0)))
+
+        # __truediv__
+        v.Set(4, 4, 4)
+        v = v / Vector3d(1, 2, 4)
+        self.assertEqual(v, Vector3d(4, 2, 1))
+
+        # __truediv__ int
+        v = v / 2
+        self.assertEqual(v, Vector3d(2, 1, 0.5))
+
+        # __mul__
+        v = v * Vector3d(2, 3, 4)
+        self.assertEqual(v, Vector3d(4, 3, 2))
+
+        v.Set(1.23, 2.35, 3.654321)
+        v.Round(1)
+        self.assertEqual(v, Vector3d(1.2, 2.4, 3.7))
+
+        # Abs
+        v.Set(-1, -2, -3)
+        self.assertEqual(v.Abs(), Vector3d(1, 2, 3))
+
+        # __truediv__
+        v.Set(1, 2, 4)
+        v /= Vector3d(1, 4, 4)
+        self.assertEqual(v, Vector3d(1, 0.5, 1))
+
+        # __mul__
+        v.Set(1, 2, 4)
+        v *= Vector3d(2, 0.5, 0.1)
+        self.assertTrue(v.Equal(Vector3d(2, 1, 0.4)))
+
+        # Test the static defines.
+        self.assertEqual(Vector3d.Zero, Vector3d(0, 0, 0))
+
+        self.assertEqual(Vector3d.One, Vector3d(1, 1, 1))
+
+        self.assertEqual(Vector3d.UnitX, Vector3d(1, 0, 0))
+
+        self.assertEqual(Vector3d.UnitY, Vector3d(0, 1, 0))
+
+        self.assertEqual(Vector3d.UnitZ, Vector3d(0, 0, 1))
+
+    def test_distance(self):
+        vec1 = Vector3d(0, 0, 0)
+        vec2 = Vector3d(1, 2, 3)
+
+        dist = vec1.Distance(vec2)
+        self.assertTrue(abs(dist - 3.74165738677) < 1e-6)
+
+        dist2 = vec1.Distance(1, 2, 3)
+        self.assertEqual(dist, dist2)
+
+    def test_sum(self):
+        vec1 = Vector3d(0, 0, 0)
+        vec2 = Vector3d(1, 2, 3)
+
+        sum1 = vec1.Sum()
+        sum2 = vec2.Sum()
+
+        self.assertEqual(sum1, 0)
+        self.assertEqual(sum2, 6)
+
+    def test_squared_length(self):
+        vec1 = Vector3d(0, 0, 0)
+        vec2 = Vector3d(1, 2, 3)
+
+        sum1 = vec1.SquaredLength()
+        sum2 = vec2.SquaredLength()
+
+        self.assertEqual(sum1, 0)
+        self.assertEqual(sum2, 14)
+
+    def test_length(self):
+        # Zero vector
+        self.assertEqual(Vector3d.Zero.Length(), 0.0)
+        self.assertEqual(Vector3d.Zero.SquaredLength(), 0.0)
+
+        # UnitXYZ vectorsIgnition::
+        self.assertEqual(Vector3d.UnitX.Length(), 1.0)
+        self.assertEqual(Vector3d.UnitY.Length(), 1.0)
+        self.assertEqual(Vector3d.UnitZ.Length(), 1.0)
+        self.assertEqual(Vector3d.UnitX.SquaredLength(), 1.0)
+        self.assertEqual(Vector3d.UnitY.SquaredLength(), 1.0)
+        self.assertEqual(Vector3d.UnitZ.SquaredLength(), 1.0)
+
+        # One vector
+        self.assertTrue(Vector3d.One.Length() -
+                        abs(math.sqrt(3.0)) < 1e-10)
+
+        self.assertEqual(Vector3d.One.SquaredLength(), 3.0)
+
+        # Arbitrary vector
+        v = Vector3d(0.1, -4.2, 2.5)
+        self.assertTrue(abs(v.Length() - 4.88876262463) < 1e-10)
+
+        self.assertTrue(abs(v.SquaredLength() - 23.9) < 1e-10)
+
+    def test_normalize(self):
+        vec1 = Vector3d(0, 0, 0)
+        vec2 = Vector3d(1, 2, 3)
+
+        vec3 = vec1.Normalize()
+        self.assertEqual(vec3, vec1)
+        self.assertEqual(vec1, Vector3d.Zero)
+
+        vec3 = vec2.Normalize()
+        self.assertEqual(vec3, vec2)
+        self.assertEqual(vec2, Vector3d(0.267261, 0.534522, 0.801784))
+
+        vecConst = Vector3d(1, 2, 3)
+        self.assertEqual(vecConst.Normalized(), vec3)
+        self.assertEqual(vecConst, Vector3d(1, 2, 3))
+
+    def test_ge_normal(self):
+        vec1 = Vector3d(0, 0, 0)
+        vec2 = Vector3d(0, 1, 0)
+        vec3 = Vector3d(1, 1, 0)
+
+        norm = Vector3d.Normal(vec1, vec2, vec3)
+        self.assertEqual(norm, Vector3d(0, 0, -1))
+
+    def test_perpendicular(self):
+        vec1 = Vector3d(1, 1, 0)
+        vec2 = Vector3d(0, 1, 1)
+        vec3 = Vector3d(1e-7, 1e-7, 1e-7)
+        vec4 = Vector3d(1, 0, 0)
+
+        self.assertEqual(vec1.Perpendicular(), Vector3d(0, 0, -1))
+        self.assertEqual(vec2.Perpendicular(), Vector3d(0, 1, -1))
+        self.assertEqual(vec3.Perpendicular(), Vector3d(0, 0, 0))
+        self.assertEqual(vec4.Perpendicular(), Vector3d(0, 0, 1))
+
+    def test_max(self):
+        vec1 = Vector3d(0.1, 0.2, 0.3)
+        vec2 = Vector3d(0.2, 0.3, 0.4)
+        vec3 = Vector3d(0.1, 0.2, 0.3)
+
+        self.assertTrue(abs(vec1.Max() - 0.3) < 1e-10)
+
+        vec1.Max(vec2)
+        self.assertEqual(vec1, Vector3d(0.2, 0.3, 0.4))
+
+        vec1.Max(vec3)
+        self.assertEqual(vec1, Vector3d(0.2, 0.3, 0.4))
+
+    def test_min(self):
+        vec1 = Vector3d(0.1, 0.2, 0.3)
+        vec2 = Vector3d(0.2, 0.3, 0.4)
+        vec3 = Vector3d(0.05, 0.1, 0.2)
+
+        self.assertTrue(abs(vec1.Min() - 0.1) < 1e-10)
+
+        vec1.Min(vec2)
+        self.assertEqual(vec1, Vector3d(0.1, 0.2, 0.3))
+
+        vec1.Min(vec3)
+        self.assertEqual(vec1, Vector3d(0.05, 0.1, 0.2))
+
+    def test_add(self):
+        vec1 = Vector3d(0.1, 0.2, 0.4)
+        vec2 = Vector3d(1.1, 2.2, 3.4)
+
+        vec3 = vec1
+        vec3 += vec2
+
+        self.assertEqual(vec1 + vec2, Vector3d(1.2, 2.4, 3.8))
+        self.assertEqual(vec3, Vector3d(1.2, 2.4, 3.8))
+
+        # Add zeros
+        # Scalar right
+        self.assertEqual(vec1 + 0, vec1)
+
+        # Vector left and right
+        self.assertEqual(Vector3d.Zero + vec1, vec1)
+        self.assertEqual(vec1 + Vector3d.Zero, vec1)
+
+        # Addition assignment
+        vec4 = vec1
+        vec4 += 0
+        self.assertEqual(vec4, vec1)
+        vec4 += Vector3d.Zero
+        self.assertEqual(vec4, vec1)
+
+        # Add non-trivial scalar values left and right
+        self.assertEqual(vec1 + 2.5, Vector3d(2.6, 2.7, 2.9))
+
+        vec1 = vec4
+        vec4 += 2.5
+        self.assertEqual(vec4, Vector3d(2.6, 2.7, 2.9))
+
+    def test_sub(self):
+        vec1 = Vector3d(0.1, 0.2, 0.4)
+        vec2 = Vector3d(1.1, 2.2, 3.4)
+
+        vec3 = vec2
+        vec3 -= vec1
+
+        self.assertEqual(vec2 - vec1, Vector3d(1.0, 2.0, 3.0))
+        self.assertEqual(vec3, Vector3d(1.0, 2.0, 3.0))
+
+        # Subtraction with zeros
+        # Scalar right
+        self.assertEqual(vec1 - 0, vec1)
+
+        # Vector left and right
+        self.assertEqual(Vector3d.Zero - vec1, -vec1)
+        self.assertEqual(vec1 - Vector3d.Zero, vec1)
+
+        # Subtraction assignment
+        vec4 = vec1
+        vec4 -= 0
+        self.assertEqual(vec4, vec1)
+        vec4 -= Vector3d.Zero
+        self.assertEqual(vec4, vec1)
+
+        # Subtract non-trivial scalar values left and right
+        self.assertEqual(vec1 - 2.5, -Vector3d(2.4, 2.3, 2.1))
+
+        vec4 = vec1
+        vec4 -= 2.5
+        self.assertEqual(vec4, -Vector3d(2.4, 2.3, 2.1))
+
+    def test_divide(self):
+        vec1 = Vector3d(0.1, 0.2, 0.4)
+
+        vec3 = vec1 / 2.0
+        self.assertEqual(vec3, Vector3d(0.05, 0.1, 0.2))
+
+        vec3 /= 4.0
+        self.assertEqual(vec3, Vector3d(0.0125, 0.025, 0.05))
+
+    def test_multiply(self):
+        v = Vector3d(0.1, 0.2, 0.3)
+
+        vec2 = v * 2.0
+        self.assertEqual(vec2, Vector3d(0.2, 0.4, 0.6))
+
+        vec2 *= 4.0
+        self.assertEqual(vec2, Vector3d(0.8, 1.6, 2.4))
+
+        # Multiply by zero
+        # Scalar right
+        self.assertEqual(v * 0, Vector3d.Zero)
+
+        # Element-wise vector multiplication
+        self.assertEqual(v * Vector3d.Zero, Vector3d.Zero)
+
+        # Multiply by one
+        # Scalar right
+        self.assertEqual(v * 1, v)
+
+        # Element-wise vector multiplication
+        self.assertEqual(v * Vector3d.One, v)
+
+        # Multiply by non-trivial scalar value
+        scalar = 2.5
+        expect = Vector3d(0.25, 0.5, 0.75)
+        self.assertEqual(v * scalar, expect)
+
+        # Multiply by itself element-wise
+        self.assertEqual(v*v, Vector3d(0.01, 0.04, 0.09))
+
+    def test_not_equal(self):
+        vec1 = Vector3d(0.1, 0.2, 0.3)
+        vec2 = Vector3d(0.2, 0.2, 0.3)
+        vec3 = Vector3d(0.1, 0.2, 0.3)
+
+        self.assertTrue(vec1 != vec2)
+        self.assertTrue(not(vec1 != vec3))
+
+    def test_equal(self):
+        self.assertTrue(not Vector3d.Zero.Equal(
+          Vector3d.One, 1e-6))
+        self.assertTrue(not Vector3d.Zero.Equal(
+          Vector3d.One, 1e-3))
+        self.assertTrue(not Vector3d.Zero.Equal(
+          Vector3d.One, 1e-1))
+
+        self.assertTrue(Vector3d.Zero.Equal(
+            Vector3d.One, 1))
+        self.assertTrue(Vector3d.Zero.Equal(
+            Vector3d.One, 1.1))
+
+    def test_finite(self):
+        vec1 = Vector3d(0.1, 0.2, 0.3)
+
+        self.assertTrue(vec1.IsFinite())
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/src/Vector4_TEST.py b/src/Vector4_TEST.py
new file mode 100644
index 000000000..ca4a6fa60
--- /dev/null
+++ b/src/Vector4_TEST.py
@@ -0,0 +1,261 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+
+import unittest
+import math
+from ignition.math import Vector4d
+
+
+class TestVector4(unittest.TestCase):
+
+    def test_construction(self):
+        v = Vector4d()
+        self.assertAlmostEqual(0.0, v.X())
+        self.assertAlmostEqual(0.0, v.Y())
+        self.assertAlmostEqual(0.0, v.Z())
+        self.assertAlmostEqual(0.0, v.W())
+
+        vec = Vector4d(1, 0, 0, 0)
+        self.assertEqual(vec.X(), 1)
+        self.assertEqual(vec.Y(), 0)
+        self.assertEqual(vec.Z(), 0)
+        self.assertEqual(vec.W(), 0)
+
+        vec2 = Vector4d(vec)
+        self.assertEqual(vec2, vec)
+
+        # Copy
+        vec3 = vec
+        self.assertEqual(vec3, vec)
+
+        # Inequality
+        vec4 = Vector4d()
+        self.assertNotEqual(vec, vec4)
+
+    def test_vector4(self):
+        v = Vector4d()
+        # Distance, Length()
+        v.Set(1, 2, 3, 4)
+        self.assertEqual(v.Length(), v.Distance(Vector4d.Zero))
+
+        # __truediv__
+        v.Set(4, 4, 4, 4)
+        v = v / Vector4d(1, 2, 2, 4)
+        self.assertEqual(v, Vector4d(4, 2, 2, 1))
+
+        # __truediv__ int
+        v = v / 2
+        self.assertEqual(v, Vector4d(2, 1, 1, 0.5))
+
+        # __mul__
+        v = v * Vector4d(2, 3, 3, 4)
+        self.assertEqual(v, Vector4d(4, 3, 3, 2))
+
+        # __truediv__
+        v.Set(1, 2, 2, 4)
+        v /= Vector4d(1, 4, 8, 4)
+        self.assertEqual(v, Vector4d(1, 0.5, 0.25, 1))
+
+        # __mul__
+        v.Set(1, 2, 2, 4)
+        v *= Vector4d(2, 0.5, 0.25, 0.1)
+        self.assertEqual(v, Vector4d(2, 1, 0.5, 0.4))
+
+        # Test the static defines.
+        self.assertEqual(Vector4d.Zero,
+                         Vector4d(0, 0, 0, 0))
+
+        self.assertEqual(Vector4d.One,
+                         Vector4d(1, 1, 1, 1))
+
+    def test_distance(self):
+        vec1 = Vector4d(0, 0, 0, 0)
+        vec2 = Vector4d(1, 2, 3, 4)
+
+        dist = vec1.Distance(vec2)
+        self.assertTrue(abs(dist - 5.47722557505) < 1e-6)
+
+    def test_squared_length(self):
+        vec1 = Vector4d(0, 0, 0, 0)
+        vec2 = Vector4d(1, 2, 3, 4)
+
+        sum1 = vec1.SquaredLength()
+        sum2 = vec2.SquaredLength()
+
+        self.assertEqual(sum1, 0)
+        self.assertEqual(sum2, 30)
+
+    def test_length(self):
+        # Zero vector
+        self.assertEqual(Vector4d.Zero.Length(), 0.0)
+        self.assertEqual(Vector4d.Zero.SquaredLength(), 0.0)
+
+        # One vector
+        self.assertTrue(abs(Vector4d.One.Length() - math.sqrt(4.0)) < 1e-10)
+
+        self.assertEqual(Vector4d.One.SquaredLength(), 4.0)
+
+        # Arbitrary vector
+        v = Vector4d(0.1, -4.2, 2.5, -1.2)
+        self.assertTrue(abs(v.Length() - 5.03388517946) < 1e-10)
+
+        self.assertTrue(abs(v.SquaredLength() - 25.34) < 1e-10)
+
+    def test_normalize(self):
+        vec1 = Vector4d(0, 0, 0, 0)
+        vec2 = Vector4d(1, 2, 3, 4)
+
+        vec3 = vec1
+        vec3.Normalize()
+        self.assertEqual(vec3, vec1)
+        self.assertEqual(vec1, Vector4d.Zero)
+
+        vec3 = vec2
+        vec2.Normalize()
+        self.assertTrue(vec2.Equal(Vector4d(0.182575,
+                        0.365150, 0.547725, 0.730300), 1e-5))
+
+    def test_add(self):
+        vec1 = Vector4d(0.1, 0.2, 0.4, 0.8)
+        vec2 = Vector4d(1.1, 2.2, 3.4, 4.3)
+
+        vec3 = vec1
+        vec3 += vec2
+
+        self.assertEqual(vec1 + vec2, Vector4d(1.2, 2.4, 3.8, 5.1))
+        self.assertEqual(vec3, Vector4d(1.2, 2.4, 3.8, 5.1))
+
+        # Addition with zeros
+        # Scalar right
+        self.assertEqual(vec1 + 0, vec1)
+
+        # Vector left and right
+        self.assertEqual(Vector4d.Zero + vec1, vec1)
+        self.assertEqual(vec1 + Vector4d.Zero, vec1)
+
+        # Addition assignment
+        vec4 = vec1
+        vec4 += 0
+        self.assertEqual(vec4, vec1)
+        vec4 += Vector4d.Zero
+        self.assertEqual(vec4, vec1)
+
+        # Add non-trivial scalar values left and right
+        self.assertEqual(vec1 + 2.5, Vector4d(2.6, 2.7, 2.9, 3.3))
+
+        vec1 = vec4
+        vec4 += 2.5
+        self.assertEqual(vec4, Vector4d(2.6, 2.7, 2.9, 3.3))
+
+    def test_sub(self):
+        vec1 = Vector4d(0.1, 0.2, 0.4, 0.8)
+        vec2 = Vector4d(1.1, 2.2, 3.4, 4.3)
+
+        vec3 = vec2
+        vec3 -= vec1
+
+        self.assertEqual(vec2 - vec1, Vector4d(1.0, 2.0, 3.0, 3.5))
+        self.assertEqual(vec3, Vector4d(1.0, 2.0, 3.0, 3.5))
+
+        # Subtraction with zeros
+        # Scalar right
+        self.assertEqual(vec1 - 0, vec1)
+
+        # Vector left and right
+        self.assertEqual(Vector4d.Zero - vec1, -vec1)
+        self.assertEqual(vec1 - Vector4d.Zero, vec1)
+
+        # Subtraction assignment
+        vec4 = vec1
+        vec4 -= 0
+        self.assertEqual(vec4, vec1)
+        vec4 -= Vector4d.Zero
+        self.assertEqual(vec4, vec1)
+
+        # Subtract non-trivial scalar values left and right
+        self.assertEqual(vec1 - 2.5, -Vector4d(2.4, 2.3, 2.1, 1.7))
+
+        vec4 = vec1
+        vec4 -= 2.5
+        self.assertEqual(vec4, -Vector4d(2.4, 2.3, 2.1, 1.7))
+
+    def test_divide(self):
+        vec1 = Vector4d(0.1, 0.2, 0.4, 0.8)
+
+        vec3 = vec1 / 2.0
+        self.assertEqual(vec3, Vector4d(0.05, 0.1, 0.2, 0.4))
+
+        vec3 /= 4.0
+        self.assertEqual(vec3, Vector4d(0.0125, 0.025, 0.05, 0.1))
+
+    def test_multiply(self):
+        v = Vector4d(0.1, 0.2, 0.3, 0.4)
+
+        vec3 = v * 2.0
+        self.assertEqual(vec3, Vector4d(0.2, 0.4, 0.6, 0.8))
+
+        vec3 *= 4.0
+        self.assertEqual(vec3, Vector4d(0.8, 1.6, 2.4, 3.2))
+
+        # Multiply by zero
+        # Scalar right
+        self.assertEqual(v * 0, Vector4d.Zero)
+
+        # Element-wise vector multiplication
+        self.assertEqual(v * Vector4d.Zero, Vector4d.Zero)
+
+        # Multiply by one
+        # Scalar right
+        self.assertEqual(v * 1, v)
+
+        # Element-wise vector multiplication
+        self.assertEqual(v * Vector4d.One, v,)
+
+        # Multiply by non-trivial scalar value
+        scalar = 2.5
+        expect = Vector4d(0.25, 0.5, 0.75, 1.0)
+        self.assertEqual(v * scalar, expect)
+
+        # Multiply by itself element-wise
+        self.assertEqual(v*v, Vector4d(0.01, 0.04, 0.09, 0.16))
+
+    def test_not_equal(self):
+        vec1 = Vector4d(0.1, 0.2, 0.3, 0.4)
+        vec2 = Vector4d(0.2, 0.2, 0.3, 0.4)
+        vec3 = Vector4d(0.1, 0.2, 0.3, 0.4)
+
+        self.assertTrue(vec1 != vec2)
+        self.assertTrue(not(vec1 != vec3))
+
+    def test_equal(self):
+        self.assertTrue(not Vector4d.Zero.Equal(
+                        Vector4d.One, 1e-6))
+        self.assertTrue(not Vector4d.Zero.Equal(
+                        Vector4d.One, 1e-3))
+        self.assertTrue(not Vector4d.Zero.Equal(
+                        Vector4d.One, 1e-1))
+
+        self.assertTrue(Vector4d.Zero.Equal(
+            Vector4d.One, 1))
+        self.assertTrue(Vector4d.Zero.Equal(
+            Vector4d.One, 1.1))
+
+    def test_finite(self):
+        vec1 = Vector4d(0.1, 0.2, 0.3, 0.4)
+        self.assertTrue(vec1.IsFinite())
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/src/ign_math.i b/src/ign_math.i
deleted file mode 100644
index 651220b6d..000000000
--- a/src/ign_math.i
+++ /dev/null
@@ -1,6 +0,0 @@
-%include "Angle.i"
-%include "GaussMarkovProcess.i"
-%include "Rand.i"
-%include "Vector2.i"
-%include "Vector3.i"
-%include "Vector4.i"
diff --git a/src/python/python.i b/src/python/python.i
new file mode 100644
index 000000000..e80c81f87
--- /dev/null
+++ b/src/python/python.i
@@ -0,0 +1,7 @@
+%module "math"
+%include ../Angle.i
+%include ../GaussMarkovProcess.i
+%include ../Rand.i
+%include ../Vector2.i
+%include ../Vector3.i
+%include ../Vector4.i
diff --git a/src/python_TEST.py b/src/python_TEST.py
new file mode 100644
index 000000000..86b6372bb
--- /dev/null
+++ b/src/python_TEST.py
@@ -0,0 +1,32 @@
+# Copyright (C) 2021 Open Source Robotics Foundation
+
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+
+#       http://www.apache.org/licenses/LICENSE-2.0
+
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+
+import unittest
+import ignition.math
+
+
+class TestPythonInterface(unittest.TestCase):
+
+    def test_construction(self):
+        angle1 = ignition.math.Angle()
+        self.assertEqual(angle1.Radian(), 0.0)
+        v1 = ignition.math.Vector3d(0, 0, 0)
+        self.assertEqual(v1, ignition.math.Vector3d.Zero)
+        v2 = ignition.math.Vector2d(1, 2)
+        self.assertEqual(v2.X(), 1)
+        self.assertEqual(v2.Y(), 2)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/src/ruby.i b/src/ruby.i
deleted file mode 100644
index 689be4030..000000000
--- a/src/ruby.i
+++ /dev/null
@@ -1,2 +0,0 @@
-%module "ignition::math"
-%include ign_math.i
diff --git a/src/ruby/ruby.i b/src/ruby/ruby.i
new file mode 100644
index 000000000..0364d1d74
--- /dev/null
+++ b/src/ruby/ruby.i
@@ -0,0 +1,7 @@
+%module "ignition::math"
+%include ../Angle.i
+%include ../GaussMarkovProcess.i
+%include ../Rand.i
+%include ../Vector2.i
+%include ../Vector3.i
+%include ../Vector4.i
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f59df24af..035530294 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -23,7 +23,7 @@ include_directories(${GTEST_INCLUDE_DIRS})
 #============================================================================
 # Do a fake install of ign-math in order to test the examples.
 #============================================================================
-set(FAKE_INSTALL_PREFIX "${CMAKE_BINARY_DIR}/fake/install")
+# install to FAKE_INSTALL_PREFIX defined in root CMakeLists.txt
 
 file(MAKE_DIRECTORY ${FAKE_INSTALL_PREFIX})