From f693761d43f9be3d53ab1f60f7488e15152237de Mon Sep 17 00:00:00 2001 From: "W. Trevor King" Date: Mon, 11 Feb 2013 10:42:51 -0500 Subject: [PATCH] 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). --- one_gaussian_bump.c | 114 ++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 110 insertions(+), 4 deletions(-) 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 . */ -#include +#include /* for malloc(), free */ +#include /* for printf() */ +#include /* for assert() */ +#include /* 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; } -- 2.26.2