one_gaussian_bump.c

   1 /*
   2 Particle collision with a Gaussian bump.
   3
   4 Copyright (C) 2013  W. Trevor King
   5
   6 This program is free software: you can redistribute it and/or modify
   7 it under the terms of the GNU General Public License as published by
   8 the Free Software Foundation, either version 3 of the License, or
   9 (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful,
  12 but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14 GNU General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program.  If not, see <http://www.gnu.org/licenses/>.
  18 */
  19
  20 #include <stdlib.h>       /* for exit(), EXIT_*, malloc(), free() */
  21 #include <stdio.h>        /* for stderr, *printf() */
  22 #include <getopt.h>       /* for getopt_long() */
  23 #include <assert.h>       /* for assert() */
  24 #include <math.h>         /* for exp(), sqrt() */
  25 #include "config.h"       /* enabling safe_strtoc() */
  26 #include "safe-strto.h"   /* for safe_strtol() and safe_strtoc() */
  27
  28 static struct option long_options[] = {
  29         {"max-time",       required_argument, 0, 't' },
  30         {"time-step",      required_argument, 0, 'd' },
  31         {"x0",             required_argument, 0, 'x' },
  32         {"y0",             required_argument, 0, 'y' },
  33         {"vx0",            required_argument, 0, 'X' },
  34         {"vx0-energy",     required_argument, 0, 'E' },
  35         {"vy0",            required_argument, 0, 'Y' },
  36         {"mass",           required_argument, 0, 'm' },
  37         {"bump-magnitude", required_argument, 0, 'k' },
  38         {"bump-width",     required_argument, 0, 'b' },
  39         {}
  40 };
  41
  42 static char *short_options = "t:d:x:y:X:E:Y:m:k:b:";
  43
  44 /* dynamic system parameters */
  45 typedef struct state_struct {
  46         double t;   /* time */
  47         size_t n;   /* number of phase space parameters */
  48         double *x;  /* position in phase space [x, y, v_x, v_y] */
  49 } state_t;
  50
  51 /* function type for phase space derivatives dx_i/dt */
  52 typedef void (*dx_dt_fn_t)(state_t *state, void *system, double *dx_dt);
  53
  54 /* static potential parameters */
  55 typedef struct gaussian_state_struct {
  56         double k;  /* magnitude */
  57         double b;  /* width */
  58 } gaussian_state_t;
  59
  60 /* static system parameters */
  61 typedef struct system_struct {
  62         double m;             /* particle mass */
  63         dx_dt_fn_t dx_dt_fn;  /* &singe_gaussian_dx_dt, etc. */
  64         void *dx_dt_state;    /* gaussian_state_t, etc. */
  65 } system_t;
  66
  67 /* Potential:    V(r) = 1/2 kb^2 exp(-r^2/b^2)
  68  * Force:        F(r) = -Divergence(V)
  69  *                    = -1/2 kb^2 * -2r/b^2 * exp(-r^2/b^2)
  70  *                    = kr exp(-r^2/b^2)
  71  *                    = k exp(-r^2/b^2) * (x ihat + y jhat)
  72  * Acceleration: a(r) = F(r)/m
  73  *                    = k/m * exp(-r^2/b^2) * (x ihat + y jhat)
  74  */
  75 void gaussian_acceleration(gaussian_state_t *state, double m,
  76                                                                                                          double x, double y,
  77                                                                                                          double *dx_dt, double *dy_dt)
  78 {
  79         double k = state->k;
  80         double b = state->b;
  81         double mag = k/m * exp(-(x*x+y*y)/b*b);
  82         *dx_dt = mag*x;
  83         *dy_dt = mag*y;
  84         return;
  85 }
  86
  87 void single_gaussian_dx_dt(state_t *state, void *system, double *dx_dt)
  88 {
  89         system_t *s = (system_t *)system;
  90
  91         assert(state->n == 4);
  92         dx_dt[0] = state->x[2];  /* dx/dt = v_x */
  93         dx_dt[1] = state->x[3];  /* dy/dt = v_y */
  94   /* dv_x/dt = a_x
  95    * dv_y/dt = a_y
  96    */
  97         gaussian_acceleration((gaussian_state_t *)s->dx_dt_state, s->m,
  98                               state->x[0], state->x[1],
  99                               &dx_dt[2], &dx_dt[3]);
 100 }
 101
 102 void step(state_t *state, system_t *system, double dt)
 103 {
 104         double *dx_dt = NULL;
 105         size_t i;
 106
 107         dx_dt = malloc(sizeof(double) * state->n);
 108         assert(dx_dt);  /* check for successful malloc() */
 109         (*system->dx_dt_fn)(state, (void *)system, dx_dt);
 110         for (i = 0; i < state->n; i++)
 111                 state->x[i] += dx_dt[i] * dt;
 112         state->t += dt;
 113         free(dx_dt);
 114 }
 115
 116 /* prints "t\tx[0]\tx[1]\t...\tx[n-1]\n" */
 117 void print_state(state_t *state)
 118 {
 119         size_t i;
 120
 121         printf("%g", state->t);
 122         for (i = 0; i < state->n; i++)
 123                 printf("\t%g", state->x[i]);
 124         printf("\n");
 125 }
 126
 127 void parse_args(int argc, char **argv, double *dt, double *t_max,
 128                 state_t *state, system_t *system)
 129 {
 130         int ch;
 131         gaussian_state_t *gstate = (gaussian_state_t *)system->dx_dt_state;
 132
 133         while (1) {
 134                 ch = getopt_long(argc, argv, short_options, long_options, NULL);
 135                 if (ch == -1)
 136                         break;
 137                 switch(ch) {
 138                 case 't':
 139                         *t_max = safe_strtod(optarg, "max-time");
 140                         break;
 141                 case 'd':
 142                         *dt = safe_strtod(optarg, "time-step");
 143                         break;
 144                 case 'x':
 145                         state->x[0] = safe_strtod(optarg, "x0");
 146                         break;
 147                 case 'y':
 148                         state->x[1] = safe_strtod(optarg, "y0");
 149                         break;
 150                 case 'X':
 151                         state->x[2] = safe_strtod(optarg, "vx0");
 152                         break;
 153                 case 'E':
 154                         state->x[2] = -sqrt(2*safe_strtod(optarg, "vx0-energy")/system->m);
 155                         break;
 156                 case 'Y':
 157                         state->x[3] = safe_strtod(optarg, "vy0");
 158                         break;
 159                 case 'm':
 160                         system->m = safe_strtod(optarg, "mass");
 161                         break;
 162                 case 'k':
 163                         gstate->k = safe_strtod(optarg, "bump-magnitude");
 164                         break;
 165                 case 'b':
 166                         gstate->b = safe_strtod(optarg, "bump-width");
 167                         break;
 168                 case '?':
 169                         exit(EXIT_FAILURE);
 170                 default:
 171                         printf("?? getopt returned character code 0%o ??\n", ch);
 172                 }
 173         }
 174         return;
 175 }
 176
 177 int main(int argc, char **argv) {
 178         double t_max = 7;
 179         double dt = 0.001;
 180         state_t state = {};
 181         double x[4] = {3.5, 0, 0, 0};
 182         system_t system = {};
 183         gaussian_state_t gaussian = {};
 184
 185         state.x = (double *)x;
 186         state.n = sizeof(x) / sizeof(double);
 187         system.m = 1;
 188         system.dx_dt_fn = &single_gaussian_dx_dt;
 189         system.dx_dt_state = (void *)&gaussian;
 190         gaussian.k = -5;
 191         gaussian.b = 1;
 192         state.x[2] = -sqrt(2*0.5/system.m);  /* v = -sqrt(2E/m),  with E = 0.5 */
 193
 194         parse_args(argc, argv, &dt, &t_max, &state, &system);
 195
 196         print_state(&state);
 197         while (state.t < t_max) {
 198                 step(&state, &system, dt);
 199                 print_state(&state);
 200         }
 201         return 0;
 202 }