Merge pull request #299 from berquist/verlet-cpp

Update C++ implementations for Verlet algorithms

Author: Gathros

contents/verlet_integration/code/c++/verlet.cpp

@@ -1,6 +1,11 @@
+#include <iomanip>
 #include <iostream>
 
-// Simple function for velocity-verlet
+struct timestep {
+  double time;
+  double vel;
+};
+
 double verlet(double pos, double acc, double dt) {
 
   double prev_pos = pos;
@@ -16,8 +21,7 @@ double verlet(double pos, double acc, double dt) {
   return time;
 }
 
-// Simple function for stormer-verlet
-double stormer_verlet(double pos, double acc, double dt) {
+timestep stormer_verlet(double pos, double acc, double dt) {
 
   double prev_pos = pos;
   double time = 0;
@@ -28,14 +32,15 @@ double stormer_verlet(double pos, double acc, double dt) {
     prev_pos = pos;
     pos = next_pos;
 
-    // The acceleration is constant, so the velocity is straightforward
+    // The acceleration is constant, so the velocity is
+    // straightforward
     vel += acc * dt;
   }
 
-  return time;
+  return timestep { time, vel };
 }
 
-double velocity_verlet(double pos, double acc, double dt) {
+timestep velocity_verlet(double pos, double acc, double dt) {
 
   double time = 0;
   double vel = 0;
@@ -45,18 +50,35 @@ double velocity_verlet(double pos, double acc, double dt) {
     vel += acc * dt;
   }
 
-  return time;
+  return timestep { time, vel };
 }
 
 int main() {
+  std::cout << std::fixed << std::setprecision(8);
+
+  // Note that depending on the simulation, you might want to have the
+  // Verlet loop outside.
+
+  // For example, if your acceleration chages as a function of time,
+  // you might need to also change the acceleration to be read into
+  // each of these functions.
+
+  double time = verlet(5.0, -10, 0.01);
+  std::cout << "Time for Verlet integration is: " \
+            << time << std::endl;
+
+  timestep timestep_sv = stormer_verlet(5.0, -10, 0.01);
+  std::cout << "Time for Stormer Verlet integration is: " \
+            << timestep_sv.time << std::endl;
+  std::cout << "Velocity for Stormer Verlet integration is: " \
+            << timestep_sv.vel << std::endl;
+
+  timestep timestep_vv = velocity_verlet(5.0, -10, 0.01);
+  std::cout << "Time for velocity Verlet integration is: " \
+            << timestep_vv.time << std::endl;
+  std::cout << "Velocity for velocity Verlet integration is: " \
+            << timestep_vv.vel << std::endl;
 
-  // Note that depending on the simulation, you might want to have the verlet
-  // loop outside.
+  return 0;
 
-  // For example, if your acceleration chages as a function of time, you might
-  // need to also change the acceleration to be read into each of these
-  // functions
-  std::cout << verlet(5.0, -10, 0.01) << std::endl;
-  std::cout << stormer_verlet(5.0, -10, 0.01) << std::endl;
-  std::cout << velocity_verlet(5.0, -10, 0.01) << std::endl;
 }

contents/verlet_integration/verlet_integration.md

@@ -33,7 +33,7 @@ Here is what it looks like in code:
 {% sample lang="jl" %}
 [import:1-13, lang:"julia"](code/julia/verlet.jl)
 {% sample lang="cpp" %}
-[import:4-17, lang:"c_cpp"](code/c++/verlet.cpp)
+[import:9-22, lang:"c_cpp"](code/c++/verlet.cpp)
 {% sample lang="c" %}
 [import:3-14, lang:"c_cpp"](code/c/verlet.c)
 {% sample lang="java" %}
@@ -78,7 +78,7 @@ However, the error for this is $$\mathcal{O}(\Delta t)$$, which is quite poor, b
 {% sample lang="jl" %}
 [import:15-31, lang:"julia"](code/julia/verlet.jl)
 {% sample lang="cpp" %}
-[import:20-36, lang:"c_cpp"](code/c++/verlet.cpp)
+[import:24-41, lang:"c_cpp"](code/c++/verlet.cpp)
 {% sample lang="c" %}
 [import:16-31, lang:"c_cpp"](code/c/verlet.c)
 {% sample lang="java" %}
@@ -137,7 +137,7 @@ Here is the velocity Verlet method in code:
 {% sample lang="jl" %}
 [import:33-45, lang:"julia"](code/julia/verlet.jl)
 {% sample lang="cpp" %}
-[import:38-49, lang:"c_cpp"](code/c++/verlet.cpp)
+[import:43-54, lang:"c_cpp"](code/c++/verlet.cpp)
 {% sample lang="c" %}
 [import:33-43, lang:"c_cpp"](code/c/verlet.c)
 {% sample lang="java" %}