From 864080689f887ef0830c7c9f0d0c3555b60a30eb Mon Sep 17 00:00:00 2001
From: "W. Trevor King" <wking@tremily.us>
Date: Tue, 12 Feb 2013 13:03:27 -0500
Subject: [PATCH] one_gaussian_bump.c: Print total energy for error checks

---
 one_gaussian_bump.c | 47 ++++++++++++++++++++++++++++++++++++++-------
 1 file changed, 40 insertions(+), 7 deletions(-)

diff --git a/one_gaussian_bump.c b/one_gaussian_bump.c
index 4adb80a..792ca64 100644
--- a/one_gaussian_bump.c
+++ b/one_gaussian_bump.c
@@ -51,6 +51,9 @@ typedef struct state_struct {
 /* function type for phase space derivatives dx_i/dt */
 typedef void (*dx_dt_fn_t)(state_t *state, void *system, double *dx_dt);
 
+/* function type for energy calculations */
+typedef double (*energy_fn_t)(state_t *state, void *system);
+
 /* static potential parameters */
 typedef struct gaussian_state_struct {
 	double k;  /* magnitude */
@@ -59,11 +62,20 @@ typedef struct gaussian_state_struct {
 
 /* static system parameters */
 typedef struct system_struct {
-	double m;             /* particle mass */
-	dx_dt_fn_t dx_dt_fn;  /* &singe_gaussian_dx_dt, etc. */
-	void *dx_dt_state;    /* gaussian_state_t, etc. */
+	double m;               /* particle mass */
+	dx_dt_fn_t dx_dt_fn;    /* &single_gaussian_dx_dt, etc. */
+	energy_fn_t energy_fn;  /* &single_gaussian_energy, etc. */
+	void *dx_dt_state;      /* gaussian_state_t, etc. */
 } system_t;
 
+/* Potential:    V(r) = 1/2 kb^2 exp(-r^2/b^2) */
+double gaussian_energy(gaussian_state_t *state, double x, double y)
+{
+	double k = state->k;
+	double b = state->b;
+	return 0.5 * k*b*b * exp(-(x*x+y*y)/b*b);
+}
+
 /* Potential:    V(r) = 1/2 kb^2 exp(-r^2/b^2)
  * Force:        F(r) = -Divergence(V)
  *                    = -1/2 kb^2 * -2r/b^2 * exp(-r^2/b^2)
@@ -99,6 +111,18 @@ void single_gaussian_dx_dt(state_t *state, void *system, double *dx_dt)
 	                      &dx_dt[2], &dx_dt[3]);
 }
 
+double single_gaussian_energy(state_t *state, void *system)
+{
+	system_t *s = (system_t *)system;
+	double kinetic, potential;
+
+	assert(state->n == 4);
+	kinetic = 0.5 * s->m * (state->x[2]*state->x[2] + state->x[3]*state->x[3]);
+	potential = gaussian_energy((gaussian_state_t *)s->dx_dt_state,
+	                            state->x[0], state->x[1]);
+	return kinetic + potential;
+}
+
 void step(state_t *state, system_t *system, double dt)
 {
 	double *dx_dt = NULL;
@@ -113,14 +137,22 @@ void step(state_t *state, system_t *system, double dt)
 	free(dx_dt);
 }
 
-/* prints "t\tx[0]\tx[1]\t...\tx[n-1]\n" */
-void print_state(state_t *state)
+/* prints "t\tx[0]\tx[1]\t...\tx[n-1]\n"
+ * If system->energy_fn is defined, also print "\tenergy" before the
+ * newline.
+ */
+void print_state(state_t *state, system_t *system)
 {
 	size_t i;
+	double energy;
 
 	printf("%g", state->t);
 	for (i = 0; i < state->n; i++)
 		printf("\t%g", state->x[i]);
+	if (system->energy_fn) {
+		energy = (*system->energy_fn)(state, (void *)system);
+		printf("\t%g", energy);
+	}
 	printf("\n");
 }
 
@@ -186,6 +218,7 @@ int main(int argc, char **argv) {
 	state.n = sizeof(x) / sizeof(double);
 	system.m = 1;
 	system.dx_dt_fn = &single_gaussian_dx_dt;
+	system.energy_fn = &single_gaussian_energy;
 	system.dx_dt_state = (void *)&gaussian;
 	gaussian.k = -5;
 	gaussian.b = 1;
@@ -193,10 +226,10 @@ int main(int argc, char **argv) {
 
 	parse_args(argc, argv, &dt, &t_max, &state, &system);
 
-	print_state(&state);
+	print_state(&state, &system);
 	while (state.t < t_max) {
 		step(&state, &system, dt);
-		print_state(&state);
+		print_state(&state, &system);
 	}
 	return 0;
 }
-- 
2.26.2