Skip to content

Instantly share code, notes, and snippets.

@harieamjari

harieamjari/diwav.c

Last active Jan 26, 2021
Embed
What would you like to do?
Wave equation
/*
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;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment