From c5d7479d00925ef8573ec640226866244f74f138 Mon Sep 17 00:00:00 2001
From: Zeno Rogue <zeno@attnam.com>
Date: Sat, 11 Apr 2020 15:01:01 +0200
Subject: [PATCH] nisot:: geodesic_step is now based on RK4 instead of midpoint

---
 nonisotropic.cpp | 22 ++++++++++++----------
 1 file changed, 12 insertions(+), 10 deletions(-)

diff --git a/nonisotropic.cpp b/nonisotropic.cpp
index 3894fd16..a56791dd 100644
--- a/nonisotropic.cpp
+++ b/nonisotropic.cpp
@@ -1970,18 +1970,20 @@ EX namespace nisot {
     #endif
     return true;
     }
+
+  hyperpoint get_acceleration(const hyperpoint& at, const hyperpoint& vel) {
+    return christoffel(at, vel, vel);
+    }
   
-  EX void geodesic_step(hyperpoint& at, hyperpoint& velocity) {
-    auto acc = christoffel(at, velocity, velocity);
+  EX void geodesic_step(hyperpoint& at, hyperpoint& vel) {
+    /* RK4 method */
+    auto acc1 = get_acceleration(at, vel);
+    auto acc2 = get_acceleration(at + vel/2, vel + acc1/2);
+    auto acc3 = get_acceleration(at + vel/2 + acc1/4, vel + acc2/2);
+    auto acc4 = get_acceleration(at + vel + acc2/2, vel + acc3);
     
-    auto at2 = at + velocity / 2;
-    auto velocity2 = velocity + acc / 2;
-    
-    auto acc2 = christoffel(at2, velocity2, velocity2);
-    
-    at = at + velocity + acc2 / 2;
-    
-    velocity = velocity + acc;
+    at += vel + (acc1+acc2+acc3)/6;
+    vel += (acc1+2*acc2+2*acc3+acc4)/6;
     }
   
   EX hyperpoint numerical_exp(hyperpoint v, int steps) {