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

Add files via upload

This commit is contained in:
nilsberglund-orleans 2021-07-18 15:06:59 +02:00 committed by GitHub
parent 279a6e8801
commit f56b8a0ab4
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 951 additions and 425 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
drop_billiard.c
View File

@ -33,32 +33,28 @@
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
#define XMIN -1.8
#define XMAX 1.8 /* x interval */
#define YMIN -0.91
#define YMAX 1.115 /* y interval for 9/16 aspect ratio */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
// #define YMIN -1.125
// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
// #define XMIN -1.8
// #define XMAX 1.8 /* x interval */
// #define YMIN -0.91
// #define YMAX 1.115 /* y interval for 9/16 aspect ratio */
#define SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */
/* Choice of the billiard table */
/* Choice of the billiard table, see global_particles.c */
#define B_DOMAIN 9 /* choice of domain shape */
#define B_DOMAIN 13 /* choice of domain shape */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_STADIUM 2 /* stadium-shaped domain */
#define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular billiard */
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_REULEAUX 9 /* Reuleaux and star shapes */
#define D_FLOWER 10 /* Bunimovich flower */
#define D_ALT_REU 11 /* alternating between star and Reuleaux */
#define CIRCLE_PATTERN 0 /* pattern of circles */
#define NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */
// #define NCX 10 /* number of circles in x direction */
// #define NCY 15 /* number of circles in y direction */
#define NCX 15 /* number of circles in x direction */
#define NCY 20 /* number of circles in y direction */
// #define LAMBDA 1.0 /* parameter controlling shape of billiard */
#define LAMBDA 1.124950941 /* sin(36°)/sin(31.5°) for 5-star shape with 45° angles */
@ -72,16 +68,17 @@
#define APOLY -1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define DRAW_BILLIARD 1 /* set to 1 to draw billiard */
#define DRAW_CONSTRUCTION_LINES 1 /* set to 1 to draw additional construction lines for billiard */
#define PERIODIC_BC 0 /* set to 1 to enforce periodic boundary conditions when drawing particles */
#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */
#define NPART 100000 /* number of particles */
#define NPART 50000 /* number of particles */
#define NPARTMAX 100000 /* maximal number of particles after resampling */
#define NSTEPS 4000 /* number of frames of movie */
#define TIME 15 /* time between movie frames, for fluidity of real-time simulation */
#define NSTEPS 5500 /* number of frames of movie */
#define TIME 40 /* time between movie frames, for fluidity of real-time simulation */
#define DPHI 0.0001 /* integration step */
#define NVID 10 /* number of iterations between images displayed on screen */
#define NVID 20 /* number of iterations between images displayed on screen */
/* Decreasing TIME accelerates the animation and the movie */
/* For constant speed of movie, TIME*DPHI should be kept constant */
@ -91,16 +88,17 @@
/* simulation parameters */
#define LMAX 0.01 /* minimal segment length triggering resampling */
#define LPERIODIC 1.0 /* lines longer than this are not drawn (useful for Sinai billiard) */
#define LCUT 1000.0 /* controls the max size of segments not considered as being cut */
#define LPERIODIC 0.1 /* lines longer than this are not drawn (useful for Sinai billiard) */
#define LCUT 2.0 /* controls the max size of segments not considered as being cut */
#define DMIN 0.02 /* minimal distance to boundary for triggering resampling */
#define CYCLE 0 /* set to 1 for closed curve (start in all directions) */
#define ORDER_COLORS 1 /* set to 1 if colors should be drawn in order */
/* color and other graphical parameters */
#define NCOLORS 10 /* number of colors */
#define NCOLORS 16 /* number of colors */
#define COLORSHIFT 200 /* hue of initial color */
#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */
#define NSEG 100 /* number of segments of boundary */
#define BILLIARD_WIDTH 4 /* width of billiard */
#define FRONT_WIDTH 3 /* width of wave front */
@ -120,6 +118,7 @@
#define DPI 6.283185307
#define PID 1.570796327
#include "global_particles.c"
#include "sub_part_billiard.c"
@ -394,8 +393,11 @@ double *configs[NPARTMAX];
if (configs[i][2]<0.0)
{
vbilliard(configs[i]);
color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS;
if (!RAINBOW_COLOR)
{
color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS;
}
}
configs[i][2] += DPHI;
@ -449,7 +451,8 @@ void animation()
for (i=0; i<NPARTMAX; i++)
configs[i] = (double *)malloc(8*sizeof(double));
init_drop_config(0.0, 0.1, 0.0, DPI, configs);
init_drop_config(-1.0 + 0.3*sqrt(2.0), -1.0 + 0.5*sqrt(2.0), 0.0, DPI, configs);
// init_drop_config(-0.5, -0.5, 0.0, DPI, configs);
// init_boundary_config(1.5, 1.5, 0.0, PI, configs);
@ -463,6 +466,9 @@ void animation()
for (i=0; i<NPARTMAX; i++) color[i] = 0;
if (RAINBOW_COLOR) /* rainbow color scheme */
for (i=0; i<NPART; i++) color[i] = (i*NCOLORS)/NPART;
sleep(SLEEP1);

33
global_particles.c Normal file
View File

@ -0,0 +1,33 @@
double circlex[NMAXCIRCLES], circley[NMAXCIRCLES], circlerad[NMAXCIRCLES]; /* position and radius of circular scatters */
short int circleactive[NMAXCIRCLES]; /* tells which circular scatters are active */
int ncircles = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */
/* some basic math */
#define PI 3.141592654
#define DPI 6.283185307
#define PID 1.570796327
/* Choice of the billiard table */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_STADIUM 2 /* stadium-shaped domain */
#define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular billiard */
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_REULEAUX 9 /* Reuleaux and star shapes */
#define D_FLOWER 10 /* Bunimovich flower */
#define D_ALT_REU 11 /* alternating between star and Reuleaux */
#define D_ANGLE 12 /* angular sector */
#define D_LSHAPE 13 /* L-shaped billiard for conical singularity */
#define D_GENUSN 14 /* polygon with identifies opposite sides */
#define D_CIRCLES 20 /* several circles */
#define C_FOUR_CIRCLES 0 /* four circles almost touching each other */
#define C_SQUARE 1 /* square grid of circles */
#define C_HEX 2 /* hexagonal/triangular grid of circles */

90
global_pdes.c Normal file
View File

@ -0,0 +1,90 @@
/* Global variables and parameters for wave_billiard, heat and schrodinger */
/* Basic math */
#define PI 3.141592654
#define DPI 6.283185307
#define PID 1.570796327
/* shape of domain */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_STADIUM 2 /* stadium-shaped domain */
#define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular billiard */
#define D_FLAT 6 /* flat interface */
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_YOUNG 9 /* Young diffraction slits */
#define D_GRATING 10 /* diffraction grating */
#define D_EHRENFEST 11 /* Ehrenfest urn type geometry */
/* The following 3 types are superseded by D_CIRCLES and can be removed later */
#define D_DISK_GRID 12 /* grid of disks */
#define D_DISK_HEX 13 /* haxagonl grid of disks */
#define D_DISK_PERCOL 14 /* random grid of percolation type */
#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_CIRCLES 20 /* several circles */
#define NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */
#define C_SQUARE 0 /* square grid of circles */
#define C_HEX 1 /* hexagonal/triangular grid of circles */
#define C_RAND_DISPLACED 2 /* randomly displaced square grid */
#define C_RAND_PERCOL 3 /* random percolation arrangement */
#define C_RAND_POISSON 4 /* random Poisson point process */
/* Billiard tables for heat equation */
#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
#define D_MENGER_HEATED 22 /* Menger gasket with different temperatures */
#define D_MENGER_H_OPEN 23 /* Menger gasket with different temperatures and larger domain */
#define D_MANDELBROT 24 /* Mandelbrot set */
#define D_JULIA 25 /* Julia set */
#define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */
/* Boundary conditions */
#define BC_DIRICHLET 0 /* Dirichlet boundary conditions */
#define BC_PERIODIC 1 /* periodic boundary conditions */
#define BC_ABSORBING 2 /* absorbing boundary conditions (beta version) */
#define BC_VPER_HABS 3 /* vertically periodic and horizontally absorbing boundary conditions */
/* 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 */
/* Plot types */
/* For wave equation */
#define P_AMPLITUDE 0 /* plot amplitude of wave */
#define P_ENERGY 1 /* plot energy of wave */
#define P_MIXED 2 /* plot amplitude in upper half, energy in lower half */
/* For Schrodinger equation */
#define P_MODULE 10 /* plot module of wave function squared */
#define P_PHASE 11 /* plot phase of wave function */
#define P_REAL 12 /* plot real part */
#define P_IMAGINARY 13 /* plot imaginary part */
/* Color schemes */
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
#define C_HUE 1 /* color scheme modifies hue */
#define C_PHASE 2 /* color scheme shows phase */
double circlex[NMAXCIRCLES], circley[NMAXCIRCLES], circlerad[NMAXCIRCLES]; /* position and radius of circular scatterers */
short int circleactive[NMAXCIRCLES]; /* tells which circular scatters are active */
int ncircles = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */
double julia_x = -0.5, julia_y = 0.5; /* parameters for Julia sets */
// double julia_x = 0.33267, julia_y = 0.06395; /* parameters for Julia sets */
// double julia_x = 0.37468, julia_y = 0.21115; /* parameters for Julia sets */

209
heat.c
View File

@ -50,46 +50,25 @@
/* 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 */
// #define XMIN -1.5
// #define XMAX 2.5 /* x interval */
// #define XMIN -2.0
// #define XMAX 2.0 /* x interval */
#define XMIN -2.5
#define XMAX 1.5 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define JULIA_SCALE 1.1 /* scaling for Julia sets */
// #define JULIA_SCALE 0.8 /* scaling for Julia sets */
#define JULIA_SCALE 1.0 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 25 /* choice of domain shape */
#define B_DOMAIN 26 /* choice of domain shape, see list in global_pdes.c */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_STADIUM 2 /* stadium-shaped domain */
#define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular billiard */
#define D_FLAT 6 /* flat interface */
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_YOUNG 9 /* Young diffraction slits */
#define D_GRATING 10 /* diffraction grating */
#define D_EHRENFEST 11 /* Ehrenfest urn type geometry */
#define CIRCLE_PATTERN 0 /* pattern of circles, see list in global_pdes.c */
#define D_MENGER 15 /* Menger-Sierpinski carpet */
#define D_JULIA_INT 16 /* interior of Julia set */
#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 */
/* Billiard tables for heat equation */
#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
#define D_MENGER_HEATED 22 /* Menger gasket with different temperatures */
#define D_MENGER_H_OPEN 23 /* Menger gasket with different temperatures and larger domain */
#define D_MANDELBROT 24 /* Mandelbrot set */
#define D_JULIA 25 /* Julia set */
#define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */
#define LAMBDA 0.7 /* parameter controlling the dimensions of domain */
#define LAMBDA -1.0 /* parameter controlling the dimensions of domain */
#define MU 0.1 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
@ -98,6 +77,8 @@
#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 NGRIDY 20 /* number of grid point for grid of disks */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard in sub_wave.c */
@ -116,13 +97,9 @@
// #define T_IN 2.0 /* inside temperature */
#define SPEED 0.0 /* speed of drift to the right */
/* Boundary conditions */
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 0
#define BC_DIRICHLET 0 /* Dirichlet boundary conditions */
#define BC_PERIODIC 1 /* periodic boundary conditions */
#define BC_ABSORBING 2 /* absorbing boundary conditions (beta version) */
#define B_COND 1
/* Parameters for length and speed of simulation */
@ -142,26 +119,22 @@
/* Field representation */
#define FIELD_REP 0
#define FIELD_REP 1
#define F_INTENSITY 0 /* color represents intensity */
#define F_GRADIENT 1 /* color represents norm of gradient */
#define DRAW_FIELD_LINES 1 /* set to 1 to draw field lines */
#define FIELD_LINE_WIDTH 1 /* width of field lines */
#define N_FIELD_LINES 200 /* number of field lines */
#define N_FIELD_LINES 50 /* number of field lines */
#define FIELD_LINE_FACTOR 100 /* factor controlling precision when computing origin of field lines */
/* Color schemes */
/* Color schemes, see list in global_pdes.c */
#define BLACK 1 /* black background */
#define COLOR_SCHEME 1 /* choice of color scheme */
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
#define C_HUE 1 /* color scheme modifies hue */
#define C_PHASE 2 /* color scheme shows phase */
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
// #define SLOPE 0.1 /* sensitivity of color on wave amplitude */
#define SLOPE 0.3 /* sensitivity of color on wave amplitude */
@ -171,20 +144,13 @@
#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */
#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */
#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */
#define HUEMEAN 280.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -110.0 /* amplitude of variation of hue for color scheme C_HUE */
#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -220.0 /* amplitude of variation of hue for color scheme C_HUE */
// #define HUEMEAN 270.0 /* mean value of hue for color scheme C_HUE */
// #define HUEAMP -130.0 /* amplitude of variation of hue for color scheme C_HUE */
/* Basic math */
#define PI 3.141592654
#define DPI 6.283185307
#define PID 1.570796327
double julia_x = 0.0, julia_y = 0.0; /* parameters for Julia sets */
#include "global_pdes.c"
#include "sub_wave.c"
double courant2; /* Courant parameter squared */
@ -320,12 +286,7 @@ short int *xy_in[NX];
x1 = x1 + delta*nabx/norm;
y1 = y1 + delta*naby/norm;
}
else
{
cont = 0;
// nablax[ij[0]][ij[1]] = 0.0;
// nablay[ij[0]][ij[1]] = 0.0;
}
else cont = 0;
if (!xy_in[ij[0]][ij[1]]) cont = 0;
@ -361,7 +322,6 @@ int time;
}
/* compute the gradient */
// if (FIELD_REP > 0)
compute_gradient(phi, nablax, nablay);
/* compute the position of origins of field lines */
@ -371,25 +331,24 @@ int time;
printf("computing linex\n");
x1 = sqrt(3.58);
y1 = 0.0;
x1 = LAMBDA + MU*1.01;
y1 = 1.0;
linex[0] = x1;
liney[0] = y1;
dangle = DPI/((double)(N_FIELD_LINES*FIELD_LINE_FACTOR));
for (i = 1; i < N_FIELD_LINES*FIELD_LINE_FACTOR; i++)
{
// angle = PID + (double)i*dangle;
angle = (double)i*dangle;
x2 = sqrt(3.58)*cos(angle);
y2 = sqrt(1.18)*sin(angle);
x2 = LAMBDA + MU*1.01*cos(angle);
y2 = 0.5 + MU*1.01*sin(angle);
linex[i] = x2;
liney[i] = y2;
distance[i-1] = module2(x2-x1,y2-y1);
x1 = x2;
y1 = y2;
}
distance[N_FIELD_LINES*FIELD_LINE_FACTOR - 1] = module2(x2-sqrt(3.58),y2);
distance[N_FIELD_LINES*FIELD_LINE_FACTOR - 1] = module2(x2-LAMBDA,y2-0.5);
}
dx = (XMAX-XMIN)/((double)NX);
@ -404,7 +363,6 @@ int time;
value = module2(nablax[i][j], nablay[i][j]);
}
// if ((phi[i][j] - T_IN)*(phi[i][j] - T_OUT) < 0.0)
if (xy_in[i][j] == 1)
{
color_scheme(COLOR_SCHEME, value, scale, time, rgb);
@ -431,18 +389,14 @@ int time;
else integral[i] = intens;
}
deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES);
// deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES + 1.0);
// deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES);
// printf("delta = %.5lg\n", deltaintens);
i = 0;
// draw_field_line(linex[0], liney[0], nablax, nablay, 0.00002, 100000);
draw_field_line(linex[0], liney[0], xy_in, nablax, nablay, 0.00002, 100000);
for (j = 1; j < N_FIELD_LINES+1; j++)
{
while ((integral[i] <= j*deltaintens)&&(i < N_FIELD_LINES*FIELD_LINE_FACTOR)) i++;
// draw_field_line(linex[i], liney[i], nablax, nablay, 0.00002, 100000);
draw_field_line(linex[i], liney[i], xy_in, nablax, nablay, 0.00002, 100000);
counter++;
}
@ -459,12 +413,101 @@ int time;
void evolve_wave(phi, xy_in)
void evolve_wave_half(phi_in, phi_out, xy_in)
/* time step of field evolution */
double *phi_in[NX], *phi_out[NX]; short int *xy_in[NX];
{
int i, j, iplus, iminus, jplus, jminus;
double delta1, delta2, x, y;
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta1,delta2,x,y)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j] == 1){
/* discretized Laplacian depending on boundary conditions */
if ((B_COND == BC_DIRICHLET)||(B_COND == BC_ABSORBING))
{
iplus = (i+1); if (iplus == NX) iplus = NX-1;
iminus = (i-1); if (iminus == -1) iminus = 0;
jplus = (j+1); if (jplus == NY) jplus = NY-1;
jminus = (j-1); if (jminus == -1) jminus = 0;
}
else if (B_COND == BC_PERIODIC)
{
iplus = (i+1) % NX;
iminus = (i-1) % NX;
if (iminus < 0) iminus += NX;
jplus = (j+1) % NY;
jminus = (j-1) % NY;
if (jminus < 0) jminus += NY;
}
delta1 = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j];
x = phi_in[i][j];
/* evolve phi */
if (B_COND != BC_ABSORBING)
{
phi_out[i][j] = x + intstep*(delta1 - SPEED*(phi_in[iplus][j] - phi_in[i][j]));
}
else /* case of absorbing b.c. - this is only an approximation of correct way of implementing */
{
/* in the bulk */
if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
{
phi_out[i][j] = x - intstep*delta2;
}
/* right border */
else if (i==NX-1)
{
phi_out[i][j] = x - intstep1*(x - phi_in[i-1][j]);
}
/* upper border */
else if (j==NY-1)
{
phi_out[i][j] = x - intstep1*(x - phi_in[i][j-1]);
}
/* left border */
else if (i==0)
{
phi_out[i][j] = x - intstep1*(x - phi_in[1][j]);
}
/* lower border */
else if (j==0)
{
phi_out[i][j] = x - intstep1*(x - phi_in[i][1]);
}
}
if (FLOOR)
{
if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
}
}
}
}
// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
}
void evolve_wave(phi, phi_tmp, xy_in)
/* time step of field evolution */
double *phi[NX], *phi_tmp[NX]; short int *xy_in[NX];
{
evolve_wave_half(phi, phi_tmp, xy_in);
evolve_wave_half(phi_tmp, phi, xy_in);
}
void old_evolve_wave(phi, xy_in)
/* time step of field evolution */
double *phi[NX]; short int *xy_in[NX];
{
int i, j, iplus, iminus, jplus, jminus;
double delta1, delta2, x, y, *newphi[NX];;
double delta1, delta2, x, y, *newphi[NX];
for (i=0; i<NX; i++) newphi[i] = (double *)malloc(NY*sizeof(double));
@ -552,9 +595,6 @@ void evolve_wave(phi, xy_in)
// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
}
double compute_variance(phi, xy_in)
/* compute the variance (total probability) of the field */
double *phi[NX]; short int * xy_in[NX];
@ -658,7 +698,6 @@ short int *xy_in[NX];
sinj = sin(jangle);
julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj);
julia_y = 0.5*sinj*(1.0-cosj) + yshift;
/* need to decrease 0.05 for i > 2000 */
// julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj);
// julia_y = 0.5*sinj*(1.0-cosj);
init_julia_set(phi, xy_in);
@ -669,7 +708,7 @@ short int *xy_in[NX];
void animation()
{
double time, scale, dx, var, jangle, cosj, sinj;
double *phi[NX];
double *phi[NX], *phi_tmp[NX];
short int *xy_in[NX];
int i, j, s;
@ -677,6 +716,7 @@ void animation()
for (i=0; i<NX; i++)
{
phi[i] = (double *)malloc(NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
@ -687,7 +727,7 @@ void animation()
// julia_x = 0.1;
// julia_y = 0.6;
set_Julia_parameters(0, phi, xy_in);
// set_Julia_parameters(0, phi, xy_in);
printf("Integration step %.3lg\n", intstep);
@ -710,7 +750,7 @@ void animation()
draw_wave(phi, xy_in, 1.0, 0);
draw_billiard();
print_Julia_parameters(i);
// print_Julia_parameters(i);
// print_level(MDEPTH);
@ -734,17 +774,17 @@ void animation()
draw_wave(phi, xy_in, scale, i);
for (j=0; j<NVID; j++) evolve_wave(phi, xy_in);
for (j=0; j<NVID; j++) evolve_wave(phi, phi_tmp, xy_in);
draw_billiard();
// print_level(MDEPTH);
print_Julia_parameters(i);
// print_Julia_parameters(i);
glutSwapBuffers();
/* modify Julia set */
set_Julia_parameters(i, phi, xy_in);
// set_Julia_parameters(i, phi, xy_in);
if (MOVIE)
{
@ -770,6 +810,7 @@ void animation()
for (i=0; i<NX; i++)
{
free(phi[i]);
free(phi_tmp[i]);
}
}

215
particle_billiard.c
View File

@ -47,23 +47,20 @@
// #define YMIN -0.91
// #define YMAX 1.115 /* y interval for 9/16 aspect ratio */
/* Choice of the billiard table */
/* Choice of the billiard table, see global_particles.c */
#define B_DOMAIN 9 /* choice of domain shape */
#define B_DOMAIN 14 /* choice of domain shape */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_STADIUM 2 /* stadium-shaped domain */
#define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular billiard */
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_REULEAUX 9 /* Reuleaux and star shapes */
#define D_FLOWER 10 /* Bunimovich flower */
#define D_ALT_REU 11 /* alternating between star and Reuleaux */
#define CIRCLE_PATTERN 0 /* pattern of circles */
#define LAMBDA -3.346065215 /* sin(60°)/sin(15°) for Reuleaux-type triangle with 90° angles */
#define NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */
// #define NCX 10 /* number of circles in x direction */
// #define NCY 15 /* number of circles in y direction */
#define NCX 15 /* number of circles in x direction */
#define NCY 20 /* number of circles in y direction */
#define LAMBDA 0.75 /* parameter controlling shape of billiard */
// #define LAMBDA -3.346065215 /* sin(60°)/sin(15°) for Reuleaux-type triangle with 90° angles */
// #define LAMBDA 3.0 /* parameter controlling shape of billiard */
// #define LAMBDA 0.6 /* parameter controlling shape of billiard */
// #define LAMBDA 0.4175295 /* sin(20°)/sin(55°) for 9-star shape with 30° angles */
@ -71,28 +68,30 @@
// #define LAMBDA 3.75738973 /* sin(36°)/sin(9°) for 5-star shape with 90° angles */
// #define LAMBDA -1.73205080756888 /* -sqrt(3) for Reuleaux triangle */
// #define LAMBDA 1.73205080756888 /* sqrt(3) for triangle tiling plane */
#define MU 0.1 /* second parameter controlling shape of billiard */
#define MU 0.035 /* second parameter controlling shape of billiard */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */
#define NPOLY 8 /* number of sides of polygon */
#define APOLY 0.25 /* angle by which to turn polygon, in units of Pi/2 */
#define DRAW_BILLIARD 1 /* set to 1 to draw billiard */
#define DRAW_CONSTRUCTION_LINES 1 /* set to 1 to draw additional construction lines for billiard */
#define PERIODIC_BC 0 /* set to 1 to enforce periodic boundary conditions when drawing particles */
#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */
#define DEBUG 0 /* draw trajectories, for debugging purposes */
/* Simulation parameters */
#define NPART 20000 /* number of particles */
#define NPART 5000 /* number of particles */
#define NPARTMAX 100000 /* maximal number of particles after resampling */
#define LMAX 0.01 /* minimal segment length triggering resampling */
#define DMIN 0.02 /* minimal distance to boundary for triggering resampling */
#define CYCLE 1 /* set to 1 for closed curve (start in all directions) */
#define NSTEPS 6000 /* number of frames of movie */
#define NSTEPS 3000 /* number of frames of movie */
#define TIME 1000 /* time between movie frames, for fluidity of real-time simulation */
// #define DPHI 0.000005 /* integration step */
#define DPHI 0.00002 /* integration step */
// #define DPHI 0.000002 /* integration step */
// #define DPHI 0.00002 /* integration step */
#define DPHI 0.000005 /* integration step */
#define NVID 150 /* number of iterations between images displayed on screen */
/* Decreasing TIME accelerates the animation and the movie */
@ -103,12 +102,13 @@
/* Colors and other graphical parameters */
#define NCOLORS -10 /* number of colors */
#define COLORSHIFT 200 /* hue of initial color */
#define NCOLORS 16 /* number of colors */
#define COLORSHIFT 0 /* hue of initial color */
#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */
#define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */
#define NSEG 100 /* number of segments of boundary */
#define LENGTH 0.04 /* length of velocity vectors */
#define BILLIARD_WIDTH 3 /* width of billiard */
#define LENGTH 0.02 /* length of velocity vectors */
#define BILLIARD_WIDTH 2 /* width of billiard */
#define PARTICLE_WIDTH 2 /* width of particles */
#define FRONT_WIDTH 3 /* width of wave front */
@ -123,10 +123,8 @@
#define SLEEP1 1 /* initial sleeping time */
#define SLEEP2 1000 /* final sleeping time */
#define PI 3.141592654
#define DPI 6.283185307
#define PID 1.570796327
#include "global_particles.c"
#include "sub_part_billiard.c"
@ -169,7 +167,8 @@ double *configs[NPARTMAX];
double conf[2], pos[2];
while (angle2 < angle1) angle2 += DPI;
dalpha = (angle2 - angle1)/((double)(NPART-1));
if (NPART > 1) dalpha = (angle2 - angle1)/((double)(NPART-1));
else dalpha = 0.0;
for (i=0; i<NPART; i++)
{
alpha = angle1 + dalpha*((double)i);
@ -253,8 +252,11 @@ double *configs[NPARTMAX];
if (configs[i][2]<0.0)
{
vbilliard(configs[i]);
color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS;
if (!RAINBOW_COLOR)
{
color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS;
}
}
configs[i][2] += DPHI;
@ -280,36 +282,39 @@ double *configs[NPARTMAX];
glEnd ();
/* taking care of boundary conditions - only needed for periodic boundary conditions */
if (SCALING_FACTOR*x2 > XMAX)
if (PERIODIC_BC)
{
glBegin(GL_LINE_STRIP);
glVertex2d(SCALING_FACTOR*(x1+XMIN-XMAX), SCALING_FACTOR*y1);
glVertex2d(SCALING_FACTOR*(x2+XMIN-XMAX), SCALING_FACTOR*y2);
glEnd ();
}
if (SCALING_FACTOR*x2 > XMAX)
{
glBegin(GL_LINE_STRIP);
glVertex2d(SCALING_FACTOR*(x1+XMIN-XMAX), SCALING_FACTOR*y1);
glVertex2d(SCALING_FACTOR*(x2+XMIN-XMAX), SCALING_FACTOR*y2);
glEnd ();
}
if (SCALING_FACTOR*x2 < XMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(SCALING_FACTOR*(x1-XMIN+XMAX), SCALING_FACTOR*y1);
glVertex2d(SCALING_FACTOR*(x2-XMIN+XMAX), SCALING_FACTOR*y2);
glEnd ();
}
if (SCALING_FACTOR*x2 < XMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(SCALING_FACTOR*(x1-XMIN+XMAX), SCALING_FACTOR*y1);
glVertex2d(SCALING_FACTOR*(x2-XMIN+XMAX), SCALING_FACTOR*y2);
glEnd ();
}
if (SCALING_FACTOR*y2 > YMAX)
{
glBegin(GL_LINE_STRIP);
glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMIN-YMAX));
glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMIN-YMAX));
glEnd ();
}
if (SCALING_FACTOR*y2 > YMAX)
{
glBegin(GL_LINE_STRIP);
glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMIN-YMAX));
glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMIN-YMAX));
glEnd ();
}
if (SCALING_FACTOR*y2 < YMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMAX-YMIN));
glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMAX-YMIN));
glEnd ();
if (SCALING_FACTOR*y2 < YMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMAX-YMIN));
glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMAX-YMIN));
glEnd ();
}
}
}
@ -343,11 +348,14 @@ double *configs[NPARTMAX];
{
if (configs[i][2]<0.0)
{
// printf("reflecting particle %i\n", i);
c = vbilliard(configs[i]);
// if (c>=0) color[i]++;
color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS;
if (!RAINBOW_COLOR)
{
color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS;
}
}
configs[i][2] += DPHI;
@ -363,7 +371,74 @@ double *configs[NPARTMAX];
}
void init_circle_config()
{
int i, j, n;
double dx, dy;
switch (CIRCLE_PATTERN) {
case (C_FOUR_CIRCLES):
{
ncircles = 4;
circlex[0] = 1.0;
circley[0] = 0.0;
circlerad[0] = 0.8;
circlex[1] = -1.0;
circley[1] = 0.0;
circlerad[1] = 0.8;
circlex[2] = 0.0;
circley[2] = 0.8;
circlerad[2] = 0.4;
circlex[3] = 0.0;
circley[3] = -0.8;
circlerad[3] = 0.4;
for (i=0; i<4; i++) circleactive[i] = 1;
break;
}
case (C_SQUARE):
{
ncircles = NCX*NCY;
dy = (YMAX - YMIN)/((double)NCY);
for (i = 0; i < NCX; i++)
for (j = 0; j < NCY; j++)
{
n = NCY*i + j;
circlex[n] = ((double)(i-NCX/2) + 0.5)*dy;
circley[n] = YMIN + ((double)j + 0.5)*dy;
circlerad[n] = MU;
circleactive[n] = 1;
}
break;
}
case (C_HEX):
{
ncircles = NCX*(NCY+1);
dy = (YMAX - YMIN)/((double)NCY);
dx = dy*0.5*sqrt(3.0);
for (i = 0; i < NCX; i++)
for (j = 0; j < NCY+1; j++)
{
n = (NCY+1)*i + j;
circlex[n] = ((double)(i-NCX/2) + 0.5)*dy;
circley[n] = YMIN + ((double)j - 0.5)*dy;
if ((i+NCX)%2 == 1) circley[n] += 0.5*dy;
circlerad[n] = MU;
circleactive[n] = 1;
}
break;
}
default:
{
printf("Function init_circle_config not defined for this pattern \n");
}
}
}
void animation()
@ -379,19 +454,24 @@ void animation()
active = malloc(sizeof(int)*(NPARTMAX));
for (i=0; i<NPARTMAX; i++)
configs[i] = (double *)malloc(8*sizeof(double));
/* init circle configuration if the domain is D_CIRCLES */
if (B_DOMAIN == D_CIRCLES) init_circle_config();
/* initialize system by putting particles in a given point with a range of velocities */
r = cos(PI/(double)NPOLY)/cos(DPI/(double)NPOLY);
// init_drop_config(0.0, 0.0, -0.2, 0.2, configs);
// init_drop_config(0.0, 0.0, 0.0, PI, configs);
// init_drop_config(0.5, 0.5, -1.0, 1.0, configs);
// init_sym_drop_config(-1.0, 0.5, -PID, PID, configs);
// init_drop_config(-0.999, 0.0, -alpha, alpha, configs);
// other possible initial conditions :
// init_line_config(-0.6, 0.2, -0.6, 0.7, 0.0, configs);
// init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0, configs);
init_line_config(0.0, -1.0, -1.0, 1.0, 0.25*PID, configs);
// init_line_config(-0.7, -0.45, -0.7, 0.45, 0.0, configs);
init_line_config(0.0, -0.3, 0.0, 0.3, 0.0, configs);
// init_line_config(-1.5, 0.1, -0.1, 1.0, -0.5*PID, configs);
blank();
glColor3f(0.0, 0.0, 0.0);
@ -425,7 +505,14 @@ void animation()
}
}
sleep(SLEEP1);
if (RAINBOW_COLOR) /* rainbow color scheme */
for (i=0; i<NPART; i++)
{
color[i] = (i*NCOLORS)/NPART;
newcolor[i] = (i*NCOLORS)/NPART;
}
sleep(SLEEP1);
for (i=0; i<=NSTEPS; i++)
{

179
schrodinger.c
View File

@ -42,8 +42,10 @@
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
#define NX 640 /* number of grid points on x axis */
#define NY 360 /* number of grid points on y axis */
#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 */
@ -55,67 +57,45 @@
#define JULIA_SCALE 1.0 /* scaling for Julia sets */
/* Choice of the billiard table */
/* Choice of the billiard table, see list in global_pdes.c */
#define B_DOMAIN 3 /* choice of domain shape */
#define B_DOMAIN 9 /* choice of domain shape */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_STADIUM 2 /* stadium-shaped domain */
#define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular billiard */
#define D_FLAT 6 /* flat interface */
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_YOUNG 9 /* Young diffraction slits */
#define D_GRATING 10 /* diffraction grating */
#define D_EHRENFEST 11 /* Ehrenfest urn type geometry */
#define CIRCLE_PATTERN 0 /* pattern of circles, see list in global_pdes.c */
#define D_MENGER 15 /* Menger-Sierpinski carpet */
#define D_JULIA_INT 16 /* interior of Julia set */
#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 */
/* Billiard tables for heat equation */
#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
#define D_MENGER_HEATED 22 /* Menger gasket with different temperatures */
#define D_MENGER_H_OPEN 23 /* Menger gasket with different temperatures and larger domain */
#define D_MANDELBROT 24 /* Mandelbrot set */
#define D_JULIA 25 /* Julia set */
#define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */
#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */
#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */
#define MU 0.05 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 2 /* depth of computation of Menger gasket */
#define MRATIO 5 /* ratio defining Menger gasket */
#define MDEPTH 3 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */
#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 NGRIDY 20 /* number of grid point for grid of disks */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard in sub_wave.c */
/* Physical patameters of wave equation */
#define DT 0.00000005
// #define DT 0.00000002
// #define DT 0.00000005
#define DT 0.000000005
// #define DT 0.000000005
#define HBAR 1.0
/* Boundary conditions */
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 1
#define BC_DIRICHLET 0 /* Dirichlet boundary conditions */
#define BC_PERIODIC 1 /* periodic boundary conditions */
#define BC_ABSORBING 2 /* absorbing boundary conditions (beta version) */
/* Parameters for length and speed of simulation */
#define NSTEPS 4500 /* number of frames of movie */
#define NVID 250 /* number of iterations between images displayed on screen */
#define NSTEPS 200 /* number of frames of movie */
#define NVID 1200 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
#define PAUSE 1000 /* number of frames after which to pause */
@ -128,25 +108,17 @@
#define VMAX 10.0 /* max value of wave amplitude */
/* Plot type */
/* Plot type, see list in global_pdes.c */
#define PLOT 0
#define PLOT 10
#define P_MODULE 0 /* plot module of wave function squared */
#define P_PHASE 1 /* plot phase of wave function */
#define P_REAL 2 /* plot real part */
#define P_IMAGINARY 3 /* plot imaginary part */
/* Color schemes */
/* Color schemes, see list in global_pdes.c */
#define BLACK 1 /* black background */
#define COLOR_SCHEME 1 /* choice of color scheme */
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
#define C_HUE 1 /* color scheme modifies hue */
#define C_PHASE 2 /* color scheme shows phase */
#define SCALE 1 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
@ -158,14 +130,7 @@
#define HUEMEAN 150.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -150.0 /* amplitude of variation of hue for color scheme C_HUE */
/* Basic math */
#define PI 3.141592654
#define DPI 6.283185307
#define PID 1.570796327
double julia_x = 0.0, julia_y = 0.0; /* parameters for Julia sets */
#include "global_pdes.c"
#include "sub_wave.c"
double courant2; /* Courant parameter squared */
@ -269,22 +234,14 @@ int time;
glEnd ();
}
void evolve_wave(phi, psi, xy_in)
void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in)
/* time step of field evolution */
/* phi is real part, psi is imaginary part */
double *phi[NX], *psi[NX];
short int *xy_in[NX];
double *phi_in[NX], *psi_in[NX], *phi_out[NX], *psi_out[NX]; short int *xy_in[NX];
{
int i, j, iplus, iminus, jplus, jminus;
double delta1, delta2, x, y, *newphi[NX], *newpsi[NX];
double delta1, delta2, x, y;
for (i=0; i<NX; i++)
{
newphi[i] = (double *)malloc(NY*sizeof(double));
newpsi[i] = (double *)malloc(NY*sizeof(double));
}
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta1,delta2,x,y)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
@ -307,81 +264,76 @@ void evolve_wave(phi, psi, xy_in)
if (jminus < 0) jminus += NY;
}
delta1 = phi[iplus][j] + phi[iminus][j] + phi[i][jplus] + phi[i][jminus] - 4.0*phi[i][j];
delta2 = psi[iplus][j] + psi[iminus][j] + psi[i][jplus] + psi[i][jminus] - 4.0*psi[i][j];
delta1 = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j];
delta2 = psi_in[iplus][j] + psi_in[iminus][j] + psi_in[i][jplus] + psi_in[i][jminus] - 4.0*psi_in[i][j];
x = phi[i][j];
y = psi[i][j];
x = phi_in[i][j];
y = psi_in[i][j];
/* evolve phi and psi */
if (B_COND != BC_ABSORBING)
{
newphi[i][j] = x - intstep*delta2;
newpsi[i][j] = y + intstep*delta1;
phi_out[i][j] = x - intstep*delta2;
psi_out[i][j] = y + intstep*delta1;
}
else /* case of absorbing b.c. - this is only an approximation of correct way of implementing */
{
/* in the bulk */
if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
{
newphi[i][j] = x - intstep*delta2;
newpsi[i][j] = y + intstep*delta1;
phi_out[i][j] = x - intstep*delta2;
psi_out[i][j] = y + intstep*delta1;
}
/* right border */
else if (i==NX-1)
{
newphi[i][j] = x - intstep1*(y - psi[i-1][j]);
newpsi[i][j] = y + intstep1*(x - phi[i-1][j]);
phi_out[i][j] = x - intstep1*(y - psi_in[i-1][j]);
psi_out[i][j] = y + intstep1*(x - phi_in[i-1][j]);
}
/* upper border */
else if (j==NY-1)
{
newphi[i][j] = x - intstep1*(y - psi[i][j-1]);
newpsi[i][j] = y + intstep1*(x - phi[i][j-1]);
phi_out[i][j] = x - intstep1*(y - psi_in[i][j-1]);
psi_out[i][j] = y + intstep1*(x - phi_in[i][j-1]);
}
/* left border */
else if (i==0)
{
newphi[i][j] = x - intstep1*(y - psi[1][j]);
newpsi[i][j] = y + intstep1*(x - phi[1][j]);
phi_out[i][j] = x - intstep1*(y - psi_in[1][j]);
psi_out[i][j] = y + intstep1*(x - phi_in[1][j]);
}
/* lower border */
else if (j==0)
{
newphi[i][j] = x - intstep1*(y - psi[i][1]);
newpsi[i][j] = y + intstep1*(x - phi[i][1]);
phi_out[i][j] = x - intstep1*(y - psi_in[i][1]);
psi_out[i][j] = y + intstep1*(x - phi_in[i][1]);
}
}
if (FLOOR)
{
if (newphi[i][j] > VMAX) newphi[i][j] = VMAX;
if (newphi[i][j] < -VMAX) newphi[i][j] = -VMAX;
if (newpsi[i][j] > VMAX) newpsi[i][j] = VMAX;
if (newpsi[i][j] < -VMAX) newpsi[i][j] = -VMAX;
if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX;
if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX;
}
}
}
}
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j] == 1) phi[i][j] = newphi[i][j];
if (xy_in[i][j] == 1) psi[i][j] = newpsi[i][j];
}
}
for (i=0; i<NX; i++)
{
free(newphi[i]);
free(newpsi[i]);
}
// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
}
void evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in)
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX]; short int *xy_in[NX];
{
evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in);
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
}
double compute_variance(phi, psi, xy_in)
/* compute the variance (total probability) of the field */
@ -428,7 +380,7 @@ short int * xy_in[NX];
void animation()
{
double time, scale, dx, var;
double *phi[NX], *psi[NX];
double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX];
short int *xy_in[NX];
int i, j, s;
@ -437,6 +389,8 @@ void animation()
{
phi[i] = (double *)malloc(NY*sizeof(double));
psi[i] = (double *)malloc(NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NY*sizeof(double));
psi_tmp[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
@ -447,7 +401,7 @@ void animation()
printf("Integration step %.3lg\n", intstep);
/* initialize wave wave function */
init_coherent_state(-1.2, 0.0, 20.0, 0.0, 0.2, phi, psi, xy_in);
init_coherent_state(-1.2, 0.0, 20.0, 0.0, 0.25, phi, psi, xy_in);
// init_coherent_state(0.0, 0.0, 0.0, 5.0, 0.03, phi, psi, xy_in);
// init_coherent_state(-0.5, 0.0, 1.0, 1.0, 0.05, phi, psi, xy_in);
@ -488,7 +442,7 @@ void animation()
// printf("Wave drawn\n");
for (j=0; j<NVID; j++) evolve_wave(phi, psi, xy_in);
for (j=0; j<NVID; j++) evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
draw_billiard();
@ -510,15 +464,18 @@ void animation()
}
if (MOVIE)
{
for (i=0; i<20; i++) save_frame();
s = system("mv wave*.tif tif_schrod/");
}
// if (MOVIE)
// {
// for (i=0; i<20; i++) save_frame();
// s = system("mv wave*.tif tif_schrod/");
// }
for (i=0; i<NX; i++)
{
free(phi[i]);
free(psi[i]);
free(phi_tmp[i]);
free(psi_tmp[i]);
free(xy_in[i]);
}
}

294
sub_wave.c
View File

@ -315,13 +315,129 @@ double x1, y1, x2, y2;
}
void draw_rotated_rectangle(x1, y1, x2, y2)
double x1, y1, x2, y2;
{
double pos[2];
double xa, ya, xb, yb, xc, yc;
xa = 0.5*(x1 - y2);
xb = 0.5*(x2 - y1);
xc = 0.5*(x1 - y1);
ya = 0.5*(x1 + y1);
yb = 0.5*(x2 + y2);
yc = 0.5*(x2 + y1);
glBegin(GL_LINE_LOOP);
xy_to_pos(xc, ya, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(xb, yc, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(xc, yb, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(xa, yc, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
}
void init_circle_config()
/* initialise the arrays circlex, circley, circlerad and circleactive */
/* for billiard shape D_CIRCLES */
{
int i, j, n;
double dx, dy, p;
switch (CIRCLE_PATTERN) {
case (C_SQUARE):
{
ncircles = NGRIDX*NGRIDY;
dy = (YMAX - YMIN)/((double)NGRIDY);
for (i = 0; i < NGRIDX; i++)
for (j = 0; j < NGRIDY; j++)
{
n = NGRIDY*i + j;
circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy;
circley[n] = YMIN + ((double)j + 0.5)*dy;
circlerad[n] = MU;
circleactive[n] = 1;
}
break;
}
case (C_HEX):
{
ncircles = NGRIDX*(NGRIDY+1);
dy = (YMAX - YMIN)/((double)NGRIDY);
dx = dy*0.5*sqrt(3.0);
for (i = 0; i < NGRIDX; i++)
for (j = 0; j < NGRIDY+1; j++)
{
n = (NGRIDY+1)*i + j;
circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy;
circley[n] = YMIN + ((double)j - 0.5)*dy;
if ((i+NGRIDX)%2 == 1) circley[n] += 0.5*dy;
circlerad[n] = MU;
circleactive[n] = 1;
}
break;
}
case (C_RAND_DISPLACED):
{
ncircles = NGRIDX*NGRIDY;
dy = (YMAX - YMIN)/((double)NGRIDY);
for (i = 0; i < NGRIDX; i++)
for (j = 0; j < NGRIDY; j++)
{
n = NGRIDY*i + j;
circlex[n] = ((double)(i-NGRIDX/2) + 0.5 + 0.5*(double)rand()/RAND_MAX)*dy;
circley[n] = YMIN + ((double)j + 0.5 + 0.5*(double)rand()/RAND_MAX)*dy;
circlerad[n] = MU*(1.0 + 0.35*(double)rand()/RAND_MAX);
circleactive[n] = 1;
}
break;
}
case (C_RAND_PERCOL):
{
ncircles = NGRIDX*NGRIDY;
dy = (YMAX - YMIN)/((double)NGRIDY);
for (i = 0; i < NGRIDX; i++)
for (j = 0; j < NGRIDY; j++)
{
n = NGRIDY*i + j;
circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy;
circley[n] = YMIN + ((double)j + 0.5)*dy;
circlerad[n] = MU;
p = (double)rand()/RAND_MAX;
if (p < P_PERCOL) circleactive[n] = 1;
else circleactive[n] = 0;
}
break;
}
case (C_RAND_POISSON):
{
ncircles = NPOISSON;
for (n = 0; n < NPOISSON; n++)
{
circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0);
circley[n] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN;
circlerad[n] = MU;
circleactive[n] = 1;
}
break;
}
default:
{
printf("Function init_circle_config not defined for this pattern \n");
}
}
}
int xy_in_billiard(x, y)
/* returns 1 if (x,y) represents a point in the billiard */
double x, y;
{
double l2, r2, r2mu, omega, c, angle, z, x1, y1, u, v, u1, v1;
int i, k, k1, k2, condition;
double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy;
int i, j, k, k1, k2, condition;
switch (B_DOMAIN) {
case D_RECTANGLE:
@ -416,6 +532,44 @@ double x, y;
{
return(((x-1.0)*(x-1.0) + y*y < LAMBDA*LAMBDA)||((x+1.0)*(x+1.0) + y*y < LAMBDA*LAMBDA)||((vabs(x) < 1.0)&&(vabs(y) < MU)));
}
case D_DISK_GRID:
{
dy = (YMAX - YMIN)/((double)NGRIDY);
for (i = -NGRIDX/2; i < NGRIDX/2; i++)
for (j = 0; j < NGRIDY; j++)
{
x1 = ((double)i + 0.5)*dy;
y1 = YMIN + ((double)j + 0.5)*dy;
if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < MU*MU) return(0);
}
return(1);
}
case D_DISK_HEX:
{
dy = (YMAX - YMIN)/((double)NGRIDY);
dx = dy*0.5*sqrt(3.0);
for (i = -NGRIDX/2; i < NGRIDX/2; i++)
for (j = -1; j < NGRIDY; j++)
{
x1 = ((double)i + 0.5)*dy;
y1 = YMIN + ((double)j + 0.5)*dy;
if ((i+NGRIDX)%2 == 1) y1 += 0.5*dy;
if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < MU*MU) return(0);
}
return(1);
}
case D_CIRCLES:
{
for (i = 0; i < ncircles; i++)
if (circleactive[i])
{
x1 = circlex[i];
y1 = circley[i];
r2 = circlerad[i]*circlerad[i];
if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(0);
}
return(1);
}
case D_MENGER:
{
x1 = 0.5*(x+1.0);
@ -443,6 +597,23 @@ double x, y;
if (u*u + v*v < MANDELLIMIT) return(1);
else return(0);
}
case D_MENGER_ROTATED:
{
x2 = 1.0*(x + y);
y2 = 1.0*(x - y);
if ((vabs(x2) < 1.0)&&(vabs(y2) < 1.0))
{
x1 = 0.5*(x2 + 1.0);
y1 = 0.5*(y2 + 1.0);
for (k=0; k<MDEPTH; k++)
{
x1 = x1*(double)MRATIO;
y1 = y1*(double)MRATIO;
if ((vabs(x)<1.0)&&(vabs(y)<1.0)&&(((int)x1 % MRATIO)==MRATIO/2)&&(((int)y1 % MRATIO)==MRATIO/2)) return(0);
}
}
return(1);
}
case D_ANNULUS_HEATED: /* returns 2 if in inner circle */
{
l2 = LAMBDA*LAMBDA;
@ -558,8 +729,8 @@ int i, j;
void draw_billiard() /* draws the billiard boundary */
{
double x0, x, y, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l;
int i, j, k1, k2, mr2;
double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l;
int i, j, k, k1, k2, mr2;
if (BLACK) glColor3f(1.0, 1.0, 1.0);
else glColor3f(0.0, 0.0, 0.0);
@ -845,6 +1016,71 @@ void draw_billiard() /* draws the billiard boundary */
glEnd ();
break;
}
case D_DISK_GRID:
{
glLineWidth(2);
for (i = -NGRIDX/2; i < NGRIDX/2; i++)
for (j = 0; j < NGRIDY; j++)
{
dy = (YMAX - YMIN)/((double)NGRIDY);
dx = dy*0.5*sqrt(3.0);
x1 = ((double)i + 0.5)*dy;
y1 = YMIN + ((double)j + 0.5)*dy;
glBegin(GL_LINE_LOOP);
for (k=0; k<=NSEG; k++)
{
phi = (double)k*DPI/(double)NSEG;
x = x1 + MU*cos(phi);
y = y1 + MU*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
}
break;
}
case D_DISK_HEX:
{
glLineWidth(2);
for (i = -NGRIDX/2; i < NGRIDX/2; i++)
for (j = -1; j < NGRIDY; j++)
{
dy = (YMAX - YMIN)/((double)NGRIDY);
x1 = ((double)i + 0.5)*dy;
y1 = YMIN + ((double)j + 0.5)*dy;
if ((i+NGRIDX)%2 == 1) y1 += 0.5*dy;
glBegin(GL_LINE_LOOP);
for (k=0; k<=NSEG; k++)
{
phi = (double)k*DPI/(double)NSEG;
x = x1 + MU*cos(phi);
y = y1 + MU*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
}
break;
}
case (D_CIRCLES):
{
glLineWidth(2);
for (i = 0; i < ncircles; i++)
if (circleactive[i])
{
glBegin(GL_LINE_LOOP);
for (k=0; k<=NSEG; k++)
{
phi = (double)k*DPI/(double)NSEG;
x = circlex[i] + circlerad[i]*cos(phi);
y = circley[i] + circlerad[i]*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
}
break;
}
case (D_MENGER):
{
glLineWidth(3);
@ -893,11 +1129,59 @@ void draw_billiard() /* draws the billiard boundary */
break;
}
case (D_JULIA_INT):
case (D_JULIA_INT):
{
/* Do nothing */
break;
}
case (D_MENGER_ROTATED):
{
glLineWidth(3);
// draw_rectangle(XMIN, -1.0, XMAX, 1.0);
/* level 1 */
if (MDEPTH > 0)
{
glLineWidth(2);
x = 1.0/((double)MRATIO);
draw_rotated_rectangle(x, x, -x, -x);
}
/* level 2 */
if (MDEPTH > 1)
{
glLineWidth(1);
mr2 = MRATIO*MRATIO;
l = 2.0/((double)mr2);
for (i=0; i<MRATIO; i++)
for (j=0; j<MRATIO; j++)
if ((i!=MRATIO/2)||(j!=MRATIO/2))
{
x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)MRATIO);
y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)MRATIO);
draw_rotated_rectangle(x, y, x+l, y+l);
}
}
/* level 3 */
if (MDEPTH > 2)
{
glLineWidth(1);
l = 2.0/((double)(mr2*MRATIO));
for (i=0; i<mr2; i++)
for (j=0; j<mr2; j++)
if ( (((i%MRATIO!=MRATIO/2))||(j%MRATIO!=MRATIO/2)) && (((i/MRATIO!=MRATIO/2))||(j/MRATIO!=MRATIO/2)) )
{
x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)mr2);
y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)mr2);
draw_rotated_rectangle(x, y, x+l, y+l);
}
}
break;
}
case (D_ANNULUS_HEATED):
{
glBegin(GL_LINE_LOOP);

288
wave_billiard.c
View File

@ -58,130 +58,105 @@
#define XMAX 2.0 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
// #define XMIN -1.8
// #define XMAX 1.8 /* x interval */
// #define YMIN -1.0125
// #define YMAX 1.0125 /* y interval for 9/16 aspect ratio */
#define JULIA_SCALE 1.1 /* scaling for Julia sets */
#define JULIA_SCALE 1.0 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 16 /* choice of domain shape */
#define B_DOMAIN 3 /* choice of domain shape, see list in global_pdes.c */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
#define D_STADIUM 2 /* stadium-shaped domain */
#define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular billiard */
#define D_FLAT 6 /* flat interface */
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_YOUNG 9 /* Young diffraction slits */
#define D_GRATING 10 /* diffraction grating */
#define D_EHRENFEST 11 /* Ehrenfest urn type geometry */
#define CIRCLE_PATTERN 4 /* pattern of circles, see list in global_pdes.c */
#define D_MENGER 15 /* Menger-Sierpinski carpet */
#define D_JULIA_INT 16 /* interior of Julia set */
#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 */
/* Billiard tables for heat equation */
#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
#define D_MENGER_HEATED 22 /* Menger gasket with different temperatures */
#define D_MENGER_H_OPEN 23 /* Menger gasket with different temperatures and larger domain */
#define D_MANDELBROT 24 /* Mandelbrot set */
#define D_JULIA 25 /* Julia set */
#define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */
#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define MU 0.05 /* parameter controlling the dimensions of domain */
#define NPOLY 8 /* number of sides of polygon */
#define LAMBDA 0.75 /* parameter controlling the dimensions of domain */
#define MU 0.025 /* 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 */
#define MRATIO 3 /* ratio defining Menger gasket */
#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 NGRIDY 20 /* number of grid point for grid of disks */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard below */
/* Physical parameters of wave equation */
#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */
#define OMEGA 0.9 /* frequency of periodic excitation */
#define COURANT 0.01 /* Courant number */
#define COURANTB 0.003 /* Courant number in medium B */
#define GAMMA 0.0 /* damping factor in wave equation */
#define GAMMA_BC 1.0e-4 /* damping factor on boundary */
// #define GAMMA 5.0e-10 /* damping factor in wave equation */
#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */
// #define KAPPA 5.0e-6 /* "elasticity" term enforcing oscillations */
// #define KAPPA 5.0e-9 /* "elasticity" term enforcing oscillations */
// #define KAPPA 5.0e-8 /* "elasticity" term enforcing oscillations */
#define KAPPA_BC 5.0e-4 /* "elasticity" term on absorbing boundary */
/* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
/* The physical damping coefficient is given by GAMMA/(DT)^2 */
/* Increasing COURANT speeds up the simulation, but decreases accuracy */
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
/* Boundary conditions */
/* Boundary conditions, see list in global_pdes.c */
#define B_COND 2
#define BC_DIRICHLET 0 /* Dirichlet boundary conditions */
#define BC_PERIODIC 1 /* periodic boundary conditions */
#define BC_ABSORBING 2 /* absorbing boundary conditions (beta version) */
/* 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 B_COND 3
/* Parameters for length and speed of simulation */
#define NSTEPS 4000 /* number of frames of movie */
#define NVID 25 /* number of iterations between images displayed on screen */
#define NSTEPS 7000 /* number of frames of movie */
#define NVID 50 /* 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 */
#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 */
/* Plot type, see list in global_pdes.c */
#define PLOT 1
/* Color schemes */
#define BLACK 1 /* background */
#define COLOR_SCHEME 1 /* choice of color scheme */
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
#define C_HUE 1 /* color scheme modifies hue */
#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.0 /* sensitivity of color on wave amplitude */
// #define SLOPE 0.5 /* sensitivity of color on wave amplitude */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
#define E_SCALE 750.0 /* scaling factor for energy representation */
#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */
#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */
#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */
#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */
#define HUEMEAN 230.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP 50.0 /* amplitude of variation of hue for color scheme C_HUE */
#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -220.0 /* amplitude of variation of hue for color scheme C_HUE */
// #define HUEMEAN 320.0 /* mean value of hue for color scheme C_HUE */
// #define HUEAMP 100.0 /* amplitude of variation of hue for color scheme C_HUE */
/* Basic math */
#define PI 3.141592654
#define DPI 6.283185307
#define PID 1.570796327
double julia_x = -0.5, julia_y = 0.5; /* parameters for Julia sets */
// double julia_x = 0.33267, julia_y = 0.06395; /* parameters for Julia sets */
// double julia_x = 0.37468, julia_y = 0.21115; /* parameters for Julia sets */
#include "global_pdes.c"
#include "sub_wave.c"
double courant2; /* Courant parameter squared */
double courant2, courantb2; /* Courant parameters squared */
void init_wave(x, y, phi, psi, xy_in)
/* initialise field with drop at (x,y) - phi is wave height, psi is phi at time t-1 */
double x, y, *phi[NX], *psi[NX]; short int * xy_in[NX];
{
int i, j;
double xy[2], dist2;
@ -192,11 +167,36 @@ void init_wave(x, y, phi, psi, xy_in)
ij_to_xy(i, j, xy);
dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
phi[i][j] = 0.2*exp(-dist2/0.001)*cos(-sqrt(dist2)/0.01);
if ((xy_in[i][j])||(TWOSPEEDS)) phi[i][j] = 0.2*exp(-dist2/0.001)*cos(-sqrt(dist2)/0.01);
else phi[i][j] = 0.0;
psi[i][j] = 0.0;
}
}
void init_planar_wave(x, y, phi, psi, xy_in)
/* initialise field with drop at (x,y) - phi is wave height, psi is phi at time t-1 */
/* beta version, works for vertical planar wave only so far */
double x, y, *phi[NX], *psi[NX]; short int * xy_in[NX];
{
int i, j;
double xy[2], dist2;
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
dist2 = (xy[0]-x)*(xy[0]-x);
xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
if ((xy_in[i][j])||(TWOSPEEDS)) phi[i][j] = 0.01*exp(-dist2/0.0005)*cos(-sqrt(dist2)/0.01);
else phi[i][j] = 0.0;
psi[i][j] = 0.0;
}
}
void init_wave_flat(phi, psi, xy_in)
/* initialise flat field - phi is wave height, psi is phi at time t-1 */
double *phi[NX], *psi[NX]; short int * xy_in[NX];
@ -232,7 +232,7 @@ double factor, x, y, *phi[NX], *psi[NX];
}
void oscillate_linear_wave(amplitude, t, x1, y1, x2, y2, phi, psi)
/* oscillating boundary condition at (x,y) */
/* oscillating boundary condition at (x,y), beta version */
double amplitude, t, x1, y1, x2, y2, *phi[NX], *psi[NX];
{
int i, j, ij1[2], ij2[2], imin, imax, jmin, jmax, d = 5;
@ -265,6 +265,30 @@ double amplitude, t, x1, y1, x2, y2, *phi[NX], *psi[NX];
/* animation part */
/*********************/
double compute_energy(phi, psi, xy_in, i, j)
double *phi[NX], *psi[NX];
short int *xy_in[NX];
int i, j;
{
double velocity, energy, gradientx2, gradienty2;
int iplus, iminus, jplus, jminus;
velocity = (phi[i][j] - psi[i][j]);
iplus = (i+1); if (iplus == NX) iplus = NX-1;
iminus = (i-1); if (iminus == -1) iminus = 0;
jplus = (j+1); if (jplus == NY) jplus = NY-1;
jminus = (j-1); if (jminus == -1) jminus = 0;
gradientx2 = (phi[iplus][j]-phi[i][j])*(phi[iplus][j]-phi[i][j])
+ (phi[i][j] - phi[iminus][j])*(phi[i][j] - phi[iminus][j]);
gradienty2 = (phi[i][jplus]-phi[i][j])*(phi[i][jplus]-phi[i][j])
+ (phi[i][j] - phi[i][jminus])*(phi[i][j] - phi[i][jminus]);
if (xy_in[i][j]) return(E_SCALE*E_SCALE*(velocity*velocity + 0.5*COURANT*COURANT*(gradientx2+gradienty2)));
else if (TWOSPEEDS) return(E_SCALE*E_SCALE*(velocity*velocity + 0.5*COURANTB*COURANTB*(gradientx2+gradienty2)));
else return(0);
}
void draw_wave(phi, psi, xy_in, scale, time)
/* draw the field */
@ -272,17 +296,28 @@ double *phi[NX], *psi[NX], scale;
short int *xy_in[NX];
int time;
{
int i, j;
double rgb[3], xy[2], x1, y1, x2, y2;
int i, j, iplus, iminus, jplus, jminus;
double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2;
static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
glBegin(GL_QUADS);
// printf("dtinverse = %.5lg\n", dtinverse);
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
if (xy_in[i][j])
if ((TWOSPEEDS)||(xy_in[i][j]))
{
color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb);
if (PLOT == P_AMPLITUDE)
color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb);
else if (PLOT == P_ENERGY)
color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb);
else if (PLOT == P_MIXED)
{
if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb);
else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb);
}
glColor3f(rgb[0], rgb[1], rgb[2]);
glVertex2i(i, j);
@ -303,14 +338,25 @@ void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in)
int i, j, iplus, iminus, jplus, jminus;
double delta, x, y, c, cc;
c = COURANT;
cc = courant2;
// c = COURANT;
// cc = courant2;
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j]){
/* discretized Laplacian */
if (xy_in[i][j])
{
c = COURANT;
cc = courant2;
}
else if (TWOSPEEDS)
{
c = COURANTB;
cc = courantb2;
}
if ((TWOSPEEDS)||(xy_in[i][j])){
/* discretized Laplacian for various boundary conditions */
if ((B_COND == BC_DIRICHLET)||(B_COND == BC_ABSORBING))
{
iplus = (i+1); if (iplus == NX) iplus = NX-1;
@ -327,33 +373,55 @@ void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in)
jminus = (j-1) % NY;
if (jminus < 0) jminus += NY;
}
else if (B_COND == BC_VPER_HABS)
{
iplus = (i+1); if (iplus == NX) iplus = NX-1;
iminus = (i-1); if (iminus == -1) iminus = 0;
jplus = (j+1) % NY;
jminus = (j-1) % NY;
if (jminus < 0) jminus += NY;
}
delta = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j];
x = phi_in[i][j];
y = psi_in[i][j];
/* evolve phi */
if (B_COND != BC_ABSORBING) phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
else
if ((B_COND == BC_PERIODIC)||(B_COND == BC_DIRICHLET))
phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
else if (B_COND == BC_ABSORBING)
{
if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA_BC*(x-y);
/* upper border */
else if (j==NY-1)
phi_out[i][j] = x - c*(x - phi_in[i][NY-2]) - KAPPA*x - GAMMA*(x-y);
phi_out[i][j] = x - c*(x - phi_in[i][NY-2]) - KAPPA_BC*x - GAMMA_BC*(x-y);
/* lower border */
else if (j==0)
phi_out[i][j] = x - c*(x - phi_in[i][1]) - KAPPA*x - GAMMA*(x-y);
phi_out[i][j] = x - c*(x - phi_in[i][1]) - KAPPA_BC*x - GAMMA_BC*(x-y);
/* right border */
if (i==NX-1)
phi_out[i][j] = x - c*(x - phi_in[NX-2][j]) - KAPPA*x - GAMMA*(x-y);
phi_out[i][j] = x - c*(x - phi_in[NX-2][j]) - KAPPA_BC*x - GAMMA_BC*(x-y);
/* left border */
else if (i==0)
phi_out[i][j] = x - c*(x - phi_in[1][j]) - KAPPA*x - GAMMA*(x-y);
phi_out[i][j] = x - c*(x - phi_in[1][j]) - KAPPA_BC*x - GAMMA_BC*(x-y);
}
else if (B_COND == BC_VPER_HABS)
{
if ((i>0)&&(i<NX-1))
phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
/* right border */
else if (i==NX-1)
phi_out[i][j] = x - c*(x - phi_in[NX-2][j]) - KAPPA_BC*x - GAMMA_BC*(x-y);
/* left border */
else if (i==0)
phi_out[i][j] = x - c*(x - phi_in[1][j]) - KAPPA_BC*x - GAMMA_BC*(x-y);
}
psi_out[i][j] = x;
@ -380,51 +448,7 @@ void evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in)
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
}
void old_evolve_wave(phi, psi, xy_in)
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
double *phi[NX], *psi[NX]; short int *xy_in[NX];
{
int i, j, iplus, iminus, jplus, jminus;
double delta, x, y;
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j]){
/* discretized Laplacian */
iplus = (i+1) % NX;
iminus = (i-1) % NX;
if (iminus < 0) iminus += NX;
jplus = (j+1) % NY;
jminus = (j-1) % NY;
if (jminus < 0) jminus += NY;
delta = phi[iplus][j] + phi[iminus][j] + phi[i][jplus] + phi[i][jminus] - 4.0*phi[i][j];
x = phi[i][j];
y = psi[i][j];
/* evolve phi */
phi[i][j] = -y + 2*x + courant2*delta - KAPPA*x - GAMMA*(x-y);
/* Old versions of the simulation used this: */
// phi[i][j] = (-psi[i][j] + 2*phi[i][j] + courant2*delta)*damping;
// where damping = 1.0 - 0.0001;
psi[i][j] = x;
if (FLOOR)
{
if (phi[i][j] > VMAX) phi[i][j] = VMAX;
if (phi[i][j] < -VMAX) phi[i][j] = -VMAX;
if (psi[i][j] > VMAX) psi[i][j] = VMAX;
if (psi[i][j] < -VMAX) psi[i][j] = -VMAX;
}
}
}
}
// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
}
double compute_variance(phi, psi, xy_in)
@ -448,6 +472,7 @@ double compute_variance(phi, psi, xy_in)
}
void animation()
{
double time, scale;
@ -464,12 +489,18 @@ void animation()
psi_tmp[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
/* initialise positions and radii of circles */
if (B_DOMAIN == D_CIRCLES) init_circle_config();
courant2 = COURANT*COURANT;
courantb2 = COURANTB*COURANTB;
/* initialize wave with a drop at one point, zero elsewhere */
// init_wave_flat(phi, psi, xy_in);
init_wave(0.0, 0.0, phi, psi, xy_in);
init_planar_wave(XMIN + 0.01, 0.0, phi, psi, xy_in);
// init_planar_wave(XMIN + 1.0, 0.0, phi, psi, xy_in);
// init_wave(-1.5, 0.0, phi, psi, xy_in);
// init_wave(0.0, 0.0, phi, psi, xy_in);
/* add a drop at another point */
// add_drop_to_wave(1.0, 0.7, 0.0, phi, psi);
@ -487,7 +518,7 @@ void animation()
sleep(SLEEP1);
for (i=0; i<=NSTEPS; i++)
for (i=0; i<=INITIAL_TIME + NSTEPS; i++)
{
//printf("%d\n",i);
/* compute the variance of the field to adjust color scheme */
@ -499,16 +530,12 @@ void animation()
}
else scale = 1.0;
// /* TO BE ADAPTED */
// scale = 1.0;
draw_wave(phi, psi, xy_in, scale, i);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
// if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
}
// for (j=0; j<NVID; j++) evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
draw_billiard();
@ -517,7 +544,8 @@ void animation()
if (MOVIE)
{
save_frame();
if (i >= INITIAL_TIME) save_frame();
else printf("Initial phase time %i of %i\n", i, INITIAL_TIME);
/* it seems that saving too many files too fast can cause trouble with the file system */
/* so this is to make a pause from time to time - parameter PAUSE may need adjusting */