diff --git a/Framework/Geometry/Point.h b/Framework/Geometry/Point.h
index 04b9ca41716fea495286b8b4f90efcd646d31c23..b6a6fff21eac45c962cbf044266ce6c0f8fe28bb 100644
--- a/Framework/Geometry/Point.h
+++ b/Framework/Geometry/Point.h
@@ -3,13 +3,13 @@
 
 #include <Geometry/BaseVector.h>
 #include <Geometry/QuantityVector.h>
+#include <Geometry/Vector.h>
 #include <Units/PhysicalUnits.h>
 
 class Point : public BaseVector<phys::units::length_d>
 {
     using Length = phys::units::quantity<phys::units::length_d, double>;
     
-    
 public:
     Point(CoordinateSystem const& pCS, QuantityVector<phys::units::length_d> pQVector) :
         BaseVector<phys::units::length_d>(pCS, pQVector)
@@ -43,6 +43,17 @@ public:
         BaseVector<phys::units::length_d>::qVector = getCoordinates(pCS);
         BaseVector<phys::units::length_d>::cs = &pCS;
     }
+    
+    Point operator+(Vector<phys::units::length_d> const& pVec) const
+    {        
+        return Point(*BaseVector<phys::units::length_d>::cs, getCoordinates() + pVec.getComponents(*BaseVector<phys::units::length_d>::cs));
+    }
+    
+    Vector<phys::units::length_d> operator-(Point const& pB) const
+    {
+        auto& cs = *BaseVector<phys::units::length_d>::cs;
+        return Vector<phys::units::length_d>(cs, getCoordinates() - pB.getCoordinates(cs));
+    }
 };
 
 #endif
diff --git a/Framework/Geometry/QuantityVector.h b/Framework/Geometry/QuantityVector.h
index 94ce89dcab3d7c94d32267ea3e21bf32494ddd14..64a312f3ddbabc4afee627f87b0f09300055bbd0 100644
--- a/Framework/Geometry/QuantityVector.h
+++ b/Framework/Geometry/QuantityVector.h
@@ -11,6 +11,7 @@ class QuantityVector
 {
 protected:
     using Quantity = phys::units::quantity<dim, double>;
+    using QuantitySquared = decltype(Quantity(phys::units::detail::magnitude_tag, 0) * Quantity(phys::units::detail::magnitude_tag, 0));
     
 public:
     Eigen::Vector3d eVector;
@@ -25,10 +26,32 @@ public:
     {
     }
     
-    Quantity operator[](size_t index) const
+    auto operator[](size_t index) const
     {
         return Quantity(phys::units::detail::magnitude_tag, eVector[index]);
     }
+    
+    auto norm() const
+    {
+        return Quantity(phys::units::detail::magnitude_tag, eVector.norm());
+    }
+    
+    auto squaredNorm() const
+    {
+        return QuantitySquared(phys::units::detail::magnitude_tag, eVector.squaredNorm());
+    }
+    
+    auto operator+(QuantityVector<dim> const& pQVec) const
+    {
+        return QuantityVector<dim>(eVector + pQVec.eVector);
+    }
+    
+    auto operator-(QuantityVector<dim> const& pQVec) const
+    {
+        return QuantityVector<dim>(eVector - pQVec.eVector);
+    }
+    
+    //auto operator*(
 };
 
 template <typename dim>
diff --git a/Framework/Geometry/Sphere.h b/Framework/Geometry/Sphere.h
new file mode 100644
index 0000000000000000000000000000000000000000..2a3ecf1af55c3099982504214dc6b328e23f3392
--- /dev/null
+++ b/Framework/Geometry/Sphere.h
@@ -0,0 +1,27 @@
+#ifndef SPHERE_H_
+#define SPHERE_H_
+
+#include "Point.h"
+#include <phys/units/quantity.hpp>
+
+class Sphere
+{
+    using Length = phys::units::quantity<phys::units::length_d, double>;
+    Point center;
+    Length const radius;
+    
+public:
+    Sphere(Point const& pCenter, Length const pRadius) :
+        center(pCenter), radius(pRadius)
+    {
+        
+    }
+        
+    auto isInside(Point const& p) const
+    {
+        return radius*radius > (center - p).squaredNorm();
+    }
+
+};
+
+#endif
diff --git a/Framework/Geometry/Vector.h b/Framework/Geometry/Vector.h
index fc4e171f18443a9b5a9f4bc871aaf7a54af7aa4d..417cb19a52cc5bae2af85c488ad09db766f7dd2d 100644
--- a/Framework/Geometry/Vector.h
+++ b/Framework/Geometry/Vector.h
@@ -10,13 +10,12 @@ class Vector : public BaseVector<dim>
 {
     using Quantity = phys::units::quantity<dim, double>;
     
+public:
     Vector(CoordinateSystem const& pCS, QuantityVector<dim> pQVector) :
-        BaseVector<Quantity>(pCS, pQVector)
+        BaseVector<dim>(pCS, pQVector)
     {
     }
-        
-    
-public:
+
     Vector(CoordinateSystem const& cs, Quantity x, Quantity y, Quantity z) :
         BaseVector<dim>(cs, QuantityVector<dim>(x, y, z))
     {
@@ -47,7 +46,29 @@ public:
     
     auto norm() const
     {
-        return Quantity(BaseVector<dim>::qVector.eVector.norm());
+        return BaseVector<dim>::qVector.norm();
+    }
+    
+    auto squaredNorm() const
+    {
+        return BaseVector<dim>::qVector.squaredNorm();
+    }
+        
+    template <typename dim2>
+    auto parallelProjectionOnto(BaseVector<dim2> const& pVec, CoordinateSystem const& pCS) const
+    {
+        auto const ourCompVec = getComponents(pCS);
+        auto const otherCompVec = pVec.getComponents(pVec);
+        auto const& a = ourCompVec.eVector;
+        auto const& b = otherCompVec.eVector;
+        
+        return Vector<dim>(pCS, (a * b) / b.squaredNorm() * b);
+    }
+    
+    template <typename dim2>
+    auto parallelProjectionOnto(BaseVector<dim2> const& pVec)
+    {
+        return parallelProjectionOnto<dim2>(pVec, *BaseVector<dim>::cs);
     }
     
     //~ template <typename dim2>
diff --git a/Main/geometry_example.cc b/Main/geometry_example.cc
new file mode 100644
index 0000000000000000000000000000000000000000..8c057ed54d68fe9fdb73cf279c7e05c9adbdc47a
--- /dev/null
+++ b/Main/geometry_example.cc
@@ -0,0 +1,40 @@
+/*
+  from within ngc/build/ directory:
+  g++ -O3 -I ../third_party/ -I<path to eigen> -I ../ -Wall -pedantic \
+   -std=c++14 ../Main/geometry_example.cc \
+   ../Framework/Geometry/CoordinateSystem.cc -o geometry_example
+*/
+#include <Framework/Geometry/Vector.h>
+#include <Framework/Geometry/Sphere.h>
+#include <Framework/Geometry/Point.h>
+#include <Framework/Geometry/CoordinateSystem.h>
+#include <phys/units/quantity.hpp>
+#include <phys/units/io.hpp>
+#include <iostream>
+#include <cstdlib>
+
+int main()
+{
+    using namespace phys::units;
+    using namespace phys::units::io;
+    using namespace phys::units::literals;
+    
+    CoordinateSystem root;
+    Point const p1(root, {0_m, 0_m, 0_m});
+    CoordinateSystem cs2 = root.translate({0_m, 0_m, 1_m});
+    Point const p2(cs2, {0_m, 0_m, 0_m});
+    
+    auto const diff = p2 - p1;
+    auto const norm = diff.squaredNorm();
+    
+    std::cout << diff.getComponents() << std::endl;
+    std::cout << norm << std::endl;
+    
+    Sphere s(p1, 10_m);
+    std::cout << s.isInside(p2) << std::endl;    
+    
+    Sphere s2(p1, 3_um);
+    std::cout << s2.isInside(p2) << std::endl;    
+    
+    return EXIT_SUCCESS;
+}