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 c2423049f8ba3e94dc3161704eb9fce828895f8c..92b7b7b0dcf8d17bf322d2c1609e7fba10fed9b9 100644 (file)
--- 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;
 }