From: W. Trevor King <wking@tremily.us>
Date: Mon, 11 Feb 2013 15:42:51 +0000 (-0500)
Subject: one_gaussian_bump.c: Initial integration implementation
X-Git-Url: http://git.tremily.us/?a=commitdiff_plain;h=f693761d43f9be3d53ab1f60f7488e15152237de;p=assignment-template.git

one_gaussian_bump.c: Initial integration implementation

Currently just uses forward Euler.  Based on oscillator.c from the
phys305/oscillator branch (as of d0ebe65f, oscillator.c: disambiguate
potential state from dynamic state, 2013-01-25).
---

diff --git a/one_gaussian_bump.c b/one_gaussian_bump.c
index c242304..92b7b7b 100644
--- a/one_gaussian_bump.c
+++ b/one_gaussian_bump.c
@@ -1,7 +1,7 @@
 /*
-A simple "hello world" example in C.
+Particle collision with a Gaussian bump.
 
-Copyright (C) 2012  W. Trevor King
+Copyright (C) 2013  W. Trevor King
 
 This program is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by
@@ -17,9 +17,115 @@ You should have received a copy of the GNU General Public License
 along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */
 
-#include <stdio.h>
+#include <stdlib.h>  /* for malloc(), free */
+#include <stdio.h>   /* for printf() */
+#include <assert.h>  /* for assert() */
+#include <math.h>    /* for exp() */
+
+
+/* dynamic system parameters */
+typedef struct state_struct {
+	double t;   /* time */
+	size_t n;   /* number of phase space parameters */
+	double *x;  /* position in phase space [x, y, v_x, v_y] */
+} state_t;
+
+/* function type for phase space derivatives dx_i/dt */
+typedef void (*dx_dt_fn_t)(state_t *state, void *system, double *dx_dt);
+
+/* static potential parameters */
+typedef struct gaussian_state_struct {
+	double k;  /* magnitude */
+	double b;  /* width */
+} gaussian_state_t;
+
+/* 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. */
+} system_t;
+
+/* 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)
+ *                    = kr exp(-r^2/b^2)
+ *                    = k exp(-r^2/b^2) * (x ihat + y jhat)
+ * Acceleration: a(r) = F(r)/m
+ *                    = k/m * exp(-r^2/b^2) * (x ihat + y jhat)
+ */
+void gaussian_acceleration(gaussian_state_t *state, double m,
+													 double x, double y,
+													 double *dx_dt, double *dy_dt)
+{
+	double k = state->k;
+	double b = state->b;
+	double mag = k/m * exp(-(x*x+y*y)/b*b);
+	*dx_dt = mag*x;
+	*dy_dt = mag*y;
+	return;
+}
+
+void single_gaussian_dx_dt(state_t *state, void *system, double *dx_dt)
+{
+	system_t *s = (system_t *)system;
+
+	assert(state->n == 4);
+	dx_dt[0] = state->x[2];  /* dx/dt = v_x */
+	dx_dt[1] = state->x[3];  /* dy/dt = v_y */
+  /* dv_x/dt = a_x
+   * dv_y/dt = a_y
+   */
+	gaussian_acceleration((gaussian_state_t *)s->dx_dt_state, s->m,
+	                      state->x[0], state->x[1],
+	                      &dx_dt[2], &dx_dt[3]);
+}
+
+void step(state_t *state, system_t *system, double dt)
+{
+	double *dx_dt = NULL;
+	size_t i;
+
+	dx_dt = malloc(sizeof(double) * state->n);
+	assert(dx_dt);  /* check for successful malloc() */
+	(*system->dx_dt_fn)(state, (void *)system, dx_dt);
+	for (i = 0; i < state->n; i++)
+		state->x[i] += dx_dt[i] * dt;
+	state->t += dt;
+	free(dx_dt);
+}
+
+/* prints "t\tx[0]\tx[1]\t...\tx[n-1]\n" */
+void print_state(state_t *state)
+{
+	size_t i;
+
+	printf("%g", state->t);
+	for (i = 0; i < state->n; i++)
+		printf("\t%g", state->x[i]);
+	printf("\n");
+}
 
 int main() {
-	printf("hello, world\n");
+	double t_max = 1;
+	double dt = 0.1;
+	state_t state = {};
+	double x[4] = {};
+	system_t system = {};
+	gaussian_state_t gaussian = {};
+
+	state.x = (double *)x;
+	state.n = sizeof(x) / sizeof(double);
+	system.m = 1;
+	system.dx_dt_fn = &single_gaussian_dx_dt;
+	system.dx_dt_state = (void *)&gaussian;
+	gaussian.k = 1;
+	gaussian.b = 1;
+
+	print_state(&state);
+	while (state.t < t_max) {
+		step(&state, &system, dt);
+		print_state(&state);
+	}
 	return 0;
 }