Wave equation

diwav.c

/*
DO WHAT THE FUCK YOU WANT TO BUT IT'S NOT MY FAULT PUBLIC LICENSE
Version 1, October 2013

Copyright (c) 2020 Al-buharie Amjari <healer.harie@gmail.com>

Everyone is permitted to copy and distribute verbatim or modified copies
of this license document, and changing it is allowed as long as the name
is changed.
DO WHAT THE FUCK YOU WANT TO BUT IT'S NOT MY FAULT PUBLIC LICENSE TERMS
AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION:

0. You just DO WHAT THE FUCK YOU WANT TO.

1. Do not hold the author(s), creator(s), developer(s) or distributor(s)
liable for anything that happens or goes wrong with your use of the work.
*/

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

double init_fluid(double x){
   return -pow(x-15.0,2.0)+25.0;
}

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;
      else u[i] = init_fluid((double)i*30.0/400.0);
   }

   for (int ctr = 0; ctr < 500*4160; ctr++){
         for (int po = 0; po < 400; po++)
            v[po] += delta_x*c*c*((u[((po-1)<0)?0:po-1] + u[((po+1)>400)?400:po+1]-2.0*u[po])/h);
         for (int po = 0; po < 400; po++) u[po] += delta_x*v[po];

      if (ctr%4160==0){
              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\" > $down/fluid/frame-%d.png", count++);
      puts(str);
      system(str);}

   }
   clean:
   free(u);free(v);
   return 0;
}