daviddoji/YouTube-simulations
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Add files via upload

Browse Source
This commit is contained in:
nilsberglund-orleans 2021-09-01 17:03:02 +02:00 committed by GitHub
parent 7daa596903
commit 7e2166b1d9
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
5 changed files with 402 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
global_pdes.c
View File

@ -29,6 +29,8 @@
#define D_MENGER 15 /* Menger-Sierpinski carpet */
#define D_JULIA_INT 16 /* interior of Julia set */
#define D_MENGER_ROTATED 17 /* rotated Menger-Sierpinski carpet */
#define D_PARABOLA 18 /* parabolic domain */
#define D_TWO_PARABOLAS 19 /* two facing parabolic antennas */
#define D_CIRCLES 20 /* several circles */
@ -42,6 +44,8 @@
#define C_CLOAK 5 /* invisibility cloak */
#define C_CLOAK_A 6 /* first optimized invisibility cloak */
#define C_POISSON_DISC 8 /* Poisson disc sampling */
#define C_GOLDEN_MEAN 10 /* pattern based on vertical shifts by golden mean */
#define C_GOLDEN_SPIRAL 11 /* spiral pattern based on golden mean */
#define C_SQUARE_HEX 12 /* alternating between square and hexagonal/triangular */

189
sub_wave.c
View File

@ -122,7 +122,7 @@ void color_scheme(int scheme, double value, double scale, int time, double rgb[3
/* saturation = r, luminosity = y */
switch (scheme) {
case C_LUM:
case (C_LUM):
{
hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS;
if (hue < 0.0) hue += 360.0;
@ -135,7 +135,7 @@ void color_scheme(int scheme, double value, double scale, int time, double rgb[3
hsl_to_rgb(hue, r, y, rgb);
break;
}
case C_HUE:
case (C_HUE):
{
r = 0.9;
amplitude = color_amplitude(value, scale, time);
@ -315,6 +315,33 @@ void erase_area(double x, double y, double dx, double dy)
glEnd();
}
void erase_area_rgb(double x, double y, double dx, double dy, double rgb[3])
{
double pos[2];
glColor3f(rgb[0], rgb[1], rgb[2]);
glBegin(GL_QUADS);
xy_to_pos(x - dx, y - dy, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(x + dx, y - dy, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(x + dx, y + dy, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(x - dx, y + dy, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
}
void erase_area_hsl(double x, double y, double dx, double dy, double h, double s, double l)
{
double pos[2], rgb[3];
hsl_to_rgb(h, s, l, rgb);
erase_area_rgb(x, y, dx, dy, rgb);
}
void draw_rectangle(double x1, double y1, double x2, double y2)
{
double pos[2];
@ -398,8 +425,9 @@ void init_circle_config()
/* initialise the arrays circlex, circley, circlerad and circleactive */
/* for billiard shape D_CIRCLES */
{
int i, j, n, ncirc0;
double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma;
int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted;
double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = 3.25*MU;
short int active_poisson[NMAXCIRCLES], far;
switch (CIRCLE_PATTERN) {
case (C_SQUARE):
@ -517,6 +545,65 @@ void init_circle_config()
}
break;
}
case (C_POISSON_DISC):
{
printf("Generating Poisson disc sample\n");
/* generate first circle */
circlex[0] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0);
circley[0] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN;
active_poisson[0] = 1;
n_p_active = 1;
ncircles = 1;
while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES))
{
/* randomly select an active circle */
i = rand()%(ncircles);
while (!active_poisson[i]) i = rand()%(ncircles);
// printf("Starting from circle %i at (%.3f,%.3f)\n", i, circlex[i], circley[i]);
/* generate new candidates */
naccepted = 0;
for (j=0; j<ncandidates; j++)
{
r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0);
phi = DPI*(double)rand()/RAND_MAX;
x = circlex[i] + r*cos(phi);
y = circley[i] + r*sin(phi);
// printf("Testing new circle at (%.3f,%.3f)\t", x, y);
far = 1;
for (k=0; k<ncircles; k++) if ((k!=i))
{
/* new circle is far away from circle k */
far = far*((x - circlex[k])*(x - circlex[k]) + (y - circley[k])*(y - circley[k]) >= dpoisson*dpoisson);
/* new circle is in domain */
far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > YMIN);
}
if (far) /* accept new circle */
{
printf("New circle at (%.3f,%.3f) accepted\n", x, y);
circlex[ncircles] = x;
circley[ncircles] = y;
circlerad[ncircles] = MU;
circleactive[ncircles] = 1;
active_poisson[ncircles] = 1;
ncircles++;
n_p_active++;
naccepted++;
}
// else printf("Rejected\n");
}
if (naccepted == 0) /* inactivate circle i */
{
// printf("No candidates work, inactivate circle %i\n", i);
active_poisson[i] = 0;
n_p_active--;
}
printf("%i active circles\n", n_p_active);
}
printf("Generated %i circles\n", ncircles);
break;
}
case (C_GOLDEN_MEAN):
{
ncircles = 300;
@ -590,7 +677,7 @@ void init_circle_config()
if ((circley[i] < YMAX + MU)&&(circley[i] > YMIN - MU)) circleactive[i] = 1;
// printf("i = %i, circlex = %.3lg, circley = %.3lg\n", i, circlex[i], circley[i]);
}
break;
break;
}
case (C_SQUARE_HEX):
{
@ -648,11 +735,11 @@ void init_circle_config()
}
int xy_in_billiard(x, y)
int xy_in_billiard(double x, double y)
/* returns 1 if (x,y) represents a point in the billiard */
double x, y;
// double x, y;
{
double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy;
double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width;
int i, j, k, k1, k2, condition;
switch (B_DOMAIN) {
@ -774,6 +861,21 @@ double x, y;
}
return(1);
}
case (D_PARABOLA):
{
return(x > 0.25*y*y/LAMBDA - LAMBDA);
}
case (D_TWO_PARABOLAS):
{
x1 = 0.25*y*y/MU - MU - LAMBDA;
x2 = -x1;
width = 0.25*MU;
if (width > 0.2) width = 0.2;
if (vabs(y) > 1.5*MU) return(1);
else if ((x < x1 - width)||(x > x2 + width)) return(1);
else if ((x > x1)&&(x < x2)) return(1);
else return(0);
}
case (D_CIRCLES):
{
for (i = 0; i < ncircles; i++)
@ -944,7 +1046,7 @@ int ij_in_billiard(int i, int j)
void draw_billiard() /* draws the billiard boundary */
{
double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l;
double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width;
int i, j, k, k1, k2, mr2;
if (BLACK) glColor3f(1.0, 1.0, 1.0);
@ -1205,6 +1307,75 @@ void draw_billiard() /* draws the billiard boundary */
}
break;
}
case (D_PARABOLA):
{
dy = (YMAX - YMIN)/(double)NSEG;
glBegin(GL_LINE_STRIP);
for (i = 0; i < NSEG+1; i++)
{
y = YMIN + dy*(double)i;
x = 0.25*y*y/LAMBDA - LAMBDA;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
draw_circle(0.0, 0.0, r, NSEG);
}
break;
}
case (D_TWO_PARABOLAS):
{
dy = 3.0*MU/(double)NSEG;
width = 0.25*MU;
if (width > 0.2) width = 0.2;
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -1.5*MU + dy*(double)i;
x = 0.25*y*y/MU - MU - LAMBDA;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i = 0; i < NSEG+1; i++)
{
y = 1.5*MU - dy*(double)i;
x = 0.25*y*y/MU - (MU + width) - LAMBDA;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
glBegin(GL_LINE_LOOP);
for (i = 0; i < NSEG+1; i++)
{
y = -1.5*MU + dy*(double)i;
x = LAMBDA + MU - 0.25*y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i = 0; i < NSEG+1; i++)
{
y = 1.5*MU - dy*(double)i;
x = LAMBDA + (MU + width) - 0.25*y*y/MU;
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
draw_circle(-LAMBDA, 0.0, r, NSEG);
draw_circle(LAMBDA, 0.0, r, NSEG);
}
break;
}
case (D_CIRCLES):
{
glLineWidth(BOUNDARY_WIDTH);

174
sub_wave_comp.c
View File

@ -7,9 +7,10 @@ void init_circle_config_half(int pattern, int top)
/* initialise the arrays circlex, circley, circlerad and circleactive */
/* for billiard shape D_CIRCLES */
{
int i, j, n, ncirc0;
double dx, dy, p, phi, r, ra[5], sa[5], height, y = 0.0, gamma, ymean, ytop, ybottom;
int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted, nnew;
double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, ymean, ytop, ybottom, dpoisson = 3.05*MU;
short int active_poisson[NMAXCIRCLES], far;
ymean = 0.5*(YMIN + YMAX);
switch (pattern) {
case (C_SQUARE):
@ -164,6 +165,80 @@ void init_circle_config_half(int pattern, int top)
ncircles += 1;
break;
}
case (C_POISSON_DISC):
{
printf("Generating Poisson disc sample\n");
/* generate first circle */
n = ncircles;
circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0);
if (top) y = ymean + (YMAX-ymean)*(double)rand()/RAND_MAX;
else y = ymean + (YMIN-ymean)*(double)rand()/RAND_MAX;
circley[n] = y;
circlerad[n] = MU;
circleactive[n] = 1;
circletop[n] = top;
active_poisson[n] = 1;
n_p_active = 1;
ncirc0 = 1;
while ((n_p_active > 0)&&(ncircles + ncirc0 < NMAXCIRCLES))
{
/* randomly select an active circle */
i = rand()%(ncirc0);
n = ncircles + i;
while (!active_poisson[ncircles + i]) i = rand()%(ncirc0);
// printf("Starting from circle %i at (%.3f,%.3f)\n", i, circlex[i], circley[i]);
/* generate new candidates */
naccepted = 0;
for (j=0; j<ncandidates; j++)
{
r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0);
phi = DPI*(double)rand()/RAND_MAX;
x = circlex[n] + r*cos(phi);
y = circley[n] + r*sin(phi);
// printf("Testing new circle at (%.3f,%.3f)\t", x, y);
far = 1;
for (k=ncircles; k<ncircles + ncirc0; k++) if ((k!=n))
{
/* new circle is far away from circle k */
far = far*((x - circlex[k])*(x - circlex[k]) + (y - circley[k])*(y - circley[k]) >= dpoisson*dpoisson);
/* new circle is in domain */
if (top) far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > 0.0);
else far = far*(vabs(x) < LAMBDA)*(y > YMIN)*(y < 0.0);
}
if (far) /* accept new circle */
{
printf("New circle at (%.3f,%.3f) accepted\n", x, y);
nnew = ncircles + ncirc0;
circlex[nnew] = x;
circley[nnew] = y;
circlerad[nnew] = MU;
circleactive[nnew] = 1;
active_poisson[nnew] = 1;
circleactive[nnew] = 1;
circletop[nnew] = top;
ncirc0++;
n_p_active++;
naccepted++;
}
// else printf("Rejected\n");
}
if (naccepted == 0) /* inactivate circle i */
{
// printf("No candidates work, inactivate circle %i\n", ncircles + i);
active_poisson[ncircles + i] = 0;
n_p_active--;
}
printf("%i active circles\n", n_p_active);
// sleep(1);
}
printf("Already existing: %i circles\n", ncircles);
ncircles += ncirc0;
printf("Generated %i circles\n", ncirc0);
printf("Total: %i circles\n", ncircles);
break;
}
case (C_GOLDEN_MEAN):
{
gamma = (sqrt(5.0) - 1.0)*0.5; /* golden mean */
@ -224,6 +299,43 @@ void init_circle_config_half(int pattern, int top)
break;
}
case (C_GOLDEN_SPIRAL):
{
circlex[ncircles] = 0.0;
circley[ncircles] = 0.0;
circlerad[ncircles] = MU;
circleactive[ncircles] = 1;
circletop[ncircles] = top;
gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */
phi = 0.0;
r0 = 2.0*MU;
r = r0 + MU;
for (i=0; i<1000; i++)
{
x = r*cos(phi);
y = r*sin(phi);
phi += gamma;
r += MU*r0/r;
if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU))
{
circlex[ncircles] = x;
circley[ncircles] = y;
circlerad[ncircles] = MU;
if (((top)&&(circley[ncircles] < YMAX + MU)&&(circley[ncircles] > ymean - MU))
||((!top)&&(circley[ncircles] < ymean + MU)&&(circley[ncircles] > YMIN - MU)))
{
circleactive[ncircles] = 1;
circletop[ncircles] = top;
ncircles++;
}
}
}
break;
}
case (C_ONE):
{
circlex[ncircles] = 0.0;
@ -276,6 +388,13 @@ void init_circle_config_comp()
init_circle_config_half(CIRCLE_PATTERN_B, 0);
}
void init_circle_config_energy()
/* initialise the arrays circlex, circley, circlerad and circleactive */
/* for billiard shape D_CIRCLES */
{
ncircles = 0;
init_circle_config_half(CIRCLE_PATTERN, 0);
}
int xy_in_billiard_half(double x, double y, int domain, int pattern, int top)
/* returns 1 if (x,y) represents a point in the billiard */
@ -551,9 +670,12 @@ void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX],
/* compute total energy in top/bottom left/right boxes */
{
int i, j, ij[2];
double energy = 0.0, pos, xleft = XMAX, xright = XMIN;
double energy = 0.0, rescale, pos, xleft = XMAX, xright = XMIN, emax = 1.2;
double energy_ij[NX][NY];
static short int first = 1;
static int ileft, iright, jmid = NY/2;
static double sqremax;
if (first) /* compute box limits */
{
@ -572,48 +694,74 @@ void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX],
printf("xleft = %.3lg, xright = %.3lg", xleft, xright);
}
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
energy_ij[i][j] = compute_energy(phi, psi, xy_in, i, j);
/* prevent local energy from growing too large */
if (FLOOR)
{
for (i=10; i<NX; i++)
for (j=0; j<NY; j++)
if ((xy_in[i][j])&&(energy_ij[i][j] > emax))
{
rescale = sqrt(emax/energy_ij[i][j]);
if (j%100 == 0) printf("Rescaling at (%i,%i) by %.5lg\n", i, j, rescale);
phi[i][j] = phi[i][j]*rescale;
psi[i][j] = psi[i][j]*rescale;
}
else if (energy_ij[i][j] > 0.1*emax)
{
rescale = sqrt(0.1*emax/energy_ij[i][j]);
if (j%10 == 0) printf("Rescaling at (%i,%i) by %.5lg\n", i, j, rescale);
phi[i][j] = phi[i][j]*rescale;
psi[i][j] = psi[i][j]*rescale;
}
}
/* top left box */
for (i=0; i<ileft; i++)
for (j=jmid; j<NY; j++)
energy += compute_energy(phi, psi, xy_in, i, j);
energy += energy_ij[i][j];
energies[0] = energy;
/* top middle box */
energy = 0.0;
for (i=ileft; i<iright; i++)
for (j=jmid; j<NY; j++)
energy += compute_energy(phi, psi, xy_in, i, j);
energy += energy_ij[i][j];
energies[1] = energy;
/* top right box */
energy = 0.0;
for (i=iright; i<NX; i++)
for (j=jmid; j<NY; j++)
energy += compute_energy(phi, psi, xy_in, i, j);
energy += energy_ij[i][j];
energies[2] = energy;
/* bottom left box */
energy = 0.0;
for (i=0; i<ileft; i++)
for (j=0; j<jmid; j++)
energy += compute_energy(phi, psi, xy_in, i, j);
energy += energy_ij[i][j];
energies[3] = energy;
/* bottom middle box */
energy = 0.0;
for (i=ileft; i<iright; i++)
for (j=0; j<jmid; j++)
energy += compute_energy(phi, psi, xy_in, i, j);
energy += energy_ij[i][j];
energies[4] = energy;
/* bottom right box */
energy = 0.0;
for (i=iright; i<NX; i++)
for (j=0; j<jmid; j++)
energy += compute_energy(phi, psi, xy_in, i, j);
energy += energy_ij[i][j];
energies[5] = energy;
printf("Energies: %.5lg, %.5lg, %.5lg\n %.5lg, %.5lg, %.5lg\n", energies[0], energies[1], energies[2], energies[3], energies[4], energies[5]);
// printf("Energies: %.5lg, %.5lg, %.5lg\n %.5lg, %.5lg, %.5lg\n", energies[0], energies[1], energies[2], energies[3], energies[4], energies[5]);
}
void print_energies(double energies[6], double top_energy, double bottom_energy)
@ -622,8 +770,8 @@ void print_energies(double energies[6], double top_energy, double bottom_energy)
double ytop, ybot, pos[2], centerx = -0.075;
ytop = YMAX - 0.1;
ybot = -0.1;
// ybot = YMIN + 0.05;
// ybot = -0.1;
ybot = YMIN + 0.05;
erase_area(XMIN + 0.175, ytop + 0.025, 0.1, 0.05);

52
wave_billiard.c
View File

@ -42,8 +42,8 @@
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <omp.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */
#define MOVIE 1 /* set to 1 to generate movie */
#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */
/* General geometrical parameters */
@ -52,11 +52,6 @@
#define NX 1280 /* number of grid points on x axis */
#define NY 720 /* number of grid points on y axis */
// #define NX 640 /* number of grid points on x axis */
// #define NY 360 /* number of grid points on y axis */
/* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */
/* but will multiply run time by 4 */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
@ -67,15 +62,15 @@
/* Choice of the billiard table */
#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 19 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 11 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_pdes.c */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
#define LAMBDA 0.85 /* parameter controlling the dimensions of domain */
#define MU 0.03 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.0 /* parameter controlling the dimensions of domain */
#define MU 1.25 /* parameter controlling the dimensions of domain */
#define NPOLY 3 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 4 /* depth of computation of Menger gasket */
@ -83,7 +78,7 @@
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 15 /* number of grid point for grid of disks */
#define NGRIDX 16 /* number of grid point for grid of disks */
#define NGRIDY 20 /* number of grid point for grid of disks */
/* You can add more billiard tables by adapting the functions */
@ -102,7 +97,7 @@
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 1.0e-6 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
#define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */
#define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */
#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */
#define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */
#define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */
@ -113,26 +108,28 @@
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 3
#define B_COND 2
/* Parameters for length and speed of simulation */
#define NSTEPS 3000 /* number of frames of movie */
#define NVID 20 /* number of iterations between images displayed on screen */
#define NSTEPS 5000 /* number of frames of movie */
#define NVID 25 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
#define INITIAL_TIME 250 /* time after which to start saving frames */
#define INITIAL_TIME 0 /* time after which to start saving frames */
#define BOUNDARY_WIDTH 2 /* width of billiard boundary */
#define PAUSE 1000 /* number of frames after which to pause */
#define PAUSE 1000 /* number of frames after which to pause */
#define PSLEEP 1 /* sleep time during pause */
#define SLEEP1 1 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */
#define SLEEP2 1 /* final sleeping time */
#define MID_FRAMES 20 /* number of still frames between parts of two-part movie */
#define END_FRAMES 100 /* number of still frames at end of movie */
/* Plot type, see list in global_pdes.c */
#define PLOT 0
#define PLOT 1
#define PLOT_B 1 /* plot type for second movie */
#define PLOT_B 0 /* plot type for second movie */
/* Color schemes */
@ -141,9 +138,9 @@
#define COLOR_SCHEME 1 /* choice of color scheme, see list in global_pdes.c */
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 1.5 /* sensitivity of color on wave amplitude */
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 2000.0 /* scaling factor for energy representation */
#define E_SCALE 500.0 /* scaling factor for energy representation */
// #define E_SCALE 2500.0 /* scaling factor for energy representation */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
@ -184,7 +181,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
// c = COURANT;
// cc = courant2;
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y,c,cc,gamma)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j])
@ -332,7 +329,8 @@ void animation()
/* initialize wave with a drop at one point, zero elsewhere */
// init_wave_flat(phi, psi, xy_in);
init_planar_wave(XMIN + 0.015, 0.0, phi, psi, xy_in);
init_wave(-LAMBDA, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.015, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.02, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 0.8, 0.0, phi, psi, xy_in);
// init_wave(-1.5, 0.0, phi, psi, xy_in);
@ -412,13 +410,13 @@ void animation()
draw_billiard();
glutSwapBuffers();
}
for (i=0; i<20; i++) save_frame();
for (i=0; i<MID_FRAMES; i++) save_frame();
if (DOUBLE_MOVIE)
{
draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
draw_billiard();
glutSwapBuffers();
for (i=0; i<20; i++) save_frame_counter(NSTEPS + 21 + counter + i);
for (i=0; i<END_FRAMES; i++) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
}
s = system("mv wave*.tif tif_wave/");

54
wave_comparison.c
View File

@ -42,7 +42,7 @@
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <omp.h>
#define MOVIE 0 /* set to 1 to generate movie */
#define MOVIE 1 /* set to 1 to generate movie */
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
@ -59,11 +59,11 @@
/* Choice of the billiard table */
#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */
#define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */
#define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN_B 10 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN 11 /* pattern of circles, see list in global_pdes.c */
#define CIRCLE_PATTERN_B 8 /* pattern of circles, see list in global_pdes.c */
#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */
#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */
@ -78,7 +78,7 @@
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NGRIDX 15 /* number of grid point for grid of disks */
#define NGRIDX 20 /* number of grid point for grid of disks */
#define NGRIDY 20 /* number of grid point for grid of disks */
/* You can add more billiard tables by adapting the functions */
@ -86,16 +86,22 @@
/* Physical parameters of wave equation */
#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */
#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */
#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */
#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */
#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */
#define OMEGA 0.0035 /* frequency of periodic excitation */
#define AMPLITUDE 1.0 /* amplitude of periodic excitation */
#define OMEGA 0.0 /* frequency of periodic excitation */
#define AMPLITUDE 0.025 /* amplitude of periodic excitation */
#define COURANT 0.02 /* Courant number */
#define COURANTB 0.0075 /* Courant number in medium B */
#define COURANTB 0.004 /* Courant number in medium B */
// #define COURANTB 0.005 /* Courant number in medium B */
// #define COURANTB 0.008 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMAB 1.0e-5 /* damping factor in wave equation */
// #define GAMMA 1.0e-8 /* damping factor in wave equation */
#define GAMMAB 1.0e-8 /* damping factor in wave equation */
// #define GAMMAB 1.0e-6 /* damping factor in wave equation */
// #define GAMMAB 2.0e-4 /* damping factor in wave equation */
// #define GAMMAB 2.5e-4 /* damping factor in wave equation */
#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */
#define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */
#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */
@ -112,7 +118,7 @@
/* Parameters for length and speed of simulation */
#define NSTEPS 6000 /* number of frames of movie */
#define NSTEPS 3200 /* number of frames of movie */
#define NVID 25 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
#define INITIAL_TIME 200 /* time after which to start saving frames */
@ -123,10 +129,11 @@
#define PSLEEP 1 /* sleep time during pause */
#define SLEEP1 1 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */
#define END_FRAMES 100 /* number of still frames at end of movie */
/* Plot type, see list in global_pdes.c */
#define PLOT 1
#define PLOT 0
/* Color schemes */
@ -135,7 +142,7 @@
#define COLOR_SCHEME 1 /* choice of color scheme, see list in global_pdes.c */
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 1.5 /* sensitivity of color on wave amplitude */
#define SLOPE 50.0 /* sensitivity of color on wave amplitude */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 2000.0 /* scaling factor for energy representation */
@ -148,7 +155,7 @@
/* For debugging purposes only */
#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 10.0 /* max value of wave amplitude */
#define VMAX 5.0 /* max value of wave amplitude */
#include "global_pdes.c" /* constants and global variables */
@ -174,7 +181,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
time++;
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y,c,cc,gamma)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j])
@ -228,8 +235,9 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
iminus = (i-1); if (iminus == -1) iminus = 0;
if (j < jmid) /* lower half */
{
jplus = (j+1) % jmid;
jminus = (j-1) % jmid;
jplus = (j+1);
if (jplus >= jmid) jplus -= jmid;
jminus = (j-1);
if (jminus < 0) jminus += jmid;
}
else /* upper half */
@ -293,6 +301,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
/* add oscillating boundary condition on the left */
if ((i == 0)&&(OSCILLATE_LEFT)) phi_out[i][j] = AMPLITUDE*cos((double)time*OMEGA);
psi_out[i][j] = x;
if (FLOOR)
@ -344,7 +353,8 @@ void animation()
courantb2 = COURANTB*COURANTB;
/* initialize wave with a drop at one point, zero elsewhere */
int_planar_wave_comp(XMIN + 0.015, 0.0, phi, psi, xy_in);
init_wave_flat_comp(phi, psi, xy_in);
// int_planar_wave_comp(XMIN + 0.015, 0.0, phi, psi, xy_in);
// int_planar_wave_comp(XMIN + 0.5, 0.0, phi, psi, xy_in);
printf("initializing wave\n");
// int_planar_wave_comp(XMIN + 0.1, 0.0, phi, psi, xy_in);
@ -434,7 +444,7 @@ void animation()
if (MOVIE)
{
for (i=0; i<20; i++) save_frame();
for (i=0; i<END_FRAMES; i++) save_frame();
s = system("mv wave*.tif tif_wave/");
}
for (i=0; i<NX; i++)