diff --git a/src/route/routing-table-calculator.cpp b/src/route/routing-table-calculator.cpp
index 54d486d..7aa6f48 100644
--- a/src/route/routing-table-calculator.cpp
+++ b/src/route/routing-table-calculator.cpp
@@ -483,37 +483,125 @@
     return UNKNOWN_DISTANCE;
   }
 
-  double srcTheta = srcLsa->getCorTheta();
-  double destTheta = destLsa->getCorTheta();
-
-  double diffTheta = fabs(srcTheta - destTheta);
-
-  if (diffTheta > MATH_PI) {
-    diffTheta = 2 * MATH_PI - diffTheta;
-  }
+  std::vector<double> srcTheta = srcLsa->getCorTheta();
+  std::vector<double> destTheta = destLsa->getCorTheta();
 
   double srcRadius = srcLsa->getCorRadius();
   double destRadius = destLsa->getCorRadius();
 
-  if (srcRadius == UNKNOWN_RADIUS && destRadius == UNKNOWN_RADIUS) {
+  double diffTheta = calculateAngularDistance(srcTheta, destTheta);
+
+  if (srcRadius == UNKNOWN_RADIUS || destRadius == UNKNOWN_RADIUS ||
+      diffTheta == UNKNOWN_DISTANCE) {
     return UNKNOWN_DISTANCE;
   }
 
-  if (diffTheta == 0) {
-    distance = fabs(srcRadius - destRadius);
-  }
-  else {
-    distance = acosh((cosh(srcRadius) * cosh(destRadius)) -
-                     (sinh(srcRadius) * sinh(destRadius) * cos(diffTheta)));
-  }
+  // double r_i, double r_j, double delta_theta, double zeta = 1 (default)
+  distance = calculateHyperbolicDistance(srcRadius, destRadius, diffTheta);
 
   _LOG_TRACE("Distance from " << src << " to " << dest << " is " << distance);
 
   return distance;
 }
 
-void HyperbolicRoutingCalculator::addNextHop(ndn::Name dest, std::string faceUri,
-                                             double cost, RoutingTable& rt)
+double
+HyperbolicRoutingCalculator::calculateAngularDistance(std::vector<double> angleVectorI,
+                                                      std::vector<double> angleVectorJ)
+{
+  // It is not possible for angle vector size to be zero as ensured by conf-file-processor
+
+  // https://en.wikipedia.org/wiki/N-sphere#Spherical_coordinates
+
+  // Check if two vector lengths are the same
+  if (angleVectorI.size() != angleVectorJ.size()) {
+    _LOG_ERROR("Angle vector sizes do not match");
+    return UNKNOWN_DISTANCE;
+  }
+
+  // Check if all angles are within the [0, PI] and [0, 2PI] ranges
+  if (angleVectorI.size() > 1) {
+    for (unsigned int k = 0; k < angleVectorI.size() - 1; k++) {
+      if ((angleVectorI[k] > M_PI && angleVectorI[k] < 0.0) ||
+          (angleVectorJ[k] > M_PI && angleVectorJ[k] < 0.0)) {
+        _LOG_ERROR("Angle outside [0, PI]");
+        return UNKNOWN_DISTANCE;
+      }
+    }
+  }
+  if (angleVectorI[angleVectorI.size()-1] > 2.*M_PI ||
+      angleVectorI[angleVectorI.size()-1] < 0.0) {
+    _LOG_ERROR("Angle not within [0, 2PI]");
+    return UNKNOWN_DISTANCE;
+  }
+
+  if (angleVectorI[angleVectorI.size()-1] > 2.*M_PI ||
+      angleVectorI[angleVectorI.size()-1] < 0.0) {
+    _LOG_ERROR("Angle not within [0, 2PI]");
+    return UNKNOWN_DISTANCE;
+  }
+
+  // deltaTheta = arccos(vectorI . vectorJ) -> do the inner product
+  double innerProduct = 0.0;
+
+  // Calculate x0 of the vectors
+  double x0i = std::cos(angleVectorI[0]);
+  double x0j = std::cos(angleVectorJ[0]);
+
+  // Calculate xn of the vectors
+  double xni = std::sin(angleVectorI[angleVectorI.size() - 1]);
+  double xnj = std::sin(angleVectorJ[angleVectorJ.size() - 1]);
+
+  // Do the aggregation of the (n-1) coordinates (if there is more than one angle)
+  // i.e contraction of all (n-1)-dimensional angular coordinates to one variable
+  for (unsigned int k = 0; k < angleVectorI.size() - 1; k++) {
+    xni *= std::sin(angleVectorI[k]);
+    xnj *= std::sin(angleVectorJ[k]);
+  }
+  innerProduct += (x0i * x0j) + (xni * xnj);
+
+  // If d > 1
+  if (angleVectorI.size() > 1) {
+    for (unsigned int m = 1; m < angleVectorI.size(); m++) {
+      // calculate euclidean coordinates given the angles and assuming R_sphere = 1
+      double xmi = std::cos(angleVectorI[m]);
+      double xmj = std::cos(angleVectorJ[m]);
+      for (unsigned int l = 0; l < m; l++) {
+        xmi *= std::sin(angleVectorI[l]);
+        xmj *= std::sin(angleVectorJ[l]);
+      }
+      innerProduct += xmi * xmj;
+    }
+  }
+
+  // ArcCos of the inner product gives the angular distance
+  // between two points on a d-dimensional sphere
+  return std::acos(innerProduct);
+}
+
+double
+HyperbolicRoutingCalculator::calculateHyperbolicDistance(double rI, double rJ,
+                                                         double deltaTheta)
+{
+  if (deltaTheta == UNKNOWN_DISTANCE) {
+    return UNKNOWN_DISTANCE;
+  }
+
+  // Usually, we set zeta = 1 in all experiments
+  double zeta = 1;
+
+  if (deltaTheta <= 0.0 || rI <= 0.0 || rJ <= 0.0) {
+    _LOG_ERROR("Delta theta or rI or rJ is <= 0");
+    return UNKNOWN_DISTANCE;
+  }
+
+  double xij = (1. / zeta) * std::acosh(std::cosh(zeta*rI) * std::cosh(zeta*rJ) -
+               std::sinh(zeta*rI)*std::sinh(zeta*rJ)*std::cos(deltaTheta));
+  return xij;
+}
+
+void
+HyperbolicRoutingCalculator::addNextHop(ndn::Name dest, std::string faceUri,
+                                        double cost, RoutingTable& rt)
 {
   NextHop hop(faceUri, cost);
   hop.setHyperbolic(true);