From 60e51757389abd6154f3209b230583d295eb39e2 Mon Sep 17 00:00:00 2001
From: Lukas Winkler <git@lw1.at>
Date: Tue, 13 Apr 2021 17:10:58 +0200
Subject: [PATCH] log energy every 100 years

---
 heartbeat/heartbeat.c | 10 ++++++++++
 water_sim.py          |  7 +++++--
 2 files changed, 15 insertions(+), 2 deletions(-)

diff --git a/heartbeat/heartbeat.c b/heartbeat/heartbeat.c
index 67ab81c..1300c84 100644
--- a/heartbeat/heartbeat.c
+++ b/heartbeat/heartbeat.c
@@ -18,6 +18,13 @@ int hb_sun_collision_index = 0;
 
 int needs_synchronize = 0;
 
+FILE *logfile;
+
+void init_logfile(char *filename) {
+    FILE *f = fopen(filename, "a");
+    logfile = f;
+}
+
 void heartbeat(struct reb_simulation *sim) {
     if ((sim->steps_done % 10) == 0) {
         const struct reb_particle *const particles = sim->particles;
@@ -60,4 +67,7 @@ void heartbeat(struct reb_simulation *sim) {
             }
         }
     }
+    if ((sim->steps_done % 10000) == 0) { // ~ every 100 years
+        fprintf(logfile, "%f, %f\n", sim->t, reb_tools_energy(sim));
+    }
 }
diff --git a/water_sim.py b/water_sim.py
index 49dd8b9..9af5ddb 100644
--- a/water_sim.py
+++ b/water_sim.py
@@ -1,8 +1,7 @@
 import re
 import time
-from ctypes import Structure, c_uint32, c_double, c_uint, cdll, c_int
+from ctypes import Structure, c_uint32, c_double, c_uint, cdll, c_int, create_string_buffer, c_char_p
 from dataclasses import dataclass
-from math import radians
 from pathlib import Path
 from shutil import copy
 from sys import argv
@@ -173,6 +172,10 @@ def main(fn: Path, testrun=False):
 
     check_heartbeat_needs_recompile()
     clibheartbeat = cdll.LoadLibrary("heartbeat/heartbeat.so")
+    clibheartbeat.init_logfile.argtypes = [c_char_p]
+    logfile = create_string_buffer(128)
+    logfile.value = str(fn.with_suffix(".energylog.csv")).encode()
+    clibheartbeat.init_logfile(logfile)
     sim.heartbeat = clibheartbeat.heartbeat
     innermost_semimajor_axis = third_kepler_law(
         orbital_period=sim.dt * year * MIN_TIMESTEP_PER_ORBIT