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-05-29 16:50:49 +02:00 committed by GitHub
parent c5137a51c5
commit bd6aa073a7
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
6 changed files with 2319 additions and 855 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

191
drop_billiard.c
View File

@ -33,14 +33,18 @@
#define WINWIDTH 1280 /* window width */ #define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */ #define WINHEIGHT 720 /* window height */
#define XMIN -2.0 #define XMIN -1.8
#define XMAX 2.0 /* x interval */ #define XMAX 1.8 /* x interval */
#define YMIN -1.125 #define YMIN -0.91
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #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 */
/* Choice of the billiard table */ /* Choice of the billiard table */
#define B_DOMAIN 5 /* choice of domain shape */ #define B_DOMAIN 9 /* choice of domain shape */
#define D_RECTANGLE 0 /* rectangular domain */ #define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */ #define D_ELLIPSE 1 /* elliptical domain */
@ -48,36 +52,59 @@
#define D_SINAI 3 /* Sinai billiard */ #define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */ #define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular 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 LAMBDA 1.5 /* parameter controlling shape of billiard */ // #define LAMBDA 1.0 /* parameter controlling shape of billiard */
#define LAMBDA 1.73205080756888 /* sqrt(3) for triangle tiling plane */ #define LAMBDA 1.124950941 /* sin(36°)/sin(31.5°) for 5-star shape with 45° angles */
// #define LAMBDA 1.445124904 /* sin(36°)/sin(24°) for 5-star shape with 60° angles */
// #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 FOCI 1 /* set to 1 to draw focal points of ellipse */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NPOLY 5 /* number of sides of polygon */
#define APOLY -1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */ #define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */
#define NPART 10000 /* number of particles */ #define NPART 100000 /* number of particles */
#define NPARTMAX 20000 /* maximal number of particles after resampling */ #define NPARTMAX 100000 /* maximal number of particles after resampling */
#define LMAX 0.01 /* minimal segment length triggering resampling */
#define LPERIODIC 2.0 /* lines longer than this are not drawn (useful for Sinai billiard) */
#define DMIN 0.02 /* minimal distance to boundary for triggering resampling */
#define MARGIN 1.0 /* distance above which points of curve are not drawn */
#define CYCLE 0 /* set to 1 for closed curve (start in all directions) */
#define NSTEPS 10000 /* number of frames of movie */ #define NSTEPS 4000 /* number of frames of movie */
#define TIME 100 /* time between movie frames, for fluidity of real-time simulation */ #define TIME 15 /* time between movie frames, for fluidity of real-time simulation */
#define DPHI 0.0001 /* integration step */ #define DPHI 0.0001 /* integration step */
#define NVID 50 /* number of iterations between images displayed on screen */ #define NVID 10 /* number of iterations between images displayed on screen */
#define NCOLORS 10 /* number of colors */
#define COLORSHIFT 200 /* hue of initial color */
#define NSEG 100 /* number of segments of boundary */
#define BLACK 1 /* set to 1 for black background */
/* Decreasing TIME accelerates the animation and the movie */ /* Decreasing TIME accelerates the animation and the movie */
/* For constant speed of movie, TIME*DPHI should be kept constant */ /* For constant speed of movie, TIME*DPHI should be kept constant */
/* However, increasing DPHI too much deterioriates quality of simulation */ /* However, increasing DPHI too much deterioriates quality of simulation */
/* For a good quality movie, take for instance TIME = 50, DPHI = 0.0002 */ /* For a good quality movie, take for instance TIME = 50, DPHI = 0.0002 */
/* 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 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 COLORSHIFT 200 /* hue of initial color */
#define NSEG 100 /* number of segments of boundary */
#define BILLIARD_WIDTH 4 /* width of billiard */
#define FRONT_WIDTH 3 /* width of wave front */
#define BLACK 1 /* set to 1 for black background */
#define COLOR_OUTSIDE 0 /* set to 1 for colored outside */
#define OUTER_COLOR 300.0 /* color outside billiard */
#define PAINT_INT 1 /* set to 1 to paint interior in other color (for polygon) */
#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 PSLEEP 1 /* sleep time during pause */
#define SLEEP1 1 /* initial sleeping time */ #define SLEEP1 1 /* initial sleeping time */
@ -89,6 +116,7 @@
#include "sub_part_billiard.c" #include "sub_part_billiard.c"
/*********************/ /*********************/
/* animation part */ /* animation part */
/*********************/ /*********************/
@ -126,7 +154,8 @@ double *configs[NPARTMAX];
double dalpha, alpha, pos[2]; double dalpha, alpha, pos[2];
while (angle2 < angle1) angle2 += DPI; while (angle2 < angle1) angle2 += DPI;
dalpha = (angle2 - angle1)/((double)(NPART-1)); dalpha = (angle2 - angle1)/((double)(NPART));
// dalpha = (angle2 - angle1)/((double)(NPART-1));
for (i=0; i<NPART; i++) for (i=0; i<NPART; i++)
{ {
alpha = angle1 + dalpha*((double)i); alpha = angle1 + dalpha*((double)i);
@ -227,19 +256,19 @@ double *configs[NPARTMAX];
else return(1); else return(1);
} }
void draw_config(color, configs) void draw_config(color, configs)
/* draw the wave front */ /* draw the wave front by ordering colors */
int color[NPARTMAX]; int color[NPARTMAX];
double *configs[NPARTMAX]; double *configs[NPARTMAX];
{ {
int i; int i;
double x1, y1, x2, y2, cosphi, sinphi, rgb[3], dist; double x1, y1, x2, y2, cosphi, sinphi, rgb[3], dist, dmax;
glutSwapBuffers(); glutSwapBuffers();
blank(); blank();
if (PAINT_INT) paint_billiard_interior();
glLineWidth(5.0); glLineWidth(FRONT_WIDTH);
glEnable(GL_LINE_SMOOTH); glEnable(GL_LINE_SMOOTH);
if (CYCLE) glBegin (GL_LINE_LOOP); if (CYCLE) glBegin (GL_LINE_LOOP);
@ -247,17 +276,6 @@ double *configs[NPARTMAX];
for (i=0; i<nparticles; i++) for (i=0; i<nparticles; i++)
{ {
// print_config(configs[i]);
if (configs[i][2]<0.0)
{
vbilliard(configs[i]);
color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS;
}
configs[i][2] += DPHI;
cosphi = (configs[i][6] - configs[i][4])/configs[i][3]; cosphi = (configs[i][6] - configs[i][4])/configs[i][3];
sinphi = (configs[i][7] - configs[i][5])/configs[i][3]; sinphi = (configs[i][7] - configs[i][5])/configs[i][3];
x2 = configs[i][4] + configs[i][2]*cosphi; x2 = configs[i][4] + configs[i][2]*cosphi;
@ -267,18 +285,20 @@ double *configs[NPARTMAX];
if (i>0) dist = module2(x2-x1,y2-y1); if (i>0) dist = module2(x2-x1,y2-y1);
else dist = 0.0; else dist = 0.0;
dmax = DPI*((double)global_time)*DPHI/((double)nparticles);
/* expected maximal distance between points for growing circle */
rgb_color_scheme(color[i], rgb); rgb_color_scheme(color[i], rgb);
glColor3d(rgb[0], rgb[1], rgb[2]); glColor3d(rgb[0], rgb[1], rgb[2]);
/* draw line only if it does not exceed */ /* draw line only if it does not exceed LPERIODIC and 2*dmax */
if ((xy_in_billiard(x2, y2))&&(dist < LPERIODIC)) glVertex2d(x2, y2); if ((xy_in_billiard(x2, y2))&&(dist < LPERIODIC)&&(dist < LCUT*dmax)) glVertex2d(x2, y2);
else else
{ {
glEnd(); glEnd();
glBegin (GL_LINE_STRIP); glBegin (GL_LINE_STRIP);
} }
if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3]; if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3];
/* keep track of previous point to determine segment length */ /* keep track of previous point to determine segment length */
@ -290,6 +310,67 @@ double *configs[NPARTMAX];
} }
void draw_ordered_config(color, configs)
/* draw the wave front, one color after the other */
int color[NPARTMAX];
double *configs[NPARTMAX];
{
int i, col;
double x1, y1, x2, y2, cosphi, sinphi, rgb[3], dist, dmax;
glutSwapBuffers();
blank();
if (PAINT_INT) paint_billiard_interior();
glLineWidth(FRONT_WIDTH);
glEnable(GL_LINE_SMOOTH);
for (col=0; col<NCOLORS; col++)
{
glBegin(GL_LINE_STRIP);
for (i=0; i<nparticles; i++)
{
if (color[i] == col)
{
cosphi = (configs[i][6] - configs[i][4])/configs[i][3];
sinphi = (configs[i][7] - configs[i][5])/configs[i][3];
x2 = configs[i][4] + configs[i][2]*cosphi;
y2 = configs[i][5] + configs[i][2]*sinphi;
/* determine length of segment to avoid drawing too long segments */
if (i>0) dist = module2(x2-x1,y2-y1);
else dist = 0.0;
dmax = DPI*((double)global_time)*DPHI/((double)nparticles);
/* expected maximal distance between points for growing circle */
rgb_color_scheme(color[i], rgb);
glColor3d(rgb[0], rgb[1], rgb[2]);
/* draw line only if it does not exceed LPERIODIC and 2*dmax */
if ((i>0)&&(xy_in_billiard(x2, y2))&&(dist < LPERIODIC)&&(dist < LCUT*dmax)) glVertex2d(x2, y2);
else
{
glEnd();
glBegin (GL_LINE_STRIP);
}
if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3];
/* keep track of previous point to determine segment length */
x1 = x2;
y1 = y2;
}
}
glEnd ();
}
draw_billiard(LAMBDA);
}
void graph_movie(time, color, configs) void graph_movie(time, color, configs)
/* compute next movie frame */ /* compute next movie frame */
int time, color[NPARTMAX]; int time, color[NPARTMAX];
@ -299,6 +380,7 @@ double *configs[NPARTMAX];
for (j=0; j<time; j++) for (j=0; j<time; j++)
{ {
global_time++;
for (i=0; i<nparticles; i++) for (i=0; i<nparticles; i++)
{ {
// print_config(configs[i]); // print_config(configs[i]);
@ -315,9 +397,6 @@ double *configs[NPARTMAX];
if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3]; if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3];
} }
} }
draw_config(color, configs);
} }
void graph_no_movie(time, color, configs) void graph_no_movie(time, color, configs)
@ -329,10 +408,9 @@ double *configs[NPARTMAX];
for (j=0; j<time; j++) for (j=0; j<time; j++)
{ {
global_time++;
for (i=0; i<nparticles; i++) for (i=0; i<nparticles; i++)
{ {
// print_config(configs[i]);
if (configs[i][2]<0.0) if (configs[i][2]<0.0)
{ {
vbilliard(configs[i]); vbilliard(configs[i]);
@ -347,9 +425,6 @@ double *configs[NPARTMAX];
if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3]; if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3];
} }
} }
draw_config(color, configs);
draw_billiard(color, configs);
} }
@ -368,20 +443,15 @@ void animation()
for (i=0; i<NPARTMAX; i++) for (i=0; i<NPARTMAX; i++)
configs[i] = (double *)malloc(8*sizeof(double)); configs[i] = (double *)malloc(8*sizeof(double));
init_drop_config(0.0, -0.5, 0.0, DPI, configs); init_drop_config(0.0, 0.1, 0.0, DPI, configs);
// init_drop_config(0.0, -1.0, 0.0, PI, configs);
// init_drop_config(0.95, 0.95, PI, 3.0*DPI, configs);
// init_drop_config(1.4, 0.9, DPI, configs);
// init_boundary_config(1.5, 1.5, 0.0, PI, configs); // init_boundary_config(1.5, 1.5, 0.0, PI, configs);
// other possible initial conditions :
// init_drop_config(sqrt(LAMBDA*LAMBDA-1.0) - 0.1,0.0, 0.0, DPI, configs); /* Start at focus */
// init_boundary_config(0.0, 0.0, 0.0, PI, configs);
// init_drop_config(LAMBDA-0.01, 0.0, PID, 3.0*PID, configs);
blank(); blank();
glColor3d(0.0, 0.0, 0.0); glColor3d(0.0, 0.0, 0.0);
draw_billiard(LAMBDA); draw_billiard(LAMBDA);
if (PAINT_INT) paint_billiard_interior();
glutSwapBuffers(); glutSwapBuffers();
@ -395,6 +465,11 @@ void animation()
if (MOVIE) graph_movie(TIME, color, configs); if (MOVIE) graph_movie(TIME, color, configs);
else graph_no_movie(NVID, color, configs); else graph_no_movie(NVID, color, configs);
if (ORDER_COLORS) draw_ordered_config(color, configs);
else draw_config(color, configs);
draw_billiard();
/* for the ellipse, paths passing close to the foci are stronly divergent /* for the ellipse, paths passing close to the foci are stronly divergent
* and the configurations may need to be resampled be adding extra points */ * and the configurations may need to be resampled be adding extra points */
if ((RESAMPLE)&&(i % 5 == 0)&&(nparticles < NPARTMAX)) resamp = resample(color, configs); if ((RESAMPLE)&&(i % 5 == 0)&&(nparticles < NPARTMAX)) resamp = resample(color, configs);

250
particle_billiard.c
View File

@ -39,10 +39,14 @@
#define XMAX 2.0 /* x interval */ #define XMAX 2.0 /* x interval */
#define YMIN -1.125 #define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #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 */
/* Choice of the billiard table */ /* Choice of the billiard table */
#define B_DOMAIN 3 /* choice of domain shape */ #define B_DOMAIN 9 /* choice of domain shape */
#define D_RECTANGLE 0 /* rectangular domain */ #define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */ #define D_ELLIPSE 1 /* elliptical domain */
@ -50,30 +54,35 @@
#define D_SINAI 3 /* Sinai billiard */ #define D_SINAI 3 /* Sinai billiard */
#define D_DIAMOND 4 /* diamond-shaped billiard */ #define D_DIAMOND 4 /* diamond-shaped billiard */
#define D_TRIANGLE 5 /* triangular 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 LAMBDA 0.5 /* parameter controlling shape of billiard */ // #define LAMBDA 5.0 /* parameter controlling shape of billiard */
#define LAMBDA -1.949855824 /* 7-sided Reuleaux triangle */
// #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 LAMBDA 1.73205080756888 /* sqrt(3) for triangle tiling plane */
#define MU 0.1 /* second parameter controlling shape of billiard */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */
#define NPOLY 7 /* number of sides of polygon */
#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */
#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */ #define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */
#define DEBUG 0 /* draw trajectories, for debugging purposes */
#define NPART 5000 /* number of particles */ /* Simulation parameters */
#define NPARTMAX 50000 /* maximal number of particles after resampling */
#define NPART 10000 /* number of particles */
#define NPARTMAX 100000 /* maximal number of particles after resampling */
#define LMAX 0.01 /* minimal segment length triggering resampling */ #define LMAX 0.01 /* minimal segment length triggering resampling */
#define DMIN 0.02 /* minimal distance to boundary for 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 CYCLE 1 /* set to 1 for closed curve (start in all directions) */
#define NSTEPS 5000 /* number of frames of movie */ #define NSTEPS 4000 /* number of frames of movie */
#define TIME 1500 /* time between movie frames, for fluidity of real-time simulation */ #define TIME 125 /* time between movie frames, for fluidity of real-time simulation */
#define DPHI 0.00001 /* integration step */ #define DPHI 0.00001 /* integration step */
#define NVID 750 /* number of iterations between images displayed on screen */ #define NVID 150 /* number of iterations between images displayed on screen */
#define NCOLORS 10 /* number of colors */
#define COLORSHIFT 0 /* hue of initial color */
#define NSEG 100 /* number of segments of boundary */
#define LENGTH 0.05 /* length of velocity vectors */
#define BLACK 1 /* set to 1 for black background */
/* Decreasing TIME accelerates the animation and the movie */ /* Decreasing TIME accelerates the animation and the movie */
/* For constant speed of movie, TIME*DPHI should be kept constant */ /* For constant speed of movie, TIME*DPHI should be kept constant */
@ -81,6 +90,22 @@
/* NVID tells how often a picture is drawn in the animation, increase it for faster anim */ /* NVID tells how often a picture is drawn in the animation, increase it for faster anim */
/* For a good quality movie, take for instance TIME = 400, DPHI = 0.00005, NVID = 100 */ /* For a good quality movie, take for instance TIME = 400, DPHI = 0.00005, NVID = 100 */
/* Colors and other graphical parameters */
#define NCOLORS -7 /* number of colors */
#define COLORSHIFT 220 /* hue of initial color */
#define NSEG 100 /* number of segments of boundary */
#define LENGTH 0.01 /* length of velocity vectors */
#define BILLIARD_WIDTH 4 /* width of billiard */
#define PARTICLE_WIDTH 2 /* width of particles */
#define FRONT_WIDTH 3 /* width of wave front */
#define BLACK 1 /* set to 1 for black background */
#define COLOR_OUTSIDE 0 /* set to 1 for colored outside */
#define OUTER_COLOR 270.0 /* color outside billiard */
#define PAINT_INT 1 /* set to 1 to paint interior in other color (for polygon) */
#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 PSLEEP 1 /* sleep time during pause */
#define SLEEP1 1 /* initial sleeping time */ #define SLEEP1 1 /* initial sleeping time */
@ -92,6 +117,7 @@
#include "sub_part_billiard.c" #include "sub_part_billiard.c"
/*********************/ /*********************/
/* animation part */ /* animation part */
/*********************/ /*********************/
@ -119,8 +145,6 @@ double *configs[NPARTMAX];
alpha = theta + angle; alpha = theta + angle;
vbilliard_xy(configs[i], alpha, pos); vbilliard_xy(configs[i], alpha, pos);
/*
vbilliard(configs[i], alpha, pos);*/
} }
} }
@ -146,6 +170,35 @@ double *configs[NPARTMAX];
} }
} }
void init_sym_drop_config(x0, y0, angle1, angle2, configs)
/* initialize configuration with two symmetric partial drops */
double x0, y0, angle1, angle2;
double *configs[NPARTMAX];
{
int i;
double dalpha, alpha, meanangle;
double conf[2], pos[2];
while (angle2 < angle1) angle2 += DPI;
meanangle = 0.5*(angle1 + angle2);
dalpha = (angle2 - angle1)/((double)(NPART-1));
for (i=0; i<NPART/2; i++)
{
alpha = meanangle + dalpha*((double)i);
pos[0] = x0;
pos[1] = y0;
vbilliard_xy(configs[i], alpha, pos);
}
for (i=0; i<NPART/2; i++)
{
alpha = meanangle - dalpha*((double)i);
pos[0] = x0;
pos[1] = y0;
vbilliard_xy(configs[NPART/2 + i], alpha, pos);
}
}
void init_line_config(x0, y0, x1, y1, angle, configs) /* initialize configuration: line (x0,y0)-(x1,y1) in direction alpha */ void init_line_config(x0, y0, x1, y1, angle, configs) /* initialize configuration: line (x0,y0)-(x1,y1) in direction alpha */
double x0, y0, x1, y1, angle; double x0, y0, x1, y1, angle;
double *configs[NPARTMAX]; double *configs[NPARTMAX];
@ -154,8 +207,10 @@ double *configs[NPARTMAX];
double dx, dy; double dx, dy;
double conf[2], pos[2]; double conf[2], pos[2];
dx = (x1-x0)/((double)(NPART-1)); dx = (x1-x0)/((double)(NPART));
dy = (y1-y0)/((double)(NPART-1)); dy = (y1-y0)/((double)(NPART));
// dx = (x1-x0)/((double)(NPART-1));
// dy = (y1-y0)/((double)(NPART-1));
for (i=0; i<NPART; i++) for (i=0; i<NPART; i++)
{ {
pos[0] = x0 + ((double)i)*dx; pos[0] = x0 + ((double)i)*dx;
@ -165,9 +220,9 @@ double *configs[NPARTMAX];
} }
void draw_config(color, configs) void draw_config(color, configs, active)
/* draw the particles */ /* draw the particles */
int color[NPARTMAX]; int color[NPARTMAX], active[NPARTMAX];
double *configs[NPARTMAX]; double *configs[NPARTMAX];
{ {
int i; int i;
@ -176,7 +231,7 @@ double *configs[NPARTMAX];
glutSwapBuffers(); glutSwapBuffers();
blank(); blank();
glLineWidth(3.0); glLineWidth(PARTICLE_WIDTH);
glEnable(GL_LINE_SMOOTH); glEnable(GL_LINE_SMOOTH);
@ -198,65 +253,73 @@ double *configs[NPARTMAX];
x2 = configs[i][4] + (configs[i][2] + LENGTH)*cosphi; x2 = configs[i][4] + (configs[i][2] + LENGTH)*cosphi;
y2 = configs[i][5] + (configs[i][2] + LENGTH)*sinphi; y2 = configs[i][5] + (configs[i][2] + LENGTH)*sinphi;
rgb_color_scheme(color[i], rgb); /* test whether particle does not escape billiard */
glColor3f(rgb[0], rgb[1], rgb[2]); if (active[i]) active[i] = xy_in_billiard(x1, y1);
glBegin(GL_LINE_STRIP); if (active[i])
glVertex2d(x1, y1);
glVertex2d(x2, y2);
glEnd ();
/* taking care of boundary conditions - only needed for periodic boundary conditions */
if (x2 > XMAX)
{ {
rgb_color_scheme(color[i], rgb);
glColor3f(rgb[0], rgb[1], rgb[2]);
glBegin(GL_LINE_STRIP); glBegin(GL_LINE_STRIP);
glVertex2d(x1+XMIN-XMAX, y1); glVertex2d(x1, y1);
glVertex2d(x2+XMIN-XMAX, y2); glVertex2d(x2, y2);
glEnd (); glEnd ();
/* taking care of boundary conditions - only needed for periodic boundary conditions */
if (x2 > XMAX)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1+XMIN-XMAX, y1);
glVertex2d(x2+XMIN-XMAX, y2);
glEnd ();
}
if (x2 < XMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1-XMIN+XMAX, y1);
glVertex2d(x2-XMIN+XMAX, y2);
glEnd ();
}
if (y2 > YMAX)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1, y1+YMIN-YMAX);
glVertex2d(x2, y2+YMIN-YMAX);
glEnd ();
}
if (y2 < YMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1, y1+YMAX-YMIN);
glVertex2d(x2, y2+YMAX-YMIN);
glEnd ();
}
} }
if (x2 < XMIN) /* draw trajectories, for debugging purpose */
if (DEBUG)
{ {
glBegin(GL_LINE_STRIP); glLineWidth(1.0);
glVertex2d(x1-XMIN+XMAX, y1); glBegin(GL_LINES);
glVertex2d(x2-XMIN+XMAX, y2); glVertex2d(configs[i][4], configs[i][5]);
glVertex2d(configs[i][6], configs[i][7]);
glEnd (); glEnd ();
glLineWidth(3.0);
} }
if (y2 > YMAX)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1, y1+YMIN-YMAX);
glVertex2d(x2, y2+YMIN-YMAX);
glEnd ();
}
if (y2 < YMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1, y1+YMAX-YMIN);
glVertex2d(x2, y2+YMAX-YMIN);
glEnd ();
}
/* for debugging purpose */
// glLineWidth(1.0);
// glBegin(GL_LINES);
// glVertex2d(configs[i][4], configs[i][5]);
// glVertex2d(configs[i][6], configs[i][7]);
// glEnd ();
// glLineWidth(3.0);
if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3]; if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3];
} }
draw_billiard(LAMBDA); draw_billiard();
} }
void graph_movie(time, color, configs) void graph_movie(time, color, configs, active)
/* compute next movie frame */ /* compute next movie frame */
int time, color[NPARTMAX]; int time, color[NPARTMAX], active[NPARTMAX];
double *configs[NPARTMAX]; double *configs[NPARTMAX];
{ {
int i, j, c; int i, j, c;
@ -265,13 +328,13 @@ double *configs[NPARTMAX];
{ {
for (i=0; i<nparticles; i++) for (i=0; i<nparticles; i++)
{ {
// print_config(configs[i]);
if (configs[i][2]<0.0) if (configs[i][2]<0.0)
{ {
c = vbilliard(configs[i]); c = vbilliard(configs[i]);
if (c>=0) color[i]++; // if (c>=0) color[i]++;
color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS; if (color[i] >= NCOLORS) color[i] -= NCOLORS;
} }
configs[i][2] += DPHI; configs[i][2] += DPHI;
@ -286,60 +349,38 @@ double *configs[NPARTMAX];
// draw_config(color, configs); // draw_config(color, configs);
} }
void graph_no_movie(time, color, configs)
/* plot next image without making a movie */
int time, color[NPARTMAX];
double *configs[NPARTMAX];
{
int i, j, c;
for (j=0; j<time; j++)
{
for (i=0; i<nparticles; i++)
{
// print_config(configs[i]);
if (configs[i][2]<0.0)
{
c = vbilliard(configs[i]);
if (c>=0) color[i]++;
if (color[i] >= NCOLORS) color[i] -= NCOLORS;
}
configs[i][2] += DPHI;
if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3];
}
}
}
void animation() void animation()
{ {
double time, dt; double time, dt, alpha;
double *configs[NPARTMAX]; double *configs[NPARTMAX];
int i, j, resamp = 1, s; int i, j, resamp = 1, s;
int *color, *newcolor; int *color, *newcolor, *active;
/* Since NPARTMAX can be big, it seemed wiser to use some memory allocation here */ /* Since NPARTMAX can be big, it seemed wiser to use some memory allocation here */
color = malloc(sizeof(int)*(NPARTMAX)); color = malloc(sizeof(int)*(NPARTMAX));
newcolor = malloc(sizeof(int)*(NPARTMAX)); newcolor = malloc(sizeof(int)*(NPARTMAX));
active = malloc(sizeof(int)*(NPARTMAX));
for (i=0; i<NPARTMAX; i++) for (i=0; i<NPARTMAX; i++)
configs[i] = (double *)malloc(8*sizeof(double)); configs[i] = (double *)malloc(8*sizeof(double));
/* initialize system by putting particles in a given point with a range of velocities */ /* initialize system by putting particles in a given point with a range of velocities */
init_drop_config(-1.5, 0.3, -0.1, 0.1, configs); alpha = 3.0*PI/7.0;
init_sym_drop_config(-0.99, 0.0, -alpha, alpha, configs);
// init_drop_config(-0.999, 0.0, -alpha, alpha, configs);
// init_drop_config(0.0, 0.5, 0.0, DPI, configs);
// other possible initial conditions : // other possible initial conditions :
// init_line_config(0.0, -0.05, 0.0, 0.05, 0.0, configs); // init_line_config(-0.6, 0.2, -0.6, 0.7, 0.0, configs);
// init_line_config(-0.7, -0.45, -0.7, 0.45, 0.0, configs);
// init_line_config(0.0, 0.1, 0.0, 0.7, 0.0, configs);
blank(); blank();
glColor3f(0.0, 0.0, 0.0); glColor3f(0.0, 0.0, 0.0);
draw_billiard(LAMBDA); draw_billiard();
glutSwapBuffers(); glutSwapBuffers();
@ -348,17 +389,17 @@ void animation()
{ {
color[i] = 0; color[i] = 0;
newcolor[i] = 0; newcolor[i] = 0;
active[i] = 1;
} }
sleep(SLEEP1); sleep(SLEEP1);
for (i=0; i<=NSTEPS; i++) for (i=0; i<=NSTEPS; i++)
{ {
if (MOVIE) graph_movie(TIME, newcolor, configs); graph_movie(TIME, newcolor, configs, active);
else graph_no_movie(NVID, newcolor, configs);
draw_config(color, configs); draw_config(newcolor, configs, active);
draw_billiard(color, configs); draw_billiard();
for (j=0; j<NPARTMAX; j++) color[j] = newcolor[j]; for (j=0; j<NPARTMAX; j++) color[j] = newcolor[j];
@ -384,6 +425,7 @@ void animation()
} }
free(color); free(color);
free(newcolor);
for (i=0; i<NPARTMAX; i++) free(configs[i]); for (i=0; i<NPARTMAX; i++) free(configs[i]);
} }

513
schrodinger.c Normal file
View File

@ -0,0 +1,513 @@
/*********************************************************************************/
/* */
/* Animation of Schrödinger equation in a planar domain */
/* */
/* N. Berglund, May 2021 */
/* */
/* Feel free to reuse, but if doing so it would be nice to drop a */
/* line to nils.berglund@univ-orleans.fr - Thanks! */
/* */
/* compile with */
/* gcc -o schrodinger schrodinger.c */
/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */
/* */
/* To make a video, set MOVIE to 1 and create subfolder tif_schrod */
/* It may be possible to increase parameter PAUSE */
/* */
/* create movie using */
/* ffmpeg -i wave.%05d.tif -vcodec libx264 wave.mp4 */
/* */
/*********************************************************************************/
/*********************************************************************************/
/* */
/* NB: The algorithm used to simulate the wave equation is highly paralellizable */
/* One could make it much faster by using a GPU */
/* */
/*********************************************************************************/
#include <math.h>
#include <string.h>
#include <GL/glut.h>
#include <GL/glu.h>
#include <unistd.h>
#include <sys/types.h>
#include <tiffio.h> /* Sam Leffler's libtiff library. */
#include <omp.h>
#define MOVIE 0 /* set to 1 to generate movie */
/* General geometrical parameters */
#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 */
/* 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 YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
/* Choice of the billiard table */
#define B_DOMAIN 3 /* 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 LAMBDA 0.2 /* 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 FOCI 1 /* set to 1 to draw focal points of ellipse */
/* 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.00000002
// #define DT 0.000000005
#define HBAR 1.0
/* Boundary conditions */
#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 750 /* 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 */
#define PSLEEP 1 /* sleep time during pause */
#define SLEEP1 1 /* initial sleeping time */
#define SLEEP2 1 /* final sleeping time */
/* 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 type */
#define PLOT 0
#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 */
#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 */
#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 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
#include "sub_wave.c"
double courant2; /* Courant parameter squared */
double dx2; /* spatial step size squared */
double intstep; /* integration step */
double intstep1; /* integration step used in absorbing boundary conditions */
void init_coherent_state(x, y, px, py, scalex, phi, psi, xy_in)
/* initialise field with coherent state of position (x,y) and momentum (px, py) */
/* phi is real part, psi is imaginary part */
double x, y, px, py, scalex, *phi[NX], *psi[NX];
short int * xy_in[NX];
{
int i, j;
double xy[2], dist2, module, phase, scale2;
scale2 = scalex*scalex;
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
if (xy_in[i][j])
{
dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
module = exp(-dist2/scale2);
if (module < 1.0e-15) module = 1.0e-15;
phase = (px*(xy[0]-x) + py*(xy[1]-y))/scalex;
phi[i][j] = module*cos(phase);
psi[i][j] = module*sin(phase);
}
else
{
phi[i][j] = 0.0;
psi[i][j] = 0.0;
}
}
}
/*********************/
/* animation part */
/*********************/
void schrodinger_color_scheme(phi, psi, scale, time, rgb)
double phi, psi, scale, rgb[3];
int time;
{
double phase, amp, lum;
if (PLOT == P_MODULE)
color_scheme(COLOR_SCHEME, 2.0*module2(phi, psi)-1.0, scale, time, rgb);
else if (PLOT == P_PHASE)
{
amp = module2(phi,psi);
// if (amp < 1.0e-10) amp = 1.0e-10;
phase = argument(phi/amp, psi/amp);
if (phase < 0.0) phase += DPI;
lum = (color_amplitude(amp, scale, time))*0.5;
if (lum < 0.0) lum = 0.0;
hsl_to_rgb(phase*360.0/DPI, 0.9, lum, rgb);
}
else if (PLOT == P_REAL) color_scheme(COLOR_SCHEME, phi, scale, time, rgb);
else if (PLOT == P_IMAGINARY) color_scheme(COLOR_SCHEME, psi, scale, time, rgb);
}
void draw_wave(phi, psi, xy_in, scale, time)
/* draw the field */
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, amp, phase;
glBegin(GL_QUADS);
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
if (xy_in[i][j])
{
schrodinger_color_scheme(phi[i][j],psi[i][j], scale, time, rgb);
glColor3f(rgb[0], rgb[1], rgb[2]);
glVertex2i(i, j);
glVertex2i(i+1, j);
glVertex2i(i+1, j+1);
glVertex2i(i, j+1);
}
}
glEnd ();
}
void evolve_wave(phi, psi, 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];
{
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]){
/* 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[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];
x = phi[i][j];
y = psi[i][j];
/* evolve phi and psi */
if (B_COND != BC_ABSORBING)
{
phi[i][j] = x - intstep*delta2;
psi[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))
{
phi[i][j] = x - intstep*delta2;
psi[i][j] = y + intstep*delta1;
}
/* right border */
else if (i==NX-1)
{
phi[i][j] = x - intstep1*(y - psi[i-1][j]);
psi[i][j] = y + intstep1*(x - phi[i-1][j]);
}
/* upper border */
else if (j==NY-1)
{
phi[i][j] = x - intstep1*(y - psi[i][j-1]);
psi[i][j] = y + intstep1*(x - phi[i][j-1]);
}
/* left border */
else if (i==0)
{
phi[i][j] = x - intstep1*(y - psi[1][j]);
psi[i][j] = y + intstep1*(x - phi[1][j]);
}
/* lower border */
else if (j==0)
{
phi[i][j] = x - intstep1*(y - psi[i][1]);
psi[i][j] = y + intstep1*(x - phi[i][1]);
}
}
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)
/* compute the variance (total probability) of the field */
double *phi[NX], *psi[NX]; short int * xy_in[NX];
{
int i, j, n = 0;
double variance = 0.0;
for (i=1; i<NX; i++)
for (j=1; j<NY; j++)
{
if (xy_in[i][j])
{
n++;
variance += phi[i][j]*phi[i][j] + psi[i][j]*psi[i][j];
}
}
if (n==0) n=1;
return(variance/(double)n);
}
void renormalise_field(phi, psi, xy_in, variance)
/* renormalise variance of field */
double *phi[NX], *psi[NX], variance;
short int * xy_in[NX];
{
int i, j;
double stdv;
stdv = sqrt(variance);
for (i=1; i<NX; i++)
for (j=1; j<NY; j++)
{
if (xy_in[i][j])
{
phi[i][j] = phi[i][j]/stdv;
psi[i][j] = psi[i][j]/stdv;
}
}
}
void animation()
{
double time, scale, dx, var;
double *phi[NX], *psi[NX];
short int *xy_in[NX];
int i, j, s;
/* Since NX and NY are big, it seemed wiser to use some memory allocation here */
for (i=0; i<NX; i++)
{
phi[i] = (double *)malloc(NY*sizeof(double));
psi[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
dx = (XMAX-XMIN)/((double)NX);
intstep = DT/(dx*dx*HBAR);
intstep1 = DT/(dx*HBAR);
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(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);
if (SCALE)
{
var = compute_variance(phi,psi, xy_in);
scale = sqrt(1.0 + var);
renormalise_field(phi, psi, xy_in, var);
}
blank();
glColor3f(0.0, 0.0, 0.0);
glutSwapBuffers();
sleep(SLEEP1);
for (i=0; i<=NSTEPS; i++)
{
/* compute the variance of the field to adjust color scheme */
/* the color depends on the field divided by sqrt(1 + variance) */
if (SCALE)
{
var = compute_variance(phi,psi, xy_in);
scale = sqrt(1.0 + var);
// printf("Norm: %5lg\t Scaling factor: %5lg\n", var, scale);
renormalise_field(phi, psi, xy_in, var);
}
else scale = 1.0;
draw_wave(phi, psi, xy_in, scale, i);
for (j=0; j<NVID; j++) evolve_wave(phi, psi, xy_in);
draw_billiard();
glutSwapBuffers();
if (MOVIE)
{
save_frame();
/* 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 */
if (i % PAUSE == PAUSE - 1)
{
printf("Making a short pause\n");
sleep(PSLEEP);
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]);
}
}
void display(void)
{
glPushMatrix();
blank();
glutSwapBuffers();
blank();
glutSwapBuffers();
animation();
sleep(SLEEP2);
glPopMatrix();
glutDestroyWindow(glutGetWindow());
}
int main(int argc, char** argv)
{
glutInit(&argc, argv);
glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH);
glutInitWindowSize(WINWIDTH,WINHEIGHT);
glutCreateWindow("Schrodinger equation in a planar domain");
init();
glutDisplayFunc(display);
glutMainLoop();
return 0;
}

861
sub_part_billiard.c
View File

@ -1,5 +1,74 @@
long int global_time = 0; /* counter to keep track of global time of simulation */
int nparticles=NPART; int nparticles=NPART;
/*********************/
/* some basic math */
/*********************/
double vabs(x) /* absolute value */
double x;
{
double res;
if (x<0.0) res = -x;
else res = x;
return(res);
}
double module2(x, y) /* Euclidean norm */
double x, y;
{
double m;
m = sqrt(x*x + y*y);
return(m);
}
double argument(x, y)
double x, y;
{
double alph;
if (x!=0.0)
{
alph = atan(y/x);
if (x<0.0)
alph += PI;
}
else
{
alph = PID;
if (y<0.0)
alph = PI*1.5;
}
// if (alph < 0.0) alph += DPI;
return(alph);
}
int polynome(a, b, c, r)
double a, b, c, r[2];
{
double delta, rdelta;
int im = 1;
delta = b*b - 4*a*c;
if (delta<0.0)
{
/* printf("ca deconne!");*/
rdelta = 0.0;
im = 0;
}
else rdelta = sqrt(delta);
r[0] = (-b + rdelta)/(2.0*a);
r[1] = (-b - rdelta)/(2.0*a);
return(im);
}
/*********************/ /*********************/
/* Graphics routines */ /* Graphics routines */
/*********************/ /*********************/
@ -125,7 +194,14 @@ double rgb[3];
void blank() void blank()
{ {
if (BLACK) glClearColor(0.0, 0.0, 0.0, 1.0); double rgb[3];
if (COLOR_OUTSIDE)
{
hsl_to_rgb(OUTER_COLOR, 0.9, 0.15, rgb);
glClearColor(rgb[0], rgb[1], rgb[2], 1.0);
}
else if (BLACK) glClearColor(0.0, 0.0, 0.0, 1.0);
else glClearColor(1.0, 1.0, 1.0, 1.0); else glClearColor(1.0, 1.0, 1.0, 1.0);
glClear(GL_COLOR_BUFFER_BIT); glClear(GL_COLOR_BUFFER_BIT);
} }
@ -163,15 +239,59 @@ void write_text( double x, double y, char *st)
} }
void paint_billiard_interior() /* points billiard interior, for use before draw_conf */
{
double x0, x, y, phi, r = 0.01, alpha, dphi, omega;
int i, k, c;
glLineWidth(4);
glEnable(GL_LINE_SMOOTH);
switch (B_DOMAIN) {
case (D_POLYGON):
{
omega = DPI/((double)NPOLY);
if (PAINT_INT)
{
if (BLACK) glColor3f(1.0, 1.0, 1.0);
else glColor3f(0.0, 0.0, 0.0);
glBegin(GL_TRIANGLE_FAN);
glVertex2d(0.0, 0.0);
for (i=0; i<=NPOLY; i++)
{
x = cos(i*omega + APOLY*PID);
y = sin(i*omega + APOLY*PID);
glVertex2d(x, y);
x = cos((i+1)*omega + APOLY*PID);
y = sin((i+1)*omega + APOLY*PID);
glVertex2d(x, y);
}
glEnd();
}
break;
}
default:
{
}
}
}
void draw_billiard() /* draws the billiard boundary */ void draw_billiard() /* draws the billiard boundary */
{ {
double x0, x, y, phi, r = 0.01, alpha, dphi; double x0, x, y, phi, r = 0.01, alpha, dphi, omega, x1, y1, x2, beta2, angle, s;
int i; int i, j, k, c;
if (BLACK) glColor3f(1.0, 1.0, 1.0); if (PAINT_INT) glColor3f(0.5, 0.5, 0.5);
else glColor3f(0.0, 0.0, 0.0); else
glLineWidth(4); {
if (BLACK) glColor3f(1.0, 1.0, 1.0);
else glColor3f(0.0, 0.0, 0.0);
}
glLineWidth(BILLIARD_WIDTH);
glEnable(GL_LINE_SMOOTH); glEnable(GL_LINE_SMOOTH);
@ -324,6 +444,142 @@ void draw_billiard() /* draws the billiard boundary */
glEnd(); glEnd();
break; break;
} }
case (D_ANNULUS):
{
/* color inner circle */
glColor3f(0.5, 0.5, 0.5);
glBegin(GL_TRIANGLE_FAN);
glVertex2d(MU, 0.0);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi) + MU;
y = LAMBDA*sin(phi);
glVertex2d(x, y);
}
glEnd();
/* color outer domain */
glColor3f(0.2, 0.2, 0.2);
glBegin(GL_TRIANGLE_FAN);
glVertex2d(XMAX, YMAX);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*PID/(double)NSEG;
x = cos(phi);
y = sin(phi);
glVertex2d(x, y);
}
glEnd();
glBegin(GL_TRIANGLE_FAN);
glVertex2d(XMIN, YMAX);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*PID/(double)NSEG + PID;
x = cos(phi);
y = sin(phi);
glVertex2d(x, y);
}
glEnd();
glBegin(GL_TRIANGLE_FAN);
glVertex2d(XMIN, YMIN);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*PID/(double)NSEG + PI;
x = cos(phi);
y = sin(phi);
glVertex2d(x, y);
}
glEnd();
glBegin(GL_TRIANGLE_FAN);
glVertex2d(XMAX, YMIN);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*PID/(double)NSEG + 3.0*PID;
x = cos(phi);
y = sin(phi);
glVertex2d(x, y);
}
glEnd();
glBegin(GL_TRIANGLES);
glVertex2d(XMAX, YMAX);
glVertex2d(1.0, 0.0);
glVertex2d(XMAX, YMIN);
glVertex2d(XMAX, YMIN);
glVertex2d(0.0, -1.0);
glVertex2d(XMIN, YMIN);
glVertex2d(XMIN, YMIN);
glVertex2d(-1.0, 0.0);
glVertex2d(XMIN, YMAX);
glVertex2d(XMIN, YMAX);
glVertex2d(0.0, 1.0);
glVertex2d(XMAX, YMAX);
glEnd();
/* draw circles */
if (BLACK) glColor3f(1.0, 1.0, 1.0);
else glColor3f(0.0, 0.0, 0.0);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi) + MU;
y = LAMBDA*sin(phi);
glVertex2d(x, y);
}
glEnd ();
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = cos(phi);
y = sin(phi);
glVertex2d(x, y);
}
glEnd ();
break;
}
case (D_POLYGON):
{
omega = DPI/((double)NPOLY);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NPOLY; i++)
{
x = cos(i*omega + APOLY*PID);
y = sin(i*omega + APOLY*PID);
glVertex2d(x, y);
}
glEnd ();
break;
}
case (D_REULEAUX):
{
omega = DPI/((double)NPOLY);
beta2 = asin(sin(omega*0.5)/LAMBDA);
if (LAMBDA > 0.0) x2 = cos(omega*0.5) + sqrt(LAMBDA*LAMBDA - sin(omega*0.5)*sin(omega*0.5));
else x2 = cos(omega*0.5) - sqrt(LAMBDA*LAMBDA - sin(omega*0.5)*sin(omega*0.5));
glBegin(GL_LINE_STRIP);
for (i=0; i<=NPOLY; i++)
{
for (j=0; j<NSEG; j++)
{
s = 2.0*(((double)j/(double)NSEG)-0.5)*beta2;
x1 = x2 - LAMBDA*cos(s);
y1 = LAMBDA*sin(s);
angle = i*omega + APOLY*PID;
x = cos(angle)*x1 - sin(angle)*y1;
y = sin(angle)*x1 + cos(angle)*y1;
glVertex2d(x, y);
}
}
glEnd ();
break;
}
default: default:
{ {
printf("Function draw_billiard not defined for this billiard \n"); printf("Function draw_billiard not defined for this billiard \n");
@ -332,74 +588,6 @@ void draw_billiard() /* draws the billiard boundary */
} }
/*********************/
/* some basic math */
/*********************/
double vabs(x) /* absolute value */
double x;
{
double res;
if (x<0.0) res = -x;
else res = x;
return(res);
}
double module2(x, y) /* Euclidean norm */
double x, y;
{
double m;
m = sqrt(x*x + y*y);
return(m);
}
double argument(x, y)
double x, y;
{
double alph;
if (x!=0.0)
{
alph = atan(y/x);
if (x<0.0)
alph += PI;
}
else
{
alph = PID;
if (y<0.0)
alph = PI*1.5;
}
// if (alph < 0.0) alph += DPI;
return(alph);
}
int polynome(a, b, c, r)
double a, b, c, r[2];
{
double delta, rdelta;
int im = 1;
delta = b*b - 4*a*c;
if (delta<0.0)
{
/* printf("ca deconne!");*/
rdelta = 0.0;
im = 0;
}
else rdelta = sqrt(delta);
r[0] = (-b + rdelta)/(2.0*a);
r[1] = (-b - rdelta)/(2.0*a);
return(im);
}
/*********************************/ /*********************************/
/* computation of the collisions */ /* computation of the collisions */
@ -423,9 +611,25 @@ double conf[8];
printf("s = %.3lg, u = %.3lg, t = %.3lg, L = %.3lg, x0 = %.3lg, y0 = %.3lg, x1 = %.3lg, y1 = %.3lg\n", conf[0], conf[1]/PI, conf[2], conf[3], conf[4], conf[5], conf[6], conf[7]); printf("s = %.3lg, u = %.3lg, t = %.3lg, L = %.3lg, x0 = %.3lg, y0 = %.3lg, x1 = %.3lg, y1 = %.3lg\n", conf[0], conf[1]/PI, conf[2], conf[3], conf[4], conf[5], conf[6], conf[7]);
} }
void print_config_23(conf) /* for debugging purposes */
double conf[8];
{
printf("t = %.8f, L = %.8f\n", conf[2], conf[3]);
}
void print_colors(color) /* for debugging purposes */
int color[NPARTMAX];
{
int i;
for (i=0; i<NPART; i++) printf("%i ", color[i]);
printf("\n");
}
/****************************************************************************************/
/* rectangle billiard */ /* rectangle billiard */
/****************************************************************************************/
int pos_rectangle(conf, pos, alpha) int pos_rectangle(conf, pos, alpha)
@ -696,7 +900,9 @@ double conf[8];
return(1); return(1);
} }
/* stadium billiard */ /****************************************************************************************/
/* stadium billiard */
/****************************************************************************************/
int pos_stade(conf, pos, alpha) int pos_stade(conf, pos, alpha)
/* determine position on boundary of stadium */ /* determine position on boundary of stadium */
@ -869,7 +1075,9 @@ double conf[8];
return(c); return(c);
} }
/* Sinai billiard */ /****************************************************************************************/
/* Sinai billiard */
/****************************************************************************************/
int pos_sinai(conf, pos, alpha) int pos_sinai(conf, pos, alpha)
/* determine position on boundary of Sinai billiard */ /* determine position on boundary of Sinai billiard */
@ -1066,7 +1274,9 @@ double conf[8];
} }
/* triangle billiard */ /****************************************************************************************/
/* triangle billiard */
/****************************************************************************************/
int pos_triangle(conf, pos, alpha) int pos_triangle(conf, pos, alpha)
@ -1188,8 +1398,395 @@ double conf[8];
} }
/****************************************************************************************/
/* annulus billiard */
/****************************************************************************************/
/* general billiard */ int pos_annulus(conf, pos, alpha)
/* determine position on boundary of annulus */
double conf[2], pos[2], *alpha;
{
double s, theta, psi0, psi, s1, s2, s3, s4;
s = conf[0];
theta = conf[1];
if (conf[1] < 0.0) conf[1] += DPI;
if (conf[0] < DPI) /* inner circle */
{
pos[0] = LAMBDA*cos(conf[0]) + MU;
pos[1] = LAMBDA*sin(conf[0]);
theta = PID + conf[0];
*alpha = theta - conf[1];
return(0);
}
else /* outer circle */
{
pos[0] = cos(conf[0]);
pos[1] = sin(conf[0]);
theta = argument(-pos[1],pos[0]);
*alpha = theta + conf[1];
return(1);
}
}
int vannulus_xy(config, alpha, pos)
/* determine initial configuration for start at point pos = (x,y) */
double config[8], alpha, pos[2];
{
double l, s0, c0, t, t1, x, y, x1, y1, a, b, delta, res[2], s, r;
double psi, lam2, margin = 1.0e-14, theta;
// double psi, lam2, margin = 1.0e-14, theta;
int c, intb=1, intc, i;
/* initial position and velocity */
c0 = cos(alpha);
s0 = sin(alpha);
/* intersection with inner circle, using parametric equation of line */
b = (pos[0]-MU)*c0 + pos[1]*s0;
a = (pos[0]-MU)*(pos[0]-MU) + pos[1]*pos[1] - LAMBDA*LAMBDA;
delta = b*b - a;
if ((delta > margin)&&(a > margin))
{
t = - b - sqrt(delta);
x1 = pos[0] + t*c0;
y1 = pos[1] + t*s0;
s = argument(x1-MU,y1);
while (s<0.0) s += DPI;
while (s>=DPI) s -= DPI;
config[0] = s;
config[1] = 3.0*PID - s + alpha;
c = 0;
}
else /* intersection with outer circle, using parametric equation of line */
{
b = pos[0]*c0 + pos[1]*s0;
a = pos[0]*pos[0] + pos[1]*pos[1] - 1.0;
t = (-b+sqrt(b*b - a));
x1 = pos[0] + t*c0;
y1 = pos[1] + t*s0;
/* parameter of intersection with outer circle */
config[0] = argument(x1, y1);
while (config[0] < DPI) config[0] += DPI;
while (config[0] >= 2.0*DPI) config[0] -= DPI;
/* computation of outgoing angle after collision with outer circle */
theta = argument(-y1,x1);
config[1] = theta - alpha;
// while (config[1] < 0.0) config[1] += DPI;
// while (config[1] > DPI) config[1] -= DPI;
c = 1;
}
if (config[1] < 0.0) config[1] += DPI;
// config[2] = 1.0e-12; /* running time */
config[2] = 0.0; /* running time */
config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */
config[4] = pos[0]; /* start position */
config[5] = pos[1];
config[6] = x1; /* position of collision */
config[7] = y1;
return(c);
}
int vannulus(config)
/* determine initial configuration when starting from boundary */
double config[8];
{
double pos[2], alpha;
int c;
c = pos_annulus(config, pos, &alpha);
vannulus_xy(config, alpha, pos);
return(c);
}
/****************************************************************************************/
/* polygonal billiard */
/****************************************************************************************/
int pos_polygon(conf, pos, alpha)
/* determine position on boundary of polygon */
/* conf[0] is arclength on boundary */
double conf[2], pos[2], *alpha;
{
double s, theta, omega, length, s1, angle, x, y;
int c;
s = conf[0];
theta = conf[1];
omega = DPI/((double)NPOLY);
length = 2.0*sin(0.5*omega);
c = (int)(s/length); /* side of polygon */
s1 = s - ((double)c)*length;
x = 1.0 + (cos(omega) - 1.0)*s1/length;
y = sin(omega)*s1/length;
angle = (double)c*omega + PID*APOLY;
pos[0] = x*cos(angle) - y*sin(angle);
pos[1] = x*sin(angle) + y*cos(angle);
*alpha = (0.5 + (double)c)*omega + theta + PID*(1.0 + APOLY);
return(c);
}
int vpolygon_xy(config, alpha, pos)
/* determine initial configuration for start at point pos = (x,y) */
double config[8], alpha, pos[2];
{
double s, theta, omega, length, rlength, s1, rangle, x, y, xp, yp, x1, y1, ca, sa;
int k, c, intb=1, intc, i;
/* dimensions/angles of polygon */
omega = DPI/((double)NPOLY);
length = 2.0*sin(0.5*omega);
rlength = cos(0.5*omega);
for (k=0; k<NPOLY; k++)
{
/* rotate position so that kth side is vertical */
rangle = (0.5 + (double)k)*omega + APOLY*PID;
theta = alpha - rangle;
if ((cos(theta) > 0.0)&&(intb))
{
ca = cos(rangle);
sa = sin(rangle);
x = pos[0]*ca + pos[1]*sa;
y = -pos[0]*sa + pos[1]*ca;
xp = rlength;
yp = y + (xp-x)*tan(theta);
if (vabs(yp) < 0.5*length)
{
/* rotate back */
x1 = xp*ca - yp*sa;
y1 = xp*sa + yp*ca;
intb = 0;
c = k;
config[0] = ((double)k + 0.5)*length + yp;
config[1] = PID - theta;
}
}
}
if (config[1] < 0.0) config[1] += DPI;
config[2] = 0.0; /* running time */
config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */
config[4] = pos[0]; /* start position */
config[5] = pos[1];
config[6] = x1; /* position of collision */
config[7] = y1;
return(c);
}
int vpolygon(config)
/* determine initial configuration when starting from boundary */
double config[8];
{
double pos[2], alpha;
int c;
c = pos_polygon(config, pos, &alpha);
vpolygon_xy(config, alpha, pos);
return(c);
}
/****************************************************************************************/
/* Reuleaux-type and star-shaped billiard */
/****************************************************************************************/
int pos_reuleaux(conf, pos, alpha)
/* determine position on boundary of polygon */
/* conf[0] is arclength on boundary */
double conf[2], pos[2], *alpha;
{
double s, theta, omega2, beta2, beta, s1, angle, x2, x, y;
int c;
s = conf[0];
theta = conf[1];
omega2 = PI/((double)NPOLY);
beta2 = asin(sin(omega2)/vabs(LAMBDA));
beta = beta2*2.0;
c = (int)(s/beta); /* side of shape */
s1 = s - ((double)c)*beta;
if (LAMBDA > 0.0) x2 = cos(omega2) + sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
else x2 = cos(omega2) - sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
x = x2 - LAMBDA*cos(s1 - beta2);
y = LAMBDA*sin(s1 - beta2);
angle = 2.0*((double)c)*omega2 + PID*APOLY;
pos[0] = x*cos(angle) - y*sin(angle);
pos[1] = x*sin(angle) + y*cos(angle);
*alpha = PID - s1 + beta2 + theta + 2.0*(double)c*omega2 + APOLY*PID;
// printf("alpha = %.5lg\t", *alpha);
return(c);
}
int vreuleaux_xy(config, alpha, pos)
/* determine initial configuration for start at point pos = (x,y) */
double config[8], alpha, pos[2];
{
double s, theta, omega2, beta, s1, rangle, x, y, x1[NPOLY], y1[NPOLY], xi, yi, t, x2;
double ca, sa, a, b, margin = 1.0e-14, tmin, tval[NPOLY], tempconf[NPOLY][2];
int k, c, intb=1, intc, i, nt = 0, cval[NPOLY], ntmin;
/* dimensions/angles of polygon */
omega2 = PI/((double)NPOLY);
beta = 2.0*asin(sin(omega2)/vabs(LAMBDA));
// printf("beta = %.5lg\n", beta);
if (LAMBDA > 0.0) x2 = cos(omega2) + sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
else x2 = cos(omega2) - sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
// printf("x2 = %.5lg\n", x2);
for (k=0; k<NPOLY; k++)
{
/* rotate position so that kth side is vertical */
rangle = 2.0*(double)k*omega2 + APOLY*PID;
theta = alpha - rangle;
// if ((intb)) /* check if condition is ok */
// if ((cos(theta) > 0.0)&&(intb)) /* check if condition is ok */
{
ca = cos(rangle);
sa = sin(rangle);
// printf("theta = %.5lg\n", theta);
// printf("rangle = %.5lg x0 = %.5lg y0 = %.5lg \n", rangle, pos[0], pos[1]);
x = pos[0]*ca + pos[1]*sa;
y = -pos[0]*sa + pos[1]*ca;
// printf("x = %.5lg\t y = %.5lg\n", x, y);
a = (x-x2)*cos(theta) + y*sin(theta);
b = (x-x2)*(x-x2) + y*y - LAMBDA*LAMBDA;
// printf("a = %.5lg\t b = %.5lg\n", a, b);
if (a*a - b > margin)
{
// t = vabs(a) - sqrt(a*a - b);
if (LAMBDA > 0.0) t = -a - sqrt(a*a - b);
else t = -a + sqrt(a*a - b);
xi = x + t*cos(theta);
yi = y + t*sin(theta);
// printf("t = %.5lg\t xi = %.5lg\t yi = %.5lg\n", t, xi, yi);
if ((t > margin)&&(vabs(yi) <= sin(omega2)))
{
cval[nt] = k;
tval[nt] = t;
/* rotate back */
x1[nt] = xi*ca - yi*sa;
y1[nt] = xi*sa + yi*ca;
// intb = 0;
// c = k;
tempconf[nt][0] = ((double)k + 0.5)*beta + asin(yi/LAMBDA);
tempconf[nt][1] = PID - asin(yi/LAMBDA) - theta;
// tempconf[nt][0] = ((double)k + 0.5)*beta + asin(yi/vabs(LAMBDA));
// tempconf[nt][1] = PID - asin(yi/vabs(LAMBDA)) - theta;
nt++;
}
}
}
}
// printf("nt = %i\n", nt);
/* find earliest intersection */
tmin = tval[0];
ntmin = 0;
for (i=1; i<nt; i++)
if (tval[i] < tmin)
{
tmin = tval[i];
ntmin = i;
}
config[0] = tempconf[ntmin][0];
config[1] = tempconf[ntmin][1];
c = cval[ntmin];
// printf("nt = %i\t ntmin = %i \tcmin = %i\n", nt, ntmin, c);
if (config[1] < 0.0) config[1] += DPI;
config[2] = 0.0; /* running time */
config[3] = module2(x1[ntmin]-pos[0], y1[ntmin]-pos[1]); /* distance to collision */
config[4] = pos[0]; /* start position */
config[5] = pos[1];
config[6] = x1[ntmin]; /* position of collision */
config[7] = y1[ntmin];
return(c);
}
int vreuleaux(config)
/* determine initial configuration when starting from boundary */
double config[8];
{
double pos[2], alpha;
int c;
c = pos_reuleaux(config, pos, &alpha);
vreuleaux_xy(config, alpha, pos);
return(c);
}
/****************************************************************************************/
/* general billiard */
/****************************************************************************************/
int pos_billiard(conf, pos, alpha) int pos_billiard(conf, pos, alpha)
/* determine initial configuration for start at point pos = (x,y) */ /* determine initial configuration for start at point pos = (x,y) */
@ -1221,9 +1818,24 @@ double conf[8];
return(pos_triangle(conf, pos, &alpha)); return(pos_triangle(conf, pos, &alpha));
break; break;
} }
case (D_ANNULUS):
{
return(pos_annulus(conf, pos, &alpha));
break;
}
case (D_POLYGON):
{
return(pos_polygon(conf, pos, &alpha));
break;
}
case (D_REULEAUX):
{
return(pos_reuleaux(conf, pos, &alpha));
break;
}
default: default:
{ {
printf("Function vbilliard_xy not defined for this billiard \n"); printf("Function pos_billiard not defined for this billiard \n");
} }
} }
} }
@ -1260,6 +1872,21 @@ double conf[8];
return(vtriangle_xy(config, alpha, pos)); return(vtriangle_xy(config, alpha, pos));
break; break;
} }
case (D_ANNULUS):
{
return(vannulus_xy(config, alpha, pos));
break;
}
case (D_POLYGON):
{
return(vpolygon_xy(config, alpha, pos));
break;
}
case (D_REULEAUX):
{
return(vreuleaux_xy(config, alpha, pos));
break;
}
default: default:
{ {
printf("Function vbilliard_xy not defined for this billiard \n"); printf("Function vbilliard_xy not defined for this billiard \n");
@ -1312,9 +1939,30 @@ double conf[8];
return(vtriangle(config, alpha, pos)); return(vtriangle(config, alpha, pos));
break; break;
} }
case (D_ANNULUS):
{
c = pos_annulus(config, pos, &alpha);
return(vannulus(config, alpha, pos));
break;
}
case (D_POLYGON):
{
c = pos_polygon(config, pos, &alpha);
return(vpolygon(config, alpha, pos));
break;
}
case (D_REULEAUX):
{
c = pos_reuleaux(config, pos, &alpha);
return(vreuleaux(config, alpha, pos));
break;
}
default: default:
{ {
printf("Function vbilliard_xy not defined for this billiard \n"); printf("Function vbilliard not defined for this billiard \n");
} }
} }
} }
@ -1323,7 +1971,8 @@ double conf[8];
/* returns 1 if (x,y) represents a point in the billiard */ /* returns 1 if (x,y) represents a point in the billiard */
double x, y; double x, y;
{ {
double l2, r2; double l2, r2, omega, omega2, c, angle, x1, y1, x2, co, so;
int condition, k;
switch (B_DOMAIN) { switch (B_DOMAIN) {
case D_RECTANGLE: case D_RECTANGLE:
@ -1369,6 +2018,56 @@ double conf[8];
else return(0); else return(0);
break; break;
} }
case D_ANNULUS:
{
l2 = LAMBDA*LAMBDA;
r2 = x*x + y*y;
if ((r2 > l2)&&(r2 < 1.0)) return(1);
else return(0);
break;
}
case D_POLYGON:
{
condition = 1;
omega = DPI/((double)NPOLY);
c = cos(omega*0.5);
for (k=0; k<NPOLY; k++)
{
angle = APOLY*PID + (k+0.5)*omega;
condition = condition*(x*cos(angle) + y*sin(angle) < c);
}
return(condition);
break;
}
case D_REULEAUX:
{
condition = 1;
omega2 = PI/((double)NPOLY);
co = cos(omega2);
so = sin(omega2);
if (LAMBDA > 0.0) x2 = co + sqrt(LAMBDA*LAMBDA - so*so);
else x2 = co - sqrt(LAMBDA*LAMBDA - so*so);
for (k=0; k<NPOLY; k++)
{
angle = 2.0*(double)k*omega2 + APOLY*PID;
x1 = x*cos(angle) + y*sin(angle);
y1 = -x*sin(angle) + y*cos(angle);
if (LAMBDA > 0.0) condition = condition*((x1-x2)*(x1-x2) + y1*y1 > LAMBDA*LAMBDA);
else condition = condition*((x1-x2)*(x1-x2) + y1*y1 < LAMBDA*LAMBDA);
// if (!condition)
// {
// printf("x = %.5lg \t y = %.5lg \t x1 = %.5lg \t y1 = %.5lg \t angle = %.5lg \n", x, y, x1, y1, angle);
// printf("k = %i \t condition = %i\n", k, condition);
// sleep(1);
// }
}
return(condition);
break;
}
/* D_REULEAUX : distance to all centers of arcs should be larger than LAMBDA */
default: default:
{ {
printf("Function ij_in_billiard not defined for this billiard \n"); printf("Function ij_in_billiard not defined for this billiard \n");

718
sub_wave.c Normal file
View File

@ -0,0 +1,718 @@
/*********************/
/* Graphics routines */
/*********************/
int writetiff(char *filename, char *description, int x, int y, int width, int height, int compression)
{
TIFF *file;
GLubyte *image, *p;
int i;
file = TIFFOpen(filename, "w");
if (file == NULL)
{
return 1;
}
image = (GLubyte *) malloc(width * height * sizeof(GLubyte) * 3);
/* OpenGL's default 4 byte pack alignment would leave extra bytes at the
end of each image row so that each full row contained a number of bytes
divisible by 4. Ie, an RGB row with 3 pixels and 8-bit componets would
be laid out like "RGBRGBRGBxxx" where the last three "xxx" bytes exist
just to pad the row out to 12 bytes (12 is divisible by 4). To make sure
the rows are packed as tight as possible (no row padding), set the pack
alignment to 1. */
glPixelStorei(GL_PACK_ALIGNMENT, 1);
glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
TIFFSetField(file, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT);
TIFFSetField(file, TIFFTAG_SAMPLESPERPIXEL, 3);
TIFFSetField(file, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
TIFFSetField(file, TIFFTAG_ROWSPERSTRIP, 1);
TIFFSetField(file, TIFFTAG_IMAGEDESCRIPTION, description);
p = image;
for (i = height - 1; i >= 0; i--)
{
// if (TIFFWriteScanline(file, p, height - i - 1, 0) < 0)
if (TIFFWriteScanline(file, p, i, 0) < 0)
{
free(image);
TIFFClose(file);
return 1;
}
p += width * sizeof(GLubyte) * 3;
}
TIFFClose(file);
return 0;
}
void init() /* initialisation of window */
{
glLineWidth(3);
glClearColor(0.0, 0.0, 0.0, 1.0);
glClear(GL_COLOR_BUFFER_BIT);
// glOrtho(XMIN, XMAX, YMIN, YMAX , -1.0, 1.0);
glOrtho(0.0, NX, 0.0, NY, -1.0, 1.0);
}
void hsl_to_rgb(h, s, l, rgb) /* color conversion from HSL to RGB */
/* h = hue, s = saturation, l = luminosity */
double h, s, l, rgb[3];
{
double c = 0.0, m = 0.0, x = 0.0;
c = (1.0 - fabs(2.0 * l - 1.0)) * s;
m = 1.0 * (l - 0.5 * c);
x = c * (1.0 - fabs(fmod(h / 60.0, 2) - 1.0));
if (h >= 0.0 && h < 60.0)
{
rgb[0] = c+m; rgb[1] = x+m; rgb[2] = m;
}
else if (h < 120.0)
{
rgb[0] = x+m; rgb[1] = c+m; rgb[2] = m;
}
else if (h < 180.0)
{
rgb[0] = m; rgb[1] = c+m; rgb[2] = x+m;
}
else if (h < 240.0)
{
rgb[0] = m; rgb[1] = x+m; rgb[2] = c+m;
}
else if (h < 300.0)
{
rgb[0] = x+m; rgb[1] = m; rgb[2] = c+m;
}
else if (h < 360.0)
{
rgb[0] = c+m; rgb[1] = m; rgb[2] = x+m;
}
else
{
rgb[0] = m; rgb[1] = m; rgb[2] = m;
}
}
double color_amplitude(value, scale, time)
/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */
double value, scale;
int time;
{
return(tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION)));
}
void color_scheme(scheme, value, scale, time, rgb) /* color scheme */
double value, scale;
int scheme, time;
double rgb[3];
{
double hue, y, r, amplitude;
int intpart;
/* saturation = r, luminosity = y */
switch (scheme) {
case C_LUM:
{
hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = LUMMEAN + amplitude*LUMAMP;
intpart = (int)y;
y -= (double)intpart;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case C_HUE:
{
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = 0.5;
hue = HUEMEAN + amplitude*HUEAMP;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
hsl_to_rgb(hue, r, y, rgb);
break;
}
}
}
void blank()
{
if (BLACK) glClearColor(0.0, 0.0, 0.0, 1.0);
else glClearColor(1.0, 1.0, 1.0, 1.0);
glClear(GL_COLOR_BUFFER_BIT);
}
void save_frame()
{
static int counter = 0;
char *name="wave.", n2[100];
char format[6]=".%05i";
counter++;
// printf (" p2 counter = %d \n",counter);
strcpy(n2, name);
sprintf(strstr(n2,"."), format, counter);
strcat(n2, ".tif");
printf(" saving frame %s \n",n2);
writetiff(n2, "Wave equation in a planar domain", 0, 0,
WINWIDTH, WINHEIGHT, COMPRESSION_LZW);
}
void write_text( double x, double y, char *st)
{
int l,i;
l=strlen( st ); // see how many characters are in text string.
glRasterPos2d( x, y); // location to start printing text
for( i=0; i < l; i++) // loop until i is greater then l
{
glutBitmapCharacter(GLUT_BITMAP_TIMES_ROMAN_24, st[i]); // Print a character on the screen
// glutBitmapCharacter(GLUT_BITMAP_8_BY_13, st[i]); // Print a character on the screen
}
}
/*********************/
/* some basic math */
/*********************/
double vabs(x) /* absolute value */
double x;
{
double res;
if (x<0.0) res = -x;
else res = x;
return(res);
}
double module2(x, y) /* Euclidean norm */
double x, y;
{
double m;
m = sqrt(x*x + y*y);
return(m);
}
double argument(x, y)
double x, y;
{
double alph;
if (x!=0.0)
{
alph = atan(y/x);
if (x<0.0)
alph += PI;
}
else
{
alph = PID;
if (y<0.0)
alph = PI*1.5;
}
return(alph);
}
/*********************/
/* drawing routines */
/*********************/
/* The billiard boundary is drawn in (x,y) coordinates */
/* However for the grid points, we use integer coordinates (i,j) */
/* GL would allow to always work in (x,y) coordinates but using both */
/* sets of coordinates decreases number of double computations when */
/* drawing the field */
void xy_to_ij(x, y, ij)
/* convert (x,y) position to (i,j) in table representing wave */
double x, y;
int ij[2];
{
double x1, y1;
x1 = (x - XMIN)/(XMAX - XMIN);
y1 = (y - YMIN)/(YMAX - YMIN);
ij[0] = (int)(x1 * (double)NX);
ij[1] = (int)(y1 * (double)NY);
}
void xy_to_pos(x, y, pos)
/* convert (x,y) position to double-valued position in table representing wave */
double x, y, pos[2];
{
double x1, y1;
x1 = (x - XMIN)/(XMAX - XMIN);
y1 = (y - YMIN)/(YMAX - YMIN);
pos[0] = x1 * (double)NX;
pos[1] = y1 * (double)NY;
}
void ij_to_xy(i, j, xy)
/* convert (i,j) position in table representing wave to (x,y) */
int i, j;
double xy[2];
{
double x1, y1;
xy[0] = XMIN + ((double)i)*(XMAX-XMIN)/((double)NX);
xy[1] = YMIN + ((double)j)*(YMAX-YMIN)/((double)NY);
}
int xy_in_billiard(x, y)
/* returns 1 if (x,y) represents a point in the billiard */
double x, y;
{
double l2, r2, omega, c, angle, z;
int i, k, k1, k2, condition;
switch (B_DOMAIN) {
case D_RECTANGLE:
{
if ((vabs(x) <LAMBDA)&&(vabs(y) < 1.0)) return(1);
else return(0);
break;
}
case D_ELLIPSE:
{
if (x*x/(LAMBDA*LAMBDA) + y*y < 1.0) return(1);
else return(0);
break;
}
case D_STADIUM:
{
if ((x > -0.5*LAMBDA)&&(x < 0.5*LAMBDA)&&(y > -1.0)&&(y < 1.0)) return(1);
else if (module2(x+0.5*LAMBDA, y) < 1.0) return(1);
else if (module2(x-0.5*LAMBDA, y) < 1.0) return(1);
else return(0);
break;
}
case D_SINAI:
{
if (x*x + y*y > LAMBDA*LAMBDA) return(1);
else return(0);
break;
}
case D_DIAMOND:
{
l2 = LAMBDA*LAMBDA;
r2 = l2 + (LAMBDA-1.0)*(LAMBDA-1.0);
if ((x*x + y*y < 1.0)&&((x-LAMBDA)*(x-LAMBDA) + (y-LAMBDA)*(y-LAMBDA) > r2)
&&((x-LAMBDA)*(x-LAMBDA) + (y+LAMBDA)*(y+LAMBDA) > r2)
&&((x+LAMBDA)*(x+LAMBDA) + (y-LAMBDA)*(y-LAMBDA) > r2)
&&((x+LAMBDA)*(x+LAMBDA) + (y+LAMBDA)*(y+LAMBDA) > r2)) return(1);
else return(0);
break;
}
case D_TRIANGLE:
{
if ((x>-LAMBDA)&&(y>-1.0)&&(LAMBDA*y+x<0.0)) return(1);
else return(0);
break;
}
case D_FLAT:
{
if (y > -LAMBDA) return(1);
else return(0);
break;
}
case D_ANNULUS:
{
l2 = LAMBDA*LAMBDA;
r2 = x*x + y*y;
if ((r2 > l2)&&(r2 < 1.0)) return(1);
else return(0);
}
case D_POLYGON:
{
condition = 1;
omega = DPI/((double)NPOLY);
c = cos(omega*0.5);
for (k=0; k<NPOLY; k++)
{
angle = APOLY*PID + (k+0.5)*omega;
condition = condition*(x*cos(angle) + y*sin(angle) < c);
}
// for (k=0; k<NPOLY; k++) condition = condition*(-x*sin((k+0.5)*omega) + y*cos((k+0.5)*omega) < c);
return(condition);
}
case D_YOUNG:
{
if ((x < -MU)||(x > MU)) return(1);
else if ((vabs(y-LAMBDA) < MU)||(vabs(y+LAMBDA) < MU)) return (1);
else return(0);
}
case D_GRATING:
{
k1 = -(int)((-YMIN)/LAMBDA);
k2 = (int)(YMAX/LAMBDA);
condition = 1;
for (i=k1; i<= k2; i++)
{
z = (double)i*LAMBDA;
condition = condition*(x*x + (y-z)*(y-z) > MU*MU);
}
// printf("x = %.3lg, y = %.3lg, k1 = %i, k2 = %i, condition = %i\n", x, y, k1, k2, condition);
return(condition);
}
case D_EHRENFEST:
{
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)));
}
default:
{
printf("Function ij_in_billiard not defined for this billiard \n");
return(0);
}
}
}
int ij_in_billiard(i, j)
/* returns 1 if (i,j) represents a point in the billiard */
int i, j;
{
double xy[2];
ij_to_xy(i, j, xy);
return(xy_in_billiard(xy[0], xy[1]));
}
void draw_billiard() /* draws the billiard boundary */
{
double x0, x, y, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z;
int i, j, k1, k2;
if (BLACK) glColor3f(1.0, 1.0, 1.0);
else glColor3f(0.0, 0.0, 0.0);
glLineWidth(5);
glEnable(GL_LINE_SMOOTH);
switch (B_DOMAIN) {
case (D_RECTANGLE):
{
glBegin(GL_LINE_LOOP);
xy_to_pos(LAMBDA, -1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(LAMBDA, 1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-LAMBDA, 1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-LAMBDA, -1.0, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
break;
}
case (D_ELLIPSE):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
/* draw foci */
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
x0 = sqrt(LAMBDA*LAMBDA-1.0);
glLineWidth(2);
glEnable(GL_LINE_SMOOTH);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = x0 + r*cos(phi);
y = r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = -x0 + r*cos(phi);
y = r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
}
break;
}
case D_STADIUM:
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = -PID + (double)i*PI/(double)NSEG;
x = 0.5*LAMBDA + cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = PID + (double)i*PI/(double)NSEG;
x = -0.5*LAMBDA + cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
break;
}
case D_SINAI:
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = LAMBDA*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
break;
}
case D_DIAMOND:
{
alpha = atan(1.0 - 1.0/LAMBDA);
dphi = (PID - 2.0*alpha)/(double)NSEG;
r = sqrt(LAMBDA*LAMBDA + (LAMBDA-1.0)*(LAMBDA-1.0));
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = alpha + (double)i*dphi;
x = -LAMBDA + r*cos(phi);
y = -LAMBDA + r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = alpha - PID + (double)i*dphi;
x = -LAMBDA + r*cos(phi);
y = LAMBDA + r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = alpha + PI + (double)i*dphi;
x = LAMBDA + r*cos(phi);
y = LAMBDA + r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = alpha + PID + (double)i*dphi;
x = LAMBDA + r*cos(phi);
y = -LAMBDA + r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
break;
}
case (D_TRIANGLE):
{
glBegin(GL_LINE_LOOP);
xy_to_pos(-LAMBDA, -1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(LAMBDA, -1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-LAMBDA, 1.0, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
break;
}
case (D_FLAT):
{
glBegin(GL_LINE_LOOP);
xy_to_pos(XMIN, -LAMBDA, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(XMAX, -LAMBDA, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
break;
}
case (D_ANNULUS):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = LAMBDA*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
break;
}
case (D_POLYGON):
{
omega = DPI/((double)NPOLY);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NPOLY; i++)
{
x = cos(i*omega + APOLY*PID);
y = sin(i*omega + APOLY*PID);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
break;
}
case (D_YOUNG):
{
glBegin(GL_LINE_STRIP);
xy_to_pos(-MU, YMIN, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-MU, -LAMBDA-MU, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(MU, -LAMBDA-MU, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(MU, YMIN, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
glBegin(GL_LINE_STRIP);
xy_to_pos(-MU, YMAX, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-MU, LAMBDA+MU, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(MU, LAMBDA+MU, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(MU, YMAX, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
glBegin(GL_LINE_LOOP);
xy_to_pos(-MU, -LAMBDA+MU, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-MU, LAMBDA-MU, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(MU, LAMBDA-MU, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(MU, -LAMBDA+MU, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
break;
}
case D_GRATING:
{
k1 = -(int)(-YMIN/LAMBDA);
k2 = (int)(YMAX/LAMBDA);
for (i=k1; i<= k2; i++)
{
z = (double)i*LAMBDA;
glBegin(GL_LINE_LOOP);
for (j=0; j<=NSEG; j++)
{
phi = (double)j*DPI/(double)NSEG;
x = MU*cos(phi);
y = z + MU*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
}
break;
}
case D_EHRENFEST:
{
alpha = asin(MU/LAMBDA);
x0 = 1.0 - sqrt(LAMBDA*LAMBDA - MU*MU);
dphi = 2.0*(PI-alpha)/((double)NSEG);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = -PI + alpha + (double)i*dphi;
x = 1.0 + LAMBDA*cos(phi);
y = LAMBDA*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = alpha + (double)i*dphi;
x = -1.0 + LAMBDA*cos(phi);
y = LAMBDA*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
break;
}
default:
{
printf("Function draw_billiard not defined for this billiard \n");
}
}
}

627
wave_billiard.c
View File

@ -2,7 +2,7 @@
/* */ /* */
/* Animation of wave equation in a planar domain */ /* Animation of wave equation in a planar domain */
/* */ /* */
/* N. Berglund, december 2012, april 2021 */ /* N. Berglund, december 2012, may 2021 */
/* */ /* */
/* UPDATE 24/04: distinction between damping and "elasticity" parameters */ /* UPDATE 24/04: distinction between damping and "elasticity" parameters */
/* UPDATE 27/04: new billiard shapes, bug in color scheme fixed */ /* UPDATE 27/04: new billiard shapes, bug in color scheme fixed */
@ -70,9 +70,13 @@
#define D_FLAT 6 /* flat interface */ #define D_FLAT 6 /* flat interface */
#define D_ANNULUS 7 /* annulus */ #define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */ #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 LAMBDA 0.3 /* parameter controlling the dimensions of domain */ #define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */ #define MU 0.05 /* parameter controlling the dimensions of domain */
#define NPOLY 8 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */
@ -84,7 +88,7 @@
#define COURANT 0.01 /* Courant number */ #define COURANT 0.01 /* Courant number */
#define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMA 0.0 /* damping factor in wave equation */
// #define GAMMA 5.0e-10 /* damping factor in wave equation */ // #define GAMMA 5.0e-10 /* damping factor in wave equation */
#define KAPPA 5.0e-7 /* "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-9 /* "elasticity" term enforcing oscillations */
// #define KAPPA 5.0e-8 /* "elasticity" term enforcing oscillations */ // #define KAPPA 5.0e-8 /* "elasticity" term enforcing oscillations */
/* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */ /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
@ -98,8 +102,8 @@
/* Parameters for length and speed of simulation */ /* Parameters for length and speed of simulation */
#define NSTEPS 4575 //7500 /* number of frames of movie */ #define NSTEPS 5000 /* number of frames of movie */
#define NVID 50 /* number of iterations between images displayed on screen */ #define NVID 25 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */ #define NSEG 100 /* number of segments of boundary */
#define PAUSE 1000 /* number of frames after which to pause */ #define PAUSE 1000 /* number of frames after which to pause */
@ -109,6 +113,8 @@
/* Color schemes */ /* Color schemes */
#define BLACK 1 /* background */
#define COLOR_SCHEME 1 /* choice of color scheme */ #define COLOR_SCHEME 1 /* choice of color scheme */
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */ #define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
@ -122,8 +128,8 @@
#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ #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 LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */
#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ #define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */
#define HUEMEAN 235.0 /* mean value of hue for color scheme C_HUE */ #define HUEMEAN 100.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP 60.0 /* amplitude of variation of hue for color scheme C_HUE */ #define HUEAMP 80.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 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 */ // #define HUEAMP 100.0 /* amplitude of variation of hue for color scheme C_HUE */
@ -133,372 +139,10 @@
#define DPI 6.283185307 #define DPI 6.283185307
#define PID 1.570796327 #define PID 1.570796327
#include "sub_wave.c"
double courant2; /* Courant parameter squared */ double courant2; /* Courant parameter squared */
/*********************/
/* Graphics routines */
/*********************/
int writetiff(char *filename, char *description, int x, int y, int width, int height, int compression)
{
TIFF *file;
GLubyte *image, *p;
int i;
file = TIFFOpen(filename, "w");
if (file == NULL)
{
return 1;
}
image = (GLubyte *) malloc(width * height * sizeof(GLubyte) * 3);
/* OpenGL's default 4 byte pack alignment would leave extra bytes at the
end of each image row so that each full row contained a number of bytes
divisible by 4. Ie, an RGB row with 3 pixels and 8-bit componets would
be laid out like "RGBRGBRGBxxx" where the last three "xxx" bytes exist
just to pad the row out to 12 bytes (12 is divisible by 4). To make sure
the rows are packed as tight as possible (no row padding), set the pack
alignment to 1. */
glPixelStorei(GL_PACK_ALIGNMENT, 1);
glReadPixels(x, y, width, height, GL_RGB, GL_UNSIGNED_BYTE, image);
TIFFSetField(file, TIFFTAG_IMAGEWIDTH, (uint32) width);
TIFFSetField(file, TIFFTAG_IMAGELENGTH, (uint32) height);
TIFFSetField(file, TIFFTAG_BITSPERSAMPLE, 8);
TIFFSetField(file, TIFFTAG_COMPRESSION, compression);
TIFFSetField(file, TIFFTAG_PHOTOMETRIC, PHOTOMETRIC_RGB);
TIFFSetField(file, TIFFTAG_ORIENTATION, ORIENTATION_BOTLEFT);
TIFFSetField(file, TIFFTAG_SAMPLESPERPIXEL, 3);
TIFFSetField(file, TIFFTAG_PLANARCONFIG, PLANARCONFIG_CONTIG);
TIFFSetField(file, TIFFTAG_ROWSPERSTRIP, 1);
TIFFSetField(file, TIFFTAG_IMAGEDESCRIPTION, description);
p = image;
for (i = height - 1; i >= 0; i--)
{
// if (TIFFWriteScanline(file, p, height - i - 1, 0) < 0)
if (TIFFWriteScanline(file, p, i, 0) < 0)
{
free(image);
TIFFClose(file);
return 1;
}
p += width * sizeof(GLubyte) * 3;
}
TIFFClose(file);
return 0;
}
void init() /* initialisation of window */
{
glLineWidth(3);
glClearColor(0.0, 0.0, 0.0, 1.0);
glClear(GL_COLOR_BUFFER_BIT);
// glOrtho(XMIN, XMAX, YMIN, YMAX , -1.0, 1.0);
glOrtho(0.0, NX, 0.0, NY, -1.0, 1.0);
}
void hsl_to_rgb(h, s, l, rgb) /* color conversion from HSL to RGB */
/* h = hue, s = saturation, l = luminosity */
double h, s, l, rgb[3];
{
double c = 0.0, m = 0.0, x = 0.0;
c = (1.0 - fabs(2.0 * l - 1.0)) * s;
m = 1.0 * (l - 0.5 * c);
x = c * (1.0 - fabs(fmod(h / 60.0, 2) - 1.0));
if (h >= 0.0 && h < 60.0)
{
rgb[0] = c+m; rgb[1] = x+m; rgb[2] = m;
}
else if (h < 120.0)
{
rgb[0] = x+m; rgb[1] = c+m; rgb[2] = m;
}
else if (h < 180.0)
{
rgb[0] = m; rgb[1] = c+m; rgb[2] = x+m;
}
else if (h < 240.0)
{
rgb[0] = m; rgb[1] = x+m; rgb[2] = c+m;
}
else if (h < 300.0)
{
rgb[0] = x+m; rgb[1] = m; rgb[2] = c+m;
}
else if (h < 360.0)
{
rgb[0] = c+m; rgb[1] = m; rgb[2] = x+m;
}
else
{
rgb[0] = m; rgb[1] = m; rgb[2] = m;
}
}
double color_amplitude(value, scale, time)
/* transforms the wave amplitude into a double in [-1,1] to feed into color scheme */
double value, scale;
int time;
{
return(tanh(SLOPE*value/scale)*exp(-((double)time*ATTENUATION)));
}
void color_scheme(scheme, value, scale, time, rgb) /* color scheme */
double value, scale;
int scheme, time;
double rgb[3];
{
double hue, y, r, amplitude;
int intpart;
/* saturation = r, luminosity = y */
switch (scheme) {
case C_LUM:
{
hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = LUMMEAN + amplitude*LUMAMP;
intpart = (int)y;
y -= (double)intpart;
hsl_to_rgb(hue, r, y, rgb);
break;
}
case C_HUE:
{
r = 0.9;
amplitude = color_amplitude(value, scale, time);
y = 0.5;
hue = HUEMEAN + amplitude*HUEAMP;
if (hue < 0.0) hue += 360.0;
if (hue >= 360.0) hue -= 360.0;
hsl_to_rgb(hue, r, y, rgb);
break;
}
}
}
void blank()
{
glClearColor(1.0, 1.0, 1.0, 1.0);
glClear(GL_COLOR_BUFFER_BIT);
}
void save_frame()
{
static int counter = 0;
char *name="wave.", n2[100];
char format[6]=".%05i";
counter++;
// printf (" p2 counter = %d \n",counter);
strcpy(n2, name);
sprintf(strstr(n2,"."), format, counter);
strcat(n2, ".tif");
printf(" saving frame %s \n",n2);
writetiff(n2, "Wave equation in a planar domain", 0, 0,
WINWIDTH, WINHEIGHT, COMPRESSION_LZW);
}
/*********************/
/* some basic math */
/*********************/
double vabs(x) /* absolute value */
double x;
{
double res;
if (x<0.0) res = -x;
else res = x;
return(res);
}
double module2(x, y) /* Euclidean norm */
double x, y;
{
double m;
m = sqrt(x*x + y*y);
return(m);
}
double argument(x, y)
double x, y;
{
double alph;
if (x!=0.0)
{
alph = atan(y/x);
if (x<0.0)
alph += PI;
}
else
{
alph = PID;
if (y<0.0)
alph = PI*1.5;
}
return(alph);
}
/*********************/
/* drawing routines */
/*********************/
/* The billiard boundary is drawn in (x,y) coordinates */
/* However for the grid points, we use integer coordinates (i,j) */
/* GL would allow to always work in (x,y) coordinates but using both */
/* sets of coordinates decreases number of double computations when */
/* drawing the field */
void xy_to_ij(x, y, ij)
/* convert (x,y) position to (i,j) in table representing wave */
double x, y;
int ij[2];
{
double x1, y1;
x1 = (x - XMIN)/(XMAX - XMIN);
y1 = (y - YMIN)/(YMAX - YMIN);
ij[0] = (int)(x1 * (double)NX);
ij[1] = (int)(y1 * (double)NY);
}
void xy_to_pos(x, y, pos)
/* convert (x,y) position to double-valued position in table representing wave */
double x, y, pos[2];
{
double x1, y1;
x1 = (x - XMIN)/(XMAX - XMIN);
y1 = (y - YMIN)/(YMAX - YMIN);
pos[0] = x1 * (double)NX;
pos[1] = y1 * (double)NY;
}
void ij_to_xy(i, j, xy)
/* convert (i,j) position in table representing wave to (x,y) */
int i, j;
double xy[2];
{
double x1, y1;
xy[0] = XMIN + ((double)i)*(XMAX-XMIN)/((double)NX);
xy[1] = YMIN + ((double)j)*(YMAX-YMIN)/((double)NY);
}
int xy_in_billiard(x, y)
/* returns 1 if (x,y) represents a point in the billiard */
double x, y;
{
double l2, r2, omega, c, angle;
int k, condition;
switch (B_DOMAIN) {
case D_RECTANGLE:
{
if ((vabs(x) <LAMBDA)&&(vabs(y) < 1.0)) return(1);
else return(0);
break;
}
case D_ELLIPSE:
{
if (x*x/(LAMBDA*LAMBDA) + y*y < 1.0) return(1);
else return(0);
break;
}
case D_STADIUM:
{
if ((x > -0.5*LAMBDA)&&(x < 0.5*LAMBDA)&&(y > -1.0)&&(y < 1.0)) return(1);
else if (module2(x+0.5*LAMBDA, y) < 1.0) return(1);
else if (module2(x-0.5*LAMBDA, y) < 1.0) return(1);
else return(0);
break;
}
case D_SINAI:
{
if (x*x + y*y > LAMBDA*LAMBDA) return(1);
else return(0);
break;
}
case D_DIAMOND:
{
l2 = LAMBDA*LAMBDA;
r2 = l2 + (LAMBDA-1.0)*(LAMBDA-1.0);
if ((x*x + y*y < 1.0)&&((x-LAMBDA)*(x-LAMBDA) + (y-LAMBDA)*(y-LAMBDA) > r2)
&&((x-LAMBDA)*(x-LAMBDA) + (y+LAMBDA)*(y+LAMBDA) > r2)
&&((x+LAMBDA)*(x+LAMBDA) + (y-LAMBDA)*(y-LAMBDA) > r2)
&&((x+LAMBDA)*(x+LAMBDA) + (y+LAMBDA)*(y+LAMBDA) > r2)) return(1);
else return(0);
break;
}
case D_TRIANGLE:
{
if ((x>-LAMBDA)&&(y>-1.0)&&(LAMBDA*y+x<0.0)) return(1);
else return(0);
break;
}
case D_FLAT:
{
if (y > -LAMBDA) return(1);
else return(0);
break;
}
case D_ANNULUS:
{
l2 = LAMBDA*LAMBDA;
r2 = x*x + y*y;
if ((r2 > l2)&&(r2 < 1.0)) return(1);
else return(0);
}
case D_POLYGON:
{
condition = 1;
omega = DPI/((double)NPOLY);
c = cos(omega*0.5);
for (k=0; k<NPOLY; k++)
{
angle = APOLY*PID + (k+0.5)*omega;
condition = condition*(x*cos(angle) + y*sin(angle) < c);
}
// for (k=0; k<NPOLY; k++) condition = condition*(-x*sin((k+0.5)*omega) + y*cos((k+0.5)*omega) < c);
return(condition);
}
default:
{
printf("Function ij_in_billiard not defined for this billiard \n");
return(0);
}
}
}
void init_wave(x, y, phi, psi, xy_in) 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 */ /* 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]; double x, y, *phi[NX], *psi[NX]; short int * xy_in[NX];
@ -537,233 +181,6 @@ double factor, x, y, *phi[NX], *psi[NX];
int ij_in_billiard(i, j)
/* returns 1 if (i,j) represents a point in the billiard */
int i, j;
{
double xy[2];
ij_to_xy(i, j, xy);
return(xy_in_billiard(xy[0], xy[1]));
}
void draw_billiard() /* draws the billiard boundary */
{
double x0, x, y, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega;
int i;
glColor3f(0.0, 0.0, 0.0);
glLineWidth(5);
glEnable(GL_LINE_SMOOTH);
switch (B_DOMAIN) {
case (D_RECTANGLE):
{
glBegin(GL_LINE_LOOP);
xy_to_pos(LAMBDA, -1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(LAMBDA, 1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-LAMBDA, 1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-LAMBDA, -1.0, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
break;
}
case (D_ELLIPSE):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
/* draw foci */
if (FOCI)
{
glColor3f(0.3, 0.3, 0.3);
x0 = sqrt(LAMBDA*LAMBDA-1.0);
glLineWidth(2);
glEnable(GL_LINE_SMOOTH);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = x0 + r*cos(phi);
y = r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = -x0 + r*cos(phi);
y = r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
}
break;
}
case D_STADIUM:
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = -PID + (double)i*PI/(double)NSEG;
x = 0.5*LAMBDA + cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = PID + (double)i*PI/(double)NSEG;
x = -0.5*LAMBDA + cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
break;
}
case D_SINAI:
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = LAMBDA*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
break;
}
case D_DIAMOND:
{
alpha = atan(1.0 - 1.0/LAMBDA);
dphi = (PID - 2.0*alpha)/(double)NSEG;
r = sqrt(LAMBDA*LAMBDA + (LAMBDA-1.0)*(LAMBDA-1.0));
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = alpha + (double)i*dphi;
x = -LAMBDA + r*cos(phi);
y = -LAMBDA + r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = alpha - PID + (double)i*dphi;
x = -LAMBDA + r*cos(phi);
y = LAMBDA + r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = alpha + PI + (double)i*dphi;
x = LAMBDA + r*cos(phi);
y = LAMBDA + r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
for (i=0; i<=NSEG; i++)
{
phi = alpha + PID + (double)i*dphi;
x = LAMBDA + r*cos(phi);
y = -LAMBDA + r*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd();
break;
}
case (D_TRIANGLE):
{
glBegin(GL_LINE_LOOP);
xy_to_pos(-LAMBDA, -1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(LAMBDA, -1.0, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(-LAMBDA, 1.0, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
break;
}
case (D_FLAT):
{
glBegin(GL_LINE_LOOP);
xy_to_pos(XMIN, -LAMBDA, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(XMAX, -LAMBDA, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
break;
}
case (D_ANNULUS):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = LAMBDA*cos(phi);
y = LAMBDA*sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = cos(phi);
y = sin(phi);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
break;
}
case (D_POLYGON):
{
omega = DPI/((double)NPOLY);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NPOLY; i++)
{
x = cos(i*omega + APOLY*PID);
y = sin(i*omega + APOLY*PID);
xy_to_pos(x, y, pos);
glVertex2d(pos[0], pos[1]);
}
glEnd ();
break;
}
default:
{
printf("Function draw_billiard not defined for this billiard \n");
}
}
}
/*********************/ /*********************/
/* animation part */ /* animation part */
/*********************/ /*********************/
@ -806,7 +223,7 @@ void evolve_wave(phi, psi, xy_in)
int i, j, iplus, iminus, jplus, jminus; int i, j, iplus, iminus, jplus, jminus;
double delta, x, y; double delta, x, y;
#pragma omp parallel for private(i,j,iplus,iminus,delta,x,y) #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
for (i=0; i<NX; i++){ for (i=0; i<NX; i++){
for (j=0; j<NY; j++){ for (j=0; j<NY; j++){
if (xy_in[i][j]){ if (xy_in[i][j]){
@ -824,7 +241,6 @@ void evolve_wave(phi, psi, xy_in)
/* evolve phi */ /* evolve phi */
phi[i][j] = -y + 2*x + courant2*delta - KAPPA*x - GAMMA*(x-y); phi[i][j] = -y + 2*x + courant2*delta - KAPPA*x - GAMMA*(x-y);
// phi[i][j] = -psi[i][j] + 2*x + courant2*delta - GAMMA*x;
/* Old versions of the simulation used this: */ /* Old versions of the simulation used this: */
// phi[i][j] = (-psi[i][j] + 2*phi[i][j] + courant2*delta)*damping; // phi[i][j] = (-psi[i][j] + 2*phi[i][j] + courant2*delta)*damping;
@ -940,10 +356,11 @@ void animation()
} }
if (MOVIE) for (i=0; i<20; i++){ if (MOVIE)
save_frame(); {
s = system("mv wave*.tif tif_wave/"); for (i=0; i<20; i++) save_frame();
} s = system("mv wave*.tif tif_wave/");
}
for (i=0; i<NX; i++) for (i=0; i<NX; i++)
{ {
free(phi[i]); free(phi[i]);