packet.c

#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>

double init_fluid(double x) {
  return 100.0 * sin(2.0 * M_PI * 0.4 * x) *
         pow(M_E, -pow(x - 15.0, 2.0) / (3.6 * 2.0)) / (3.6 * sqrt(2.0 * M_PI));
}

int main() {
  FILE *fp = NULL;

  const double delta_x = 0.00001;
  const double c = 100.0;
  const double h = 1;

  double *u = calloc(sizeof(double), 400);
  double *v = calloc(sizeof(double), 400);
  assert(v != NULL && u != NULL);

  char str[256];

  int count = 0;
  for (int i = 0; i < 400; i++) {
    // if (init_fluid((double)i*30.0/400.0)<0.0) u[i] = 0.0;
    u[i] = init_fluid((double)i * 30.0 / 400.0);
  }

  for (int ctr = 0; ctr < 500; ctr++) {
    fp = fopen("pl.dat", "w");
    if (fp == NULL) {
      perror("pl.dat");
      goto clean;
    }
    for (int a = 0; a < 400; a++) {
      fprintf(fp, "%f %f\n", (double)a * 30.0 / 400.0, u[a]);
    }
    fclose(fp);
    sprintf(str,
            "gnuplot -e \"set term png; plot [0:29] [-30:30] 'pl.dat' with "
            "lines\" > frame-%d.png",
            count++);
    puts(str);
    system(str);
    for (int ii = 0; ii < 4160; ii++) {
      for (int po = 0; po < 400; po++)
        v[po] += delta_x * c * c *
                 ((u[((po - 1) < 0) ? po : po - 1] +
                   u[((po + 1) == 400) ? po : po + 1] - 2.0 * u[po]) /
                  h);
      for (int po = 0; po < 400; po++)
        u[po] += delta_x * v[po];
    }
  }
clean:
  free(u);
  free(v);
  return 0;
}