one_gaussian_bump.c: Print total energy for error checks

diff --git a/one_gaussian_bump.c b/one_gaussian_bump.c
index 4adb80afc60f1776672686f668120509766ddf9a..792ca6452da51b833ca8982ac802232722b74766 100644 (file)
--- 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;
 }