From ea237fa1f39ddc7b461446ab79d52ba4305ce6ac Mon Sep 17 00:00:00 2001
From: "W. Trevor King" <wking@tremily.us>
Date: Tue, 12 Feb 2013 17:19:45 -0500
Subject: [PATCH] one_gaussian_bump.c: Add a midpoint integrator

---
 one_gaussian_bump.c | 44 +++++++++++++++++++++++++++++++++++++++++---
 1 file changed, 41 insertions(+), 3 deletions(-)

diff --git a/one_gaussian_bump.c b/one_gaussian_bump.c
index f3a4646..a6739b7 100644
--- a/one_gaussian_bump.c
+++ b/one_gaussian_bump.c
@@ -19,6 +19,7 @@ along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
 #include <stdlib.h>       /* for exit(), EXIT_*, malloc(), free() */
 #include <stdio.h>        /* for stderr, *printf() */
+#include <string.h>       /* for memcpy() */
 #include <getopt.h>       /* for getopt_long() */
 #include <assert.h>       /* for assert() */
 #include <math.h>         /* for exp(), sqrt() */
@@ -36,10 +37,12 @@ static struct option long_options[] = {
 	{"mass",           required_argument, 0, 'm' },
 	{"bump-magnitude", required_argument, 0, 'k' },
 	{"bump-width",     required_argument, 0, 'b' },
+	{"euler",          no_argument,       0, 'e' },
+	{"midpoint",       no_argument,       0, 'M' },
 	{}
 };
 
-static char *short_options = "t:d:x:y:X:E:Y:m:k:b:";
+static char *short_options = "t:d:x:y:X:E:Y:m:k:b:eM";
 
 /* dynamic system parameters */
 typedef struct state_struct {
@@ -140,6 +143,35 @@ void euler_step(state_t *state, system_t *system, double dt)
 	free(dx_dt);
 }
 
+void midpoint_step(state_t *state, system_t *system, double dt)
+{
+	double *x = NULL, *dx_dt = NULL;
+	size_t i, x_size = sizeof(double) * state->n;
+
+	/* allocate memory */
+	x = malloc(x_size);
+	assert(x);      /* check for successful malloc() */
+	dx_dt = malloc(x_size);
+	assert(dx_dt);  /* check for successful malloc() */
+
+	memcpy((void *)x, (void *)state->x, x_size);
+	(*system->dx_dt_fn)(state, (void *)system, dx_dt);
+	/* take a half step */
+	for (i = 0; i < state->n; i++)
+		state->x[i] += dx_dt[i] * dt / 2;
+	state->t += dt/2;
+	/* evaluate slopes at the midpoint */
+	(*system->dx_dt_fn)(state, (void *)system, dx_dt);
+	/* take a full step using the midpoint slopes */
+	for (i = 0; i < state->n; i++)
+		state->x[i] = x[i] + dx_dt[i] * dt;
+	state->t += dt/2;	
+
+	/* release memory */
+	free(x);
+	free(dx_dt);
+}
+
 /* prints "t\tx[0]\tx[1]\t...\tx[n-1]\n"
  * If system->energy_fn is defined, also print "\tenergy" before the
  * newline.
@@ -160,7 +192,7 @@ void print_state(state_t *state, system_t *system)
 }
 
 void parse_args(int argc, char **argv, double *dt, double *t_max,
-                state_t *state, system_t *system)
+                state_t *state, system_t *system, step_fn_t *step_fn)
 {
 	int ch;
 	gaussian_state_t *gstate = (gaussian_state_t *)system->dx_dt_state;
@@ -200,6 +232,12 @@ void parse_args(int argc, char **argv, double *dt, double *t_max,
 		case 'b':
 			gstate->b = safe_strtod(optarg, "bump-width");
 			break;
+		case 'e':
+			*step_fn = &euler_step;
+			break;
+		case 'M':
+			*step_fn = &midpoint_step;
+			break;
 		case '?':
 			exit(EXIT_FAILURE);
 		default:
@@ -228,7 +266,7 @@ int main(int argc, char **argv) {
 	gaussian.b = 1;
 	state.x[2] = -sqrt(2*0.5/system.m);  /* v = -sqrt(2E/m),  with E = 0.5 */
 
-	parse_args(argc, argv, &dt, &t_max, &state, &system);
+	parse_args(argc, argv, &dt, &t_max, &state, &system, &step_fn);
 
 	print_state(&state, &system);
 	while (state.t < t_max) {
-- 
2.26.2