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-06-20 23:31:23 +02:00 committed by GitHub
parent bd6aa073a7
commit 279a6e8801
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 2149 additions and 141 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
drop_billiard.c
View File

@ -42,6 +42,8 @@
// #define YMIN -1.125
// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */
/* Choice of the billiard table */
#define B_DOMAIN 9 /* choice of domain shape */
@ -55,6 +57,8 @@
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_REULEAUX 9 /* Reuleaux and star shapes */
#define D_FLOWER 10 /* Bunimovich flower */
#define D_ALT_REU 11 /* alternating between star and Reuleaux */
// #define LAMBDA 1.0 /* parameter controlling shape of billiard */
#define LAMBDA 1.124950941 /* sin(36°)/sin(31.5°) for 5-star shape with 45° angles */
@ -66,6 +70,8 @@
#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 DRAW_BILLIARD 1 /* set to 1 to draw billiard */
#define DRAW_CONSTRUCTION_LINES 1 /* set to 1 to draw additional construction lines for billiard */
#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */

812
heat.c Normal file
View File

@ -0,0 +1,812 @@
/*********************************************************************************/
/* */
/* Animation of heat 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 heat heat.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_heat */
/* 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 1280 /* number of grid points on x axis */
#define NY 720 /* number of grid points on y axis */
// #define NX 640 /* number of grid points on x axis */
// #define NY 360 /* number of grid points on y axis */
/* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */
/* but will multiply run time by 4 */
#define XMIN -2.0
#define XMAX 2.0 /* x interval */
// #define XMIN -1.5
// #define XMAX 2.5 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define JULIA_SCALE 1.1 /* scaling for Julia sets */
// #define JULIA_SCALE 0.8 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 25 /* 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 D_MENGER 15 /* Menger-Sierpinski carpet */
#define D_JULIA_INT 16 /* interior of Julia set */
/* Billiard tables for heat equation */
#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
#define D_MENGER_HEATED 22 /* Menger gasket with different temperatures */
#define D_MENGER_H_OPEN 23 /* Menger gasket with different temperatures and larger domain */
#define D_MANDELBROT 24 /* Mandelbrot set */
#define D_JULIA 25 /* Julia set */
#define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */
#define LAMBDA 0.7 /* parameter controlling the dimensions of domain */
#define MU 0.1 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 2 /* depth of computation of Menger gasket */
#define MRATIO 5 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
/* 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.00001
#define DT 0.000004
// #define DT 0.000002
// #define DT 0.00000002
// #define DT 0.000000005
#define VISCOSITY 10.0
#define T_OUT 2.0 /* outside temperature */
#define T_IN 0.0 /* inside temperature */
// #define T_OUT 0.0 /* outside temperature */
// #define T_IN 2.0 /* inside temperature */
#define SPEED 0.0 /* speed of drift to the right */
/* Boundary conditions */
#define B_COND 0
#define BC_DIRICHLET 0 /* Dirichlet boundary conditions */
#define BC_PERIODIC 1 /* periodic boundary conditions */
#define BC_ABSORBING 2 /* absorbing boundary conditions (beta version) */
/* Parameters for length and speed of simulation */
#define NSTEPS 4500 /* number of frames of movie */
#define NVID 50 /* number of iterations between images displayed on screen */
// #define NVID 100 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
#define PAUSE 100 /* number of frames after which to pause */
#define PSLEEP 1 /* sleep time during pause */
#define SLEEP1 2 /* 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 */
/* Field representation */
#define FIELD_REP 0
#define F_INTENSITY 0 /* color represents intensity */
#define F_GRADIENT 1 /* color represents norm of gradient */
#define DRAW_FIELD_LINES 1 /* set to 1 to draw field lines */
#define FIELD_LINE_WIDTH 1 /* width of field lines */
#define N_FIELD_LINES 200 /* number of field lines */
#define FIELD_LINE_FACTOR 100 /* factor controlling precision when computing origin of field lines */
/* 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 0 /* set to 1 to adjust color scheme to variance of field */
// #define SLOPE 0.1 /* sensitivity of color on wave amplitude */
#define SLOPE 0.3 /* sensitivity of color on wave amplitude */
#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 280.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP -110.0 /* amplitude of variation of hue for color scheme C_HUE */
// #define HUEMEAN 270.0 /* mean value of hue for color scheme C_HUE */
// #define HUEAMP -130.0 /* amplitude of variation of hue for color scheme C_HUE */
/* Basic math */
#define PI 3.141592654
#define DPI 6.283185307
#define PID 1.570796327
double julia_x = 0.0, julia_y = 0.0; /* parameters for Julia sets */
#include "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_gaussian(x, y, mean, amplitude, scalex, phi, xy_in)
/* initialise field with gaussian at position (x,y) */
double x, y, mean, amplitude, scalex, *phi[NX];
short int * xy_in[NX];
{
int i, j, in;
double xy[2], dist2, module, phase, scale2;
scale2 = scalex*scalex;
printf("Initialising field\n");
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]);
in = xy_in[i][j];
if (in == 1)
{
dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
module = amplitude*exp(-dist2/scale2);
if (module < 1.0e-15) module = 1.0e-15;
phi[i][j] = mean + module/scalex;
} /* boundary temperatures */
else if (in >= 2) phi[i][j] = T_IN*pow(0.75, (double)(in-2));
// else if (in >= 2) phi[i][j] = T_IN*pow(1.0 - 0.5*(double)(in-2), (double)(in-2));
// else if (in >= 2) phi[i][j] = T_IN*(1.0 - (double)(in-2)/((double)MDEPTH))*(1.0 - (double)(in-2)/((double)MDEPTH));
else phi[i][j] = T_OUT;
}
}
void init_julia_set(phi, xy_in)
/* change Julia set boundary condition */
double *phi[NX];
short int * xy_in[NX];
{
int i, j, in;
double xy[2], dist2, module, phase, scale2;
// printf("Changing Julia set\n");
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]);
in = xy_in[i][j];
if (in >= 2) phi[i][j] = T_IN;
}
}
/*********************/
/* animation part */
/*********************/
void compute_gradient(phi, nablax, nablay)
/* compute the gradient of the field */
double *phi[NX], *nablax[NX], *nablay[NX];
{
int i, j, iplus, iminus, jplus, jminus;
double dx;
dx = (XMAX-XMIN)/((double)NX);
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
iplus = i+1; if (iplus == NX) iplus = NX-1;
iminus = i-1; if (iminus == -1) iminus = 0;
jplus = j+1; if (jplus == NX) jplus = NY-1;
jminus = j-1; if (jminus == -1) jminus = 0;
nablax[i][j] = (phi[iplus][j] - phi[iminus][j])/dx;
nablay[i][j] = (phi[i][jplus] - phi[i][jminus])/dx;
}
}
void draw_field_line(x, y, xy_in, nablax, nablay, delta, nsteps)
// void draw_field_line(x, y, nablax, nablay, delta, nsteps)
/* draw a field line of the gradient, starting in (x,y) */
double x, y, *nablax[NX], *nablay[NX], delta;
int nsteps;
short int *xy_in[NX];
{
double x1, y1, x2, y2, pos[2], nabx, naby, norm2, norm;
int i = 0, ij[2], cont = 1;
glColor3f(1.0, 1.0, 1.0);
glLineWidth(FIELD_LINE_WIDTH);
x1 = x;
y1 = y;
// printf("Drawing field line \n");
glEnable(GL_LINE_SMOOTH);
glBegin(GL_LINE_STRIP);
xy_to_pos(x1, y1, pos);
glVertex2d(pos[0], pos[1]);
i = 0;
while ((cont)&&(i < nsteps))
{
xy_to_ij(x1, y1, ij);
if (ij[0] < 0) ij[0] = 0;
if (ij[0] > NX-1) ij[0] = NX-1;
if (ij[1] < 0) ij[1] = 0;
if (ij[1] > NY-1) ij[1] = NY-1;
nabx = nablax[ij[0]][ij[1]];
naby = nablay[ij[0]][ij[1]];
norm2 = nabx*nabx + naby*naby;
if (norm2 > 1.0e-14)
{
/* avoid too large step size */
if (norm2 < 1.0e-9) norm2 = 1.0e-9;
norm = sqrt(norm2);
x1 = x1 + delta*nabx/norm;
y1 = y1 + delta*naby/norm;
}
else
{
cont = 0;
// nablax[ij[0]][ij[1]] = 0.0;
// nablay[ij[0]][ij[1]] = 0.0;
}
if (!xy_in[ij[0]][ij[1]]) cont = 0;
/* stop if the boundary is hit */
// if (xy_in[ij[0]][ij[1]] != 1) cont = 0;
// printf("x1 = %.3lg \t y1 = %.3lg \n", x1, y1);
xy_to_pos(x1, y1, pos);
glVertex2d(pos[0], pos[1]);
i++;
}
glEnd();
}
void draw_wave(phi, xy_in, scale, time)
/* draw the field */
double *phi[NX], scale;
short int *xy_in[NX];
int time;
{
int i, j, iplus, iminus, jplus, jminus, ij[2], counter = 0;
static int first = 1;
double rgb[3], xy[2], x1, y1, x2, y2, dx, value, angle, dangle, intens, deltaintens, sum = 0.0;
double *nablax[NX], *nablay[NX];
static double linex[N_FIELD_LINES*FIELD_LINE_FACTOR], liney[N_FIELD_LINES*FIELD_LINE_FACTOR], distance[N_FIELD_LINES*FIELD_LINE_FACTOR], integral[N_FIELD_LINES*FIELD_LINE_FACTOR + 1];
for (i=0; i<NX; i++)
{
nablax[i] = (double *)malloc(NY*sizeof(double));
nablay[i] = (double *)malloc(NY*sizeof(double));
}
/* compute the gradient */
// if (FIELD_REP > 0)
compute_gradient(phi, nablax, nablay);
/* compute the position of origins of field lines */
if ((first)&&(DRAW_FIELD_LINES))
{
first = 0;
printf("computing linex\n");
x1 = sqrt(3.58);
y1 = 0.0;
linex[0] = x1;
liney[0] = y1;
dangle = DPI/((double)(N_FIELD_LINES*FIELD_LINE_FACTOR));
for (i = 1; i < N_FIELD_LINES*FIELD_LINE_FACTOR; i++)
{
// angle = PID + (double)i*dangle;
angle = (double)i*dangle;
x2 = sqrt(3.58)*cos(angle);
y2 = sqrt(1.18)*sin(angle);
linex[i] = x2;
liney[i] = y2;
distance[i-1] = module2(x2-x1,y2-y1);
x1 = x2;
y1 = y2;
}
distance[N_FIELD_LINES*FIELD_LINE_FACTOR - 1] = module2(x2-sqrt(3.58),y2);
}
dx = (XMAX-XMIN)/((double)NX);
glBegin(GL_QUADS);
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
if (FIELD_REP == F_INTENSITY) value = phi[i][j];
else if (FIELD_REP == F_GRADIENT)
{
value = module2(nablax[i][j], nablay[i][j]);
}
// if ((phi[i][j] - T_IN)*(phi[i][j] - T_OUT) < 0.0)
if (xy_in[i][j] == 1)
{
color_scheme(COLOR_SCHEME, value, scale, time, rgb);
glColor3f(rgb[0], rgb[1], rgb[2]);
}
else glColor3f(0.0, 0.0, 0.0);
glVertex2i(i, j);
glVertex2i(i+1, j);
glVertex2i(i+1, j+1);
glVertex2i(i, j+1);
}
glEnd ();
/* draw a field line */
if (DRAW_FIELD_LINES)
{
/* compute gradient norm along boundary and its integral */
for (i = 0; i < N_FIELD_LINES*FIELD_LINE_FACTOR; i++)
{
xy_to_ij(linex[i], liney[i], ij);
intens = module2(nablax[ij[0]][ij[1]], nablay[ij[0]][ij[1]])*distance[i];
if (i > 0) integral[i] = integral[i-1] + intens;
else integral[i] = intens;
}
deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES);
// deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES + 1.0);
// deltaintens = integral[N_FIELD_LINES*FIELD_LINE_FACTOR-1]/((double)N_FIELD_LINES);
// printf("delta = %.5lg\n", deltaintens);
i = 0;
// draw_field_line(linex[0], liney[0], nablax, nablay, 0.00002, 100000);
draw_field_line(linex[0], liney[0], xy_in, nablax, nablay, 0.00002, 100000);
for (j = 1; j < N_FIELD_LINES+1; j++)
{
while ((integral[i] <= j*deltaintens)&&(i < N_FIELD_LINES*FIELD_LINE_FACTOR)) i++;
// draw_field_line(linex[i], liney[i], nablax, nablay, 0.00002, 100000);
draw_field_line(linex[i], liney[i], xy_in, nablax, nablay, 0.00002, 100000);
counter++;
}
printf("%i lines\n", counter);
}
for (i=0; i<NX; i++)
{
free(nablax[i]);
free(nablay[i]);
}
}
void evolve_wave(phi, xy_in)
/* time step of field evolution */
double *phi[NX]; short int *xy_in[NX];
{
int i, j, iplus, iminus, jplus, jminus;
double delta1, delta2, x, y, *newphi[NX];;
for (i=0; i<NX; i++) newphi[i] = (double *)malloc(NY*sizeof(double));
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta1,delta2,x,y)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j] == 1){
/* discretized Laplacian depending on boundary conditions */
if ((B_COND == BC_DIRICHLET)||(B_COND == BC_ABSORBING))
{
iplus = (i+1); if (iplus == NX) iplus = NX-1;
iminus = (i-1); if (iminus == -1) iminus = 0;
jplus = (j+1); if (jplus == NY) jplus = NY-1;
jminus = (j-1); if (jminus == -1) jminus = 0;
}
else if (B_COND == BC_PERIODIC)
{
iplus = (i+1) % NX;
iminus = (i-1) % NX;
if (iminus < 0) iminus += NX;
jplus = (j+1) % NY;
jminus = (j-1) % NY;
if (jminus < 0) jminus += NY;
}
delta1 = phi[iplus][j] + phi[iminus][j] + phi[i][jplus] + phi[i][jminus] - 4.0*phi[i][j];
x = phi[i][j];
/* evolve phi */
if (B_COND != BC_ABSORBING)
{
newphi[i][j] = x + intstep*(delta1 - SPEED*(phi[iplus][j] - phi[i][j]));
}
else /* case of absorbing b.c. - this is only an approximation of correct way of implementing */
{
/* in the bulk */
if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
{
newphi[i][j] = x - intstep*delta2;
}
/* right border */
else if (i==NX-1)
{
newphi[i][j] = x - intstep1*(x - phi[i-1][j]);
}
/* upper border */
else if (j==NY-1)
{
newphi[i][j] = x - intstep1*(x - phi[i][j-1]);
}
/* left border */
else if (i==0)
{
newphi[i][j] = x - intstep1*(x - phi[1][j]);
}
/* lower border */
else if (j==0)
{
newphi[i][j] = x - intstep1*(x - phi[i][1]);
}
}
if (FLOOR)
{
if (newphi[i][j] > VMAX) phi[i][j] = VMAX;
if (newphi[i][j] < -VMAX) phi[i][j] = -VMAX;
}
}
}
}
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j] == 1) phi[i][j] = newphi[i][j];
}
}
for (i=0; i<NX; i++)
{
free(newphi[i]);
}
// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
}
double compute_variance(phi, xy_in)
/* compute the variance (total probability) of the field */
double *phi[NX]; short int * xy_in[NX];
{
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];
}
}
if (n==0) n=1;
return(variance/(double)n);
}
void renormalise_field(phi, xy_in, variance)
/* renormalise variance of field */
double *phi[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;
}
}
}
void print_level(level)
int level;
{
double pos[2];
char message[50];
glColor3f(1.0, 1.0, 1.0);
sprintf(message, "Level %i", level);
xy_to_pos(XMIN + 0.1, YMAX - 0.2, pos);
write_text(pos[0], pos[1], message);
}
void print_Julia_parameters()
{
double pos[2];
char message[50];
glColor3f(1.0, 1.0, 1.0);
if (julia_y >= 0.0) sprintf(message, "c = %.5f + %.5f i", julia_x, julia_y);
else sprintf(message, "c = %.5f %.5f i", julia_x, julia_y);
xy_to_pos(XMIN + 0.1, YMAX - 0.2, pos);
write_text(pos[0], pos[1], message);
}
void set_Julia_parameters(time, phi, xy_in)
int time;
double *phi[NX];
short int *xy_in[NX];
{
double jangle, cosj, sinj, radius = 0.15;
jangle = (double)time*DPI/(double)NSTEPS;
// jangle = (double)time*0.001;
// jangle = (double)time*0.0001;
cosj = cos(jangle);
sinj = sin(jangle);
julia_x = -0.9 + radius*cosj;
julia_y = radius*sinj;
init_julia_set(phi, xy_in);
printf("Julia set parameters : i = %i, angle = %.5lg, cx = %.5lg, cy = %.5lg \n", time, jangle, julia_x, julia_y);
}
void set_Julia_parameters_cardioid(time, phi, xy_in)
int time;
double *phi[NX];
short int *xy_in[NX];
{
double jangle, cosj, sinj, yshift;
jangle = pow(1.05 + (double)time*0.00003, 0.333);
yshift = 0.02*sin((double)time*PID*0.002);
// jangle = pow(1.0 + (double)time*0.00003, 0.333);
// jangle = pow(0.05 + (double)time*0.00003, 0.333);
// jangle = pow(0.1 + (double)time*0.00001, 0.333);
// yshift = 0.04*sin((double)time*PID*0.002);
cosj = cos(jangle);
sinj = sin(jangle);
julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj);
julia_y = 0.5*sinj*(1.0-cosj) + yshift;
/* need to decrease 0.05 for i > 2000 */
// julia_x = 0.5*(cosj*(1.0 - 0.5*cosj) + 0.5*sinj*sinj);
// julia_y = 0.5*sinj*(1.0-cosj);
init_julia_set(phi, xy_in);
printf("Julia set parameters : i = %i, angle = %.5lg, cx = %.5lg, cy = %.5lg \n", time, jangle, julia_x, julia_y);
}
void animation()
{
double time, scale, dx, var, jangle, cosj, sinj;
double *phi[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));
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
dx = (XMAX-XMIN)/((double)NX);
intstep = DT/(dx*dx*VISCOSITY);
intstep1 = DT/(dx*VISCOSITY);
// julia_x = 0.1;
// julia_y = 0.6;
set_Julia_parameters(0, phi, xy_in);
printf("Integration step %.3lg\n", intstep);
/* initialize wave wave function */
init_gaussian(-1.0, 0.0, 0.1, 0.0, 0.01, phi, xy_in);
// init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in)
if (SCALE)
{
var = compute_variance(phi, xy_in);
scale = sqrt(1.0 + var);
renormalise_field(phi, xy_in, var);
}
blank();
glColor3f(0.0, 0.0, 0.0);
glutSwapBuffers();
draw_wave(phi, xy_in, 1.0, 0);
draw_billiard();
print_Julia_parameters(i);
// print_level(MDEPTH);
glutSwapBuffers();
sleep(SLEEP1);
if (MOVIE) for (i=0; i<SLEEP1*25; i++) save_frame();
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, xy_in);
scale = sqrt(1.0 + var);
// printf("Norm: %5lg\t Scaling factor: %5lg\n", var, scale);
renormalise_field(phi, xy_in, var);
}
else scale = 1.0;
draw_wave(phi, xy_in, scale, i);
for (j=0; j<NVID; j++) evolve_wave(phi, xy_in);
draw_billiard();
// print_level(MDEPTH);
print_Julia_parameters(i);
glutSwapBuffers();
/* modify Julia set */
set_Julia_parameters(i, phi, xy_in);
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_heat/");
}
}
}
if (MOVIE)
{
for (i=0; i<20; i++) save_frame();
s = system("mv wave*.tif tif_heat/");
}
for (i=0; i<NX; i++)
{
free(phi[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("Heat equation in a planar domain");
init();
glutDisplayFunc(display);
glutMainLoop();
return 0;
}

113
particle_billiard.c
View File

@ -39,6 +39,9 @@
#define XMAX 2.0 /* x interval */
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */
// #define XMIN -1.8
// #define XMAX 1.8 /* x interval */
// #define YMIN -0.91
@ -57,31 +60,39 @@
#define D_ANNULUS 7 /* annulus */
#define D_POLYGON 8 /* polygon */
#define D_REULEAUX 9 /* Reuleaux and star shapes */
#define D_FLOWER 10 /* Bunimovich flower */
#define D_ALT_REU 11 /* alternating between star and Reuleaux */
// #define LAMBDA 5.0 /* parameter controlling shape of billiard */
#define LAMBDA -1.949855824 /* 7-sided Reuleaux triangle */
#define LAMBDA -3.346065215 /* sin(60°)/sin(15°) for Reuleaux-type triangle with 90° angles */
// #define LAMBDA 3.0 /* parameter controlling shape of billiard */
// #define LAMBDA 0.6 /* parameter controlling shape of billiard */
// #define LAMBDA 0.4175295 /* sin(20°)/sin(55°) for 9-star shape with 30° angles */
// #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 MU 0.1 /* second parameter controlling shape of billiard */
#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 NPOLY 6 /* number of sides of polygon */
#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */
#define DRAW_BILLIARD 1 /* set to 1 to draw billiard */
#define DRAW_CONSTRUCTION_LINES 1 /* set to 1 to draw additional construction lines for billiard */
#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */
#define DEBUG 0 /* draw trajectories, for debugging purposes */
/* Simulation parameters */
#define NPART 10000 /* number of particles */
#define NPART 20000 /* number of particles */
#define NPARTMAX 100000 /* maximal number of particles after resampling */
#define LMAX 0.01 /* minimal segment length triggering resampling */
#define DMIN 0.02 /* minimal distance to boundary for triggering resampling */
#define CYCLE 1 /* set to 1 for closed curve (start in all directions) */
#define NSTEPS 4000 /* number of frames of movie */
#define TIME 125 /* time between movie frames, for fluidity of real-time simulation */
#define DPHI 0.00001 /* integration step */
#define NSTEPS 6000 /* number of frames of movie */
#define TIME 1000 /* time between movie frames, for fluidity of real-time simulation */
// #define DPHI 0.000005 /* integration step */
#define DPHI 0.00002 /* integration step */
#define NVID 150 /* number of iterations between images displayed on screen */
/* Decreasing TIME accelerates the animation and the movie */
@ -92,18 +103,19 @@
/* 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 NCOLORS -10 /* number of colors */
#define COLORSHIFT 200 /* hue of initial color */
#define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */
#define NSEG 100 /* number of segments of boundary */
#define LENGTH 0.04 /* length of velocity vectors */
#define BILLIARD_WIDTH 3 /* width of billiard */
#define 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 PAINT_INT 0 /* set to 1 to paint interior in other color (for polygon/Reuleaux) */
#define PAUSE 1000 /* number of frames after which to pause */
@ -230,6 +242,7 @@ double *configs[NPARTMAX];
glutSwapBuffers();
blank();
if (PAINT_INT) paint_billiard_interior();
glLineWidth(PARTICLE_WIDTH);
@ -262,40 +275,40 @@ double *configs[NPARTMAX];
glColor3f(rgb[0], rgb[1], rgb[2]);
glBegin(GL_LINE_STRIP);
glVertex2d(x1, y1);
glVertex2d(x2, y2);
glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*y1);
glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*y2);
glEnd ();
/* taking care of boundary conditions - only needed for periodic boundary conditions */
if (x2 > XMAX)
if (SCALING_FACTOR*x2 > XMAX)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1+XMIN-XMAX, y1);
glVertex2d(x2+XMIN-XMAX, y2);
glVertex2d(SCALING_FACTOR*(x1+XMIN-XMAX), SCALING_FACTOR*y1);
glVertex2d(SCALING_FACTOR*(x2+XMIN-XMAX), SCALING_FACTOR*y2);
glEnd ();
}
if (x2 < XMIN)
if (SCALING_FACTOR*x2 < XMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1-XMIN+XMAX, y1);
glVertex2d(x2-XMIN+XMAX, y2);
glVertex2d(SCALING_FACTOR*(x1-XMIN+XMAX), SCALING_FACTOR*y1);
glVertex2d(SCALING_FACTOR*(x2-XMIN+XMAX), SCALING_FACTOR*y2);
glEnd ();
}
if (y2 > YMAX)
if (SCALING_FACTOR*y2 > YMAX)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1, y1+YMIN-YMAX);
glVertex2d(x2, y2+YMIN-YMAX);
glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMIN-YMAX));
glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMIN-YMAX));
glEnd ();
}
if (y2 < YMIN)
if (SCALING_FACTOR*y2 < YMIN)
{
glBegin(GL_LINE_STRIP);
glVertex2d(x1, y1+YMAX-YMIN);
glVertex2d(x2, y2+YMAX-YMIN);
glVertex2d(SCALING_FACTOR*x1, SCALING_FACTOR*(y1+YMAX-YMIN));
glVertex2d(SCALING_FACTOR*x2, SCALING_FACTOR*(y2+YMAX-YMIN));
glEnd ();
}
}
@ -305,15 +318,15 @@ double *configs[NPARTMAX];
{
glLineWidth(1.0);
glBegin(GL_LINES);
glVertex2d(configs[i][4], configs[i][5]);
glVertex2d(configs[i][6], configs[i][7]);
glVertex2d(SCALING_FACTOR*configs[i][4], SCALING_FACTOR*configs[i][5]);
glVertex2d(SCALING_FACTOR*configs[i][6], SCALING_FACTOR*configs[i][7]);
glEnd ();
glLineWidth(3.0);
}
if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3];
}
draw_billiard();
if (DRAW_BILLIARD) draw_billiard();
}
@ -355,9 +368,9 @@ double *configs[NPARTMAX];
void animation()
{
double time, dt, alpha;
double time, dt, alpha, r;
double *configs[NPARTMAX];
int i, j, resamp = 1, s;
int i, j, resamp = 1, s, i1, i2;
int *color, *newcolor, *active;
/* Since NPARTMAX can be big, it seemed wiser to use some memory allocation here */
@ -368,19 +381,21 @@ void animation()
configs[i] = (double *)malloc(8*sizeof(double));
/* initialize system by putting particles in a given point with a range of velocities */
alpha = 3.0*PI/7.0;
init_sym_drop_config(-0.99, 0.0, -alpha, alpha, configs);
r = cos(PI/(double)NPOLY)/cos(DPI/(double)NPOLY);
// init_drop_config(0.0, 0.0, -0.2, 0.2, configs);
// init_drop_config(0.5, 0.5, -1.0, 1.0, configs);
// init_sym_drop_config(-1.0, 0.5, -PID, PID, configs);
// init_drop_config(-0.999, 0.0, -alpha, alpha, configs);
// init_drop_config(0.0, 0.5, 0.0, DPI, configs);
// other possible initial conditions :
// 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);
init_line_config(0.0, -0.3, 0.0, 0.3, 0.0, configs);
blank();
glColor3f(0.0, 0.0, 0.0);
draw_billiard();
if (DRAW_BILLIARD) draw_billiard();
glutSwapBuffers();
@ -391,6 +406,24 @@ void animation()
newcolor[i] = 0;
active[i] = 1;
}
if (FLOWER_COLOR) /* adapt color scheme to flower configuration (beta implementation) */
{
// i1 = (int)((double)NPART*0.2538); /* the 0.27 is just a trial-and-error guess, to be improved */
// i1 = (int)((double)NPART*0.1971); /* the 0.27 is just a trial-and-error guess, to be improved */
i1 = (int)((double)NPART*0.3015); /* the 0.27 is just a trial-and-error guess, to be improved */
i2 = NPART-i1;
for (i=i1; i<i2; i++)
{
color[i] += NCOLORS/3;
newcolor[i] = NCOLORS/3;
}
for (i=i2; i<NPART; i++)
{
color[i] += 2*NCOLORS/3;
newcolor[i] = 2*NCOLORS/3;
}
}
sleep(SLEEP1);
@ -399,7 +432,7 @@ void animation()
graph_movie(TIME, newcolor, configs, active);
draw_config(newcolor, configs, active);
draw_billiard();
if (DRAW_BILLIARD) draw_billiard();
for (j=0; j<NPARTMAX; j++) color[j] = newcolor[j];

98
schrodinger.c
View File

@ -53,6 +53,8 @@
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define JULIA_SCALE 1.0 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 3 /* choice of domain shape */
@ -70,10 +72,26 @@
#define D_GRATING 10 /* diffraction grating */
#define D_EHRENFEST 11 /* Ehrenfest urn type geometry */
#define D_MENGER 15 /* Menger-Sierpinski carpet */
#define D_JULIA_INT 16 /* interior of Julia set */
/* Billiard tables for heat equation */
#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
#define D_MENGER_HEATED 22 /* Menger gasket with different temperatures */
#define D_MENGER_H_OPEN 23 /* Menger gasket with different temperatures and larger domain */
#define D_MANDELBROT 24 /* Mandelbrot set */
#define D_JULIA 25 /* Julia set */
#define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */
#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */
#define MU 0.05 /* parameter controlling the dimensions of domain */
#define NPOLY 6 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 2 /* depth of computation of Menger gasket */
#define MRATIO 5 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
/* You can add more billiard tables by adapting the functions */
@ -81,7 +99,8 @@
/* Physical patameters of wave equation */
#define DT 0.00000002
#define DT 0.00000005
// #define DT 0.00000002
// #define DT 0.000000005
#define HBAR 1.0
@ -96,7 +115,7 @@
/* 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 NVID 250 /* 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 */
@ -145,6 +164,8 @@
#define DPI 6.283185307
#define PID 1.570796327
double julia_x = 0.0, julia_y = 0.0; /* parameters for Julia sets */
#include "sub_wave.c"
double courant2; /* Courant parameter squared */
@ -154,7 +175,6 @@ 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 */
@ -252,10 +272,18 @@ int time;
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];
double *phi[NX], *psi[NX];
short int *xy_in[NX];
{
int i, j, iplus, iminus, jplus, jminus;
double delta1, delta2, x, y;
double delta1, delta2, x, y, *newphi[NX], *newpsi[NX];
for (i=0; i<NX; i++)
{
newphi[i] = (double *)malloc(NY*sizeof(double));
newpsi[i] = (double *)malloc(NY*sizeof(double));
}
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta1,delta2,x,y)
for (i=0; i<NX; i++){
@ -288,54 +316,69 @@ void evolve_wave(phi, psi, xy_in)
/* evolve phi and psi */
if (B_COND != BC_ABSORBING)
{
phi[i][j] = x - intstep*delta2;
psi[i][j] = y + intstep*delta1;
newphi[i][j] = x - intstep*delta2;
newpsi[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;
newphi[i][j] = x - intstep*delta2;
newpsi[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]);
newphi[i][j] = x - intstep1*(y - psi[i-1][j]);
newpsi[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]);
newphi[i][j] = x - intstep1*(y - psi[i][j-1]);
newpsi[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]);
newphi[i][j] = x - intstep1*(y - psi[1][j]);
newpsi[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]);
newphi[i][j] = x - intstep1*(y - psi[i][1]);
newpsi[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;
if (newphi[i][j] > VMAX) newphi[i][j] = VMAX;
if (newphi[i][j] < -VMAX) newphi[i][j] = -VMAX;
if (newpsi[i][j] > VMAX) newpsi[i][j] = VMAX;
if (newpsi[i][j] < -VMAX) newpsi[i][j] = -VMAX;
}
}
}
}
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j] == 1) phi[i][j] = newphi[i][j];
if (xy_in[i][j] == 1) psi[i][j] = newpsi[i][j];
}
}
for (i=0; i<NX; i++)
{
free(newphi[i]);
free(newpsi[i]);
}
// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
}
@ -409,12 +452,15 @@ void animation()
// 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);
@ -439,12 +485,14 @@ void animation()
else scale = 1.0;
draw_wave(phi, psi, xy_in, scale, i);
// printf("Wave drawn\n");
for (j=0; j<NVID; j++) evolve_wave(phi, psi, xy_in);
draw_billiard();
glutSwapBuffers();
glutSwapBuffers();
if (MOVIE)
{

726
sub_part_billiard.c
View File

@ -239,10 +239,45 @@ void write_text( double x, double y, char *st)
}
void paint_billiard_interior() /* points billiard interior, for use before draw_conf */
void compute_flower_parameters(omega, co, so, axis1, axis2, phimax)
/* compute parameters needed for the flower billiard in terms of LAMBDA and NPOLY */
double *omega, *co, *so, *axis1, *axis2, *phimax;
{
double omega2, co2, so2, r, a, gamma, axissquare1;
/* various angles */
*omega = DPI/((double)NPOLY);
omega2 = PI/((double)NPOLY);
co2 = cos(omega2);
so2 = sin(omega2);
*co = cos(*omega);
*so = sin(*omega);
// *co = co2*co2 - so2*so2;
// *so = 2.0*co2*so2;
/* distance of edge of ellipse to the origin */
r = LAMBDA*co2/(*co);
a = (r*co2 - *co)*(r*co2 - *co);
gamma = 0.5*r*r - r*co2*(*co) + 0.5*cos(2.0*(*omega));
axissquare1 = gamma + sqrt(gamma*gamma + a*(*so)*(*so));
/* semi-minor axis */
*axis1 = sqrt(axissquare1);
/* semi-major axis */
*axis2 = sqrt(axissquare1 + (*so)*(*so));
/* max angle in ellipse parametrization */
*phimax = asin(r*so2/(*axis2));
}
void paint_billiard_interior() /* paints billiard interior, for use before draw_conf */
{
double x0, x, y, phi, r = 0.01, alpha, dphi, omega;
int i, k, c;
double x0, x, y, phi, r = 0.01, alpha, dphi, omega, beta2, x2, s, x1, y1, angle, co, so, axis1, axis2, phimax;
int i, j, k, c;
glLineWidth(4);
@ -273,6 +308,64 @@ void paint_billiard_interior() /* points billiard interior, for use before
}
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));
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++)
{
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;
}
case (D_FLOWER):
{
compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
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++)
{
for (j=0; j<NSEG; j++)
{
s = 2.0*(((double)j/(double)NSEG)-0.5)*phimax;
x1 = co + axis1*cos(s);
y1 = axis2*sin(s);
angle = i*omega + APOLY*PID;
x = cos(angle)*x1 - sin(angle)*y1;
y = sin(angle)*x1 + cos(angle)*y1;
glVertex2d(SCALING_FACTOR*x, SCALING_FACTOR*y);
}
}
glEnd();
}
break;
}
default:
{
@ -282,7 +375,8 @@ void paint_billiard_interior() /* points billiard interior, for use before
void draw_billiard() /* draws the billiard boundary */
{
double x0, x, y, phi, r = 0.01, alpha, dphi, omega, x1, y1, x2, beta2, angle, s;
double x0, x, y, phi, r = 0.01, alpha, dphi, omega, x1, y1, x2, beta2, angle, s, x2plus, x2minus;
double omega2, co, so, axis1, axis2, phimax;
int i, j, k, c;
if (PAINT_INT) glColor3f(0.5, 0.5, 0.5);
@ -366,6 +460,20 @@ void draw_billiard() /* draws the billiard boundary */
glVertex2d(x, y);
}
glEnd();
if (DRAW_CONSTRUCTION_LINES)
{
glColor3f(0.5, 0.5, 0.5);
glBegin(GL_LINE_STRIP);
glVertex2d(-LAMBDA, -1.0);
glVertex2d(-LAMBDA, 1.0);
glEnd();
glBegin(GL_LINE_STRIP);
glVertex2d(LAMBDA, -1.0);
glVertex2d(LAMBDA, 1.0);
glEnd();
}
break;
}
case D_SINAI:
@ -580,6 +688,86 @@ void draw_billiard() /* draws the billiard boundary */
glEnd ();
break;
}
case (D_FLOWER):
{
compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
/* draw inner polygon and radial lines */
if (DRAW_CONSTRUCTION_LINES)
{
glColor3f(0.5, 0.5, 0.5);
glBegin(GL_LINE_LOOP);
for (i=0; i<=NPOLY; i++)
{
x = cos(i*omega + APOLY*PID);
y = sin(i*omega + APOLY*PID);
glVertex2d(SCALING_FACTOR*x, SCALING_FACTOR*y);
}
glEnd ();
r = LAMBDA*cos(0.5*omega)/co;
for (i=0; i<=NPOLY; i++)
{
glBegin(GL_LINE_STRIP);
glVertex2d(0.0, 0.0);
x = r*cos(((double)i + 0.5)*omega + APOLY*PID);
y = r*sin(((double)i + 0.5)*omega + APOLY*PID);
glVertex2d(SCALING_FACTOR*x, SCALING_FACTOR*y);
glEnd ();
}
}
/* draw billiard boundary */
if (!PAINT_INT)
{
if (BLACK) glColor3f(1.0, 1.0, 1.0);
else glColor3f(0.0, 0.0, 0.0);
}
glBegin(GL_LINE_STRIP);
for (i=0; i<=NPOLY; i++)
// for (i=0; i<=1; i++)
{
for (j=0; j<NSEG; j++)
{
// s = 2.0*(((double)j/(double)NSEG)-0.5)*PI;
s = 2.0*(((double)j/(double)NSEG)-0.5)*phimax;
x1 = co + axis1*cos(s);
y1 = axis2*sin(s);
angle = i*omega + APOLY*PID;
x = cos(angle)*x1 - sin(angle)*y1;
y = sin(angle)*x1 + cos(angle)*y1;
glVertex2d(SCALING_FACTOR*x, SCALING_FACTOR*y);
}
}
glEnd ();
break;
}
case (D_ALT_REU):
{
omega = DPI/((double)NPOLY);
beta2 = asin(sin(omega*0.5)/LAMBDA);
x2plus = cos(omega*0.5) + sqrt(LAMBDA*LAMBDA - sin(omega*0.5)*sin(omega*0.5));
x2minus = 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;
if (i%2==0) x1 = x2plus - LAMBDA*cos(s);
else x1 = x2minus + 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:
{
printf("Function draw_billiard not defined for this billiard \n");
@ -1688,57 +1876,38 @@ int color[NPARTMAX];
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);
x = pos[0]*ca + pos[1]*sa;
y = -pos[0]*sa + pos[1]*ca;
a = (x-x2)*cos(theta) + y*sin(theta);
b = (x-x2)*(x-x2) + y*y - LAMBDA*LAMBDA;
if (a*a - b > margin)
{
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)
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);
if ((t > margin)&&(vabs(yi) <= sin(omega2)))
{
// 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;
cval[nt] = k;
tval[nt] = t;
/* rotate back */
x1[nt] = xi*ca - yi*sa;
y1[nt] = xi*sa + yi*ca;
/* 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++;
}
tempconf[nt][0] = ((double)k + 0.5)*beta + asin(yi/LAMBDA);
tempconf[nt][1] = PID - asin(yi/LAMBDA) - theta;
nt++;
}
}
}
// printf("nt = %i\n", nt);
/* find earliest intersection */
tmin = tval[0];
@ -1783,6 +1952,390 @@ int color[NPARTMAX];
return(c);
}
/****************************************************************************************/
/* Bunimovich flower billiard */
/****************************************************************************************/
int pos_flower(conf, pos, alpha)
/* determine position on boundary of polygon */
/* conf[0] is arclength on boundary, it belongs to [0,2*NPOLY*phimax) */
double conf[2], pos[2], *alpha;
{
double s, theta, omega, co, so, axis1, axis2, phimax, s1, x, y, angle;
int c;
s = conf[0];
theta = conf[1];
compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
c = (int)(s/(2.0*phimax)); /* side of shape */
s1 = s - (((double)c)*2.0 + 1.0)*phimax;
x = co + axis1*cos(s1);
y = axis2*sin(s1);
angle = ((double)c)*omega + PID*APOLY;
// angle = 2.0*((double)c)*omega + PID*APOLY;
pos[0] = x*cos(angle) - y*sin(angle);
pos[1] = x*sin(angle) + y*cos(angle);
*alpha = argument(-axis1*sin(s1), axis2*cos(s1)) + theta + angle;
// printf("alpha = %.5lg\t", *alpha);
return(c);
}
int vflower_xy(config, alpha, pos)
/* determine initial configuration for start at point pos = (x,y) */
double config[8], alpha, pos[2];
{
double s, theta, omega, omega2, s1, rangle, x, y, x1, y1, xi, yi, t;
double ca, sa, aa, bb, cc, margin = 1.0e-14, tmin;
double co, so, co2, so2, ct, st, phimax, phi, axis1, axis2;
int k, c, intb=1, intc, i, nt = 0, ntmin, sign;
compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
for (k=0; k<NPOLY; k++) if (intb)
{
/* rotate position so that kth side is vertical */
// rangle = (double)(2*k)*omega + APOLY*PID;
rangle = (double)k*omega + APOLY*PID;
theta = alpha - rangle;
ca = cos(rangle);
sa = sin(rangle);
ct = cos(theta);
st = sin(theta);
x = pos[0]*ca + pos[1]*sa;
y = -pos[0]*sa + pos[1]*ca;
/* find intersection with elliptical arc */
aa = ct*ct/(axis1*axis1) + st*st/(axis2*axis2);
bb = (x-co)*ct/(axis1*axis1) + y*st/(axis2*axis2);
cc = (x-co)*(x-co)/(axis1*axis1) + y*y/(axis2*axis2) - 1.0;
// if (bb*bb - aa*cc > margin)
if (bb*bb - aa*cc >= 0.0)
{
t = (-bb + sqrt(bb*bb - aa*cc))/aa;
xi = x + t*cos(theta);
yi = y + t*sin(theta);
if (yi >= 0.0) phi = argument((xi - co)/axis1, yi/axis2);
else phi = -argument((xi - co)/axis1, -yi/axis2);
if ((t > margin)&&((vabs(phi) <= phimax)||(vabs(phi-DPI) <= phimax)))
{
intb = 0;
c = k;
/* rotate back */
x1 = xi*ca - yi*sa;
y1 = xi*sa + yi*ca;
config[0] = (double)(2*k + 1)*phimax + phi;
config[1] = argument(-axis1*sin(phi), axis2*cos(phi)) - theta;
}
}
}
// if (nt == 0) 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-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 old_vflower_xy(config, alpha, pos)
/* determine initial configuration for start at point pos = (x,y) */
double config[8], alpha, pos[2];
{
double s, theta, omega, omega2, s1, rangle, x, y, x1[2*NPOLY], y1[2*NPOLY], xi, yi, t;
double ca, sa, aa, bb, cc, margin = 1.0e-14, tmin, tval[2*NPOLY], tempconf[2*NPOLY][2];
double co, so, co2, so2, ct, st, phimax, phi, axis1, axis2;
int k, c, intb=1, intc, i, nt = 0, cval[2*NPOLY], ntmin, sign;
compute_flower_parameters(&omega, &co, &so, &axis1, &axis2, &phimax);
for (k=0; k<NPOLY; k++)
{
/* rotate position so that kth side is vertical */
// rangle = (double)(2*k)*omega + APOLY*PID;
rangle = (double)k*omega + APOLY*PID;
theta = alpha - rangle;
ca = cos(rangle);
sa = sin(rangle);
ct = cos(theta);
st = sin(theta);
x = pos[0]*ca + pos[1]*sa;
y = -pos[0]*sa + pos[1]*ca;
/* find intersection with elliptical arc */
aa = ct*ct/(axis1*axis1) + st*st/(axis2*axis2);
bb = (x-co)*ct/(axis1*axis1) + y*st/(axis2*axis2);
cc = (x-co)*(x-co)/(axis1*axis1) + y*y/(axis2*axis2) - 1.0;
// if (bb*bb - aa*cc > margin)
if (bb*bb - aa*cc >= 0.0)
{
t = (-bb + sqrt(bb*bb - aa*cc))/aa;
xi = x + t*cos(theta);
yi = y + t*sin(theta);
if (yi >= 0.0) phi = argument((xi - co)/axis1, yi/axis2);
else phi = -argument((xi - co)/axis1, -yi/axis2);
// phi = argument((xi - co)/axis1, yi/axis2);
// if (phi > PI) phi += -DPI;
if ((t > margin)&&((vabs(phi) <= phimax)||(vabs(phi-DPI) <= phimax)))
// if (((vabs(phi) <= phimax)||(vabs(phi-DPI) <= phimax)))
// if ((t > margin))
{
cval[nt] = k;
// cval[nt] = 2*k;
tval[nt] = t;
/* rotate back */
x1[nt] = xi*ca - yi*sa;
y1[nt] = xi*sa + yi*ca;
tempconf[nt][0] = (double)(2*k + 1)*phimax + phi;
tempconf[nt][1] = argument(-axis1*sin(phi), axis2*cos(phi)) - theta;
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];
if (nt == 0) 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 vflower(config)
/* determine initial configuration when starting from boundary */
double config[8];
{
double pos[2], alpha;
int c;
c = pos_flower(config, pos, &alpha);
vflower_xy(config, alpha, pos);
return(c);
}
/****************************************************************************************/
/* Alternating between Reuleaux-type and star-shaped billiard */
/****************************************************************************************/
int pos_alt_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, x2plus, x2minus, 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;
x2plus = cos(omega2) + sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
x2minus = cos(omega2) - sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
if (c%2 == 0) x = x2plus - LAMBDA*cos(s1 - beta2);
else x = x2minus + LAMBDA*cos(s1 - beta2);
if (c%2 == 0) y = LAMBDA*sin(s1 - beta2);
else y = -LAMBDA*sin(s1 - beta2);
/* test, to be removed */
// x = x2plus - 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 valt_reuleaux_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, x2plus, x2minus, arcsine;
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);
x2plus = cos(omega2) + sqrt(LAMBDA*LAMBDA - sin(omega2)*sin(omega2));
x2minus = 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;
ca = cos(rangle);
sa = sin(rangle);
x = pos[0]*ca + pos[1]*sa;
y = -pos[0]*sa + pos[1]*ca;
if (k%2==0)
{
a = (x-x2plus)*cos(theta) + y*sin(theta);
b = (x-x2plus)*(x-x2plus) + y*y - LAMBDA*LAMBDA;
}
else
{
a = (x-x2minus)*cos(theta) + y*sin(theta);
b = (x-x2minus)*(x-x2minus) + y*y - LAMBDA*LAMBDA;
}
if (a*a - b > margin)
{
if (k%2==0) t = -a - sqrt(a*a - b);
else t = -a + sqrt(a*a - b);
xi = x + t*cos(theta);
yi = y + t*sin(theta);
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;
if (k%2==0) arcsine = asin(yi/LAMBDA);
else arcsine = -asin(yi/LAMBDA);
tempconf[nt][0] = ((double)k + 0.5)*beta + arcsine;
tempconf[nt][1] = PID - arcsine - theta;
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 valt_reuleaux(config)
/* determine initial configuration when starting from boundary */
double config[8];
{
double pos[2], alpha;
int c;
c = pos_alt_reuleaux(config, pos, &alpha);
valt_reuleaux_xy(config, alpha, pos);
return(c);
}
/****************************************************************************************/
/* general billiard */
@ -1833,6 +2386,16 @@ int color[NPARTMAX];
return(pos_reuleaux(conf, pos, &alpha));
break;
}
case (D_FLOWER):
{
return(pos_flower(conf, pos, &alpha));
break;
}
case (D_ALT_REU):
{
return(pos_alt_reuleaux(conf, pos, &alpha));
break;
}
default:
{
printf("Function pos_billiard not defined for this billiard \n");
@ -1887,6 +2450,16 @@ int color[NPARTMAX];
return(vreuleaux_xy(config, alpha, pos));
break;
}
case (D_FLOWER):
{
return(vflower_xy(config, alpha, pos));
break;
}
case (D_ALT_REU):
{
return(valt_reuleaux_xy(config, alpha, pos));
break;
}
default:
{
printf("Function vbilliard_xy not defined for this billiard \n");
@ -1960,6 +2533,20 @@ int color[NPARTMAX];
return(vreuleaux(config, alpha, pos));
break;
}
case (D_FLOWER):
{
c = pos_flower(config, pos, &alpha);
return(vflower(config, alpha, pos));
break;
}
case (D_ALT_REU):
{
c = pos_alt_reuleaux(config, pos, &alpha);
return(valt_reuleaux(config, alpha, pos));
break;
}
default:
{
printf("Function vbilliard not defined for this billiard \n");
@ -1971,7 +2558,7 @@ int color[NPARTMAX];
/* returns 1 if (x,y) represents a point in the billiard */
double x, y;
{
double l2, r2, omega, omega2, c, angle, x1, y1, x2, co, so;
double l2, r1, r2, omega, omega2, c, angle, x1, y1, x2, co, so, x2plus, x2minus;
int condition, k;
switch (B_DOMAIN) {
@ -2021,8 +2608,9 @@ int color[NPARTMAX];
case D_ANNULUS:
{
l2 = LAMBDA*LAMBDA;
r2 = x*x + y*y;
if ((r2 > l2)&&(r2 < 1.0)) return(1);
r1 = x*x + y*y;
r2 = (x-MU)*(x-MU) + y*y;
if ((r2 > l2)&&(r1 < 1.0)) return(1);
else return(0);
break;
}
@ -2056,22 +2644,42 @@ int color[NPARTMAX];
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 */
case D_FLOWER:
{
/* TO DO */
return(1);
break;
}
case D_ALT_REU:
{
condition = 1;
omega2 = PI/((double)NPOLY);
co = cos(omega2);
so = sin(omega2);
x2plus = co + sqrt(LAMBDA*LAMBDA - so*so);
x2minus = 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 (k%2==0) condition = condition*((x1-x2plus)*(x1-x2plus) + y1*y1 > LAMBDA*LAMBDA);
else condition = condition*((x1-x2minus)*(x1-x2minus) + y1*y1 < LAMBDA*LAMBDA);
}
return(condition);
break;
}
default:
{
printf("Function ij_in_billiard not defined for this billiard \n");
return(0);
return(1);
}
}
}

331
sub_wave.c
View File

@ -198,6 +198,7 @@ void write_text( double x, double y, char *st)
/*********************/
/* some basic math */
/*********************/
@ -296,6 +297,22 @@ double xy[2];
xy[1] = YMIN + ((double)j)*(YMAX-YMIN)/((double)NY);
}
void draw_rectangle(x1, y1, x2, y2)
double x1, y1, x2, y2;
{
double pos[2];
glBegin(GL_LINE_LOOP);
xy_to_pos(x1, y1, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(x2, y1, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(x2, y2, pos);
glVertex2d(pos[0], pos[1]);
xy_to_pos(x1, y2, pos);
glVertex2d(pos[0], pos[1]);
glEnd();
}
@ -303,7 +320,7 @@ 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;
double l2, r2, r2mu, omega, c, angle, z, x1, y1, u, v, u1, v1;
int i, k, k1, k2, condition;
switch (B_DOMAIN) {
@ -399,6 +416,124 @@ double x, y;
{
return(((x-1.0)*(x-1.0) + y*y < LAMBDA*LAMBDA)||((x+1.0)*(x+1.0) + y*y < LAMBDA*LAMBDA)||((vabs(x) < 1.0)&&(vabs(y) < MU)));
}
case D_MENGER:
{
x1 = 0.5*(x+1.0);
y1 = 0.5*(y+1.0);
for (k=0; k<MDEPTH; k++)
{
x1 = x1*(double)MRATIO;
y1 = y1*(double)MRATIO;
if ((vabs(x)<1.0)&&(vabs(y)<1.0)&&(((int)x1 % MRATIO)==MRATIO/2)&&(((int)y1 % MRATIO)==MRATIO/2)) return(0);
}
return(1);
}
case D_JULIA_INT:
{
u = x/JULIA_SCALE;
v = y/JULIA_SCALE;
i = 0;
while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
{
u1 = u*u - v*v + julia_x;
v = 2.0*u*v + julia_y;
u = u1;
i++;
}
if (u*u + v*v < MANDELLIMIT) return(1);
else return(0);
}
case D_ANNULUS_HEATED: /* returns 2 if in inner circle */
{
l2 = LAMBDA*LAMBDA;
r2 = x*x + y*y;
r2mu = (x-MU)*(x-MU) + y*y;
if ((r2mu > l2)&&(r2 < 1.0)) return(1);
else if (r2mu <= l2) return(2);
else return (0);
}
case D_MENGER_HEATED:
{
if ((vabs(x) >= 1.0)||(vabs(y) >= 1.0)) return(0);
else
{
x1 = 0.5*(x+1.0);
y1 = 0.5*(y+1.0);
for (k=0; k<MDEPTH; k++)
{
x1 = x1*(double)MRATIO;
y1 = y1*(double)MRATIO;
if ((((int)x1 % MRATIO)==MRATIO/2)&&(((int)y1 % MRATIO)==MRATIO/2)) return(k+2);
}
return(1);
}
}
case D_MENGER_H_OPEN: /* returns 2 if in inner circle */
{
x1 = 0.5*(x+1.0);
y1 = 0.5*(y+1.0);
for (k=0; k<MDEPTH; k++)
{
x1 = x1*(double)MRATIO;
y1 = y1*(double)MRATIO;
if ((vabs(x)<1.0)&&(vabs(y)<1.0)&&(((int)x1 % MRATIO)==MRATIO/2)&&(((int)y1 % MRATIO)==MRATIO/2)) return(k+2);
}
return(1);
}
case D_MANDELBROT:
{
u = 0.0;
v = 0.0;
i = 0;
while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
{
u1 = u*u - v*v + x;
v = 2.0*u*v + y;
u = u1;
i++;
/* old version used */
/* u1 = u*u - v*v - x; */
/* v = 2.0*u*v - y; */
}
if (u*u + v*v < MANDELLIMIT) return(0);
else if ((x-0.5)*(x-0.5)/3.0 + y*y/1.0 > 1.2) return(2);
else return(1);
}
case D_MANDELBROT_CIRCLE:
{
u = 0.0;
v = 0.0;
i = 0;
while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
{
u1 = u*u - v*v + x;
v = 2.0*u*v + y;
u = u1;
i++;
}
if (u*u + v*v < MANDELLIMIT) return(0);
else if ((x-LAMBDA)*(x-LAMBDA) + (y-0.5)*(y-0.5) < MU*MU) return(2);
else return(1);
}
case D_JULIA:
{
u = x/JULIA_SCALE;
v = y/JULIA_SCALE;
i = 0;
while ((i<MANDELLEVEL)&&(u*u+v*v < 1000.0*MANDELLIMIT))
{
u1 = u*u - v*v + julia_x;
v = 2.0*u*v + julia_y;
u = u1;
i++;
// printf("x = %.5lg y = %.5lg i = %i r2 = %.5lg\n", x, y, i, u*u+v*v);
}
// printf("i = %i x = %.5lg y = %.5lg r2 = %.5lg\n", i, x, y, u*u+v*v);
if (u*u + v*v < MANDELLIMIT) return(0);
else if (x*x/3.0 + y*y/1.0 > 1.2) return(2);
// else if ((vabs(x) > XMAX - 0.01)||(vabs(y) > YMAX - 0.01)) return(2);
else return(1);
}
default:
{
printf("Function ij_in_billiard not defined for this billiard \n");
@ -423,8 +558,8 @@ int i, j;
void draw_billiard() /* draws the billiard boundary */
{
double x0, x, y, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z;
int i, j, k1, k2;
double x0, x, y, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l;
int i, j, k1, k2, mr2;
if (BLACK) glColor3f(1.0, 1.0, 1.0);
else glColor3f(0.0, 0.0, 0.0);
@ -710,7 +845,195 @@ void draw_billiard() /* draws the billiard boundary */
glEnd ();
break;
}
default:
case (D_MENGER):
{
glLineWidth(3);
// draw_rectangle(XMIN, -1.0, XMAX, 1.0);
/* level 1 */
if (MDEPTH > 0)
{
glLineWidth(2);
x = 1.0/((double)MRATIO);
draw_rectangle(x, x, -x, -x);
}
/* level 2 */
if (MDEPTH > 1)
{
glLineWidth(1);
mr2 = MRATIO*MRATIO;
l = 2.0/((double)mr2);
for (i=0; i<MRATIO; i++)
for (j=0; j<MRATIO; j++)
if ((i!=MRATIO/2)||(j!=MRATIO/2))
{
x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)MRATIO);
y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)MRATIO);
draw_rectangle(x, y, x+l, y+l);
}
}
/* level 3 */
if (MDEPTH > 2)
{
glLineWidth(1);
l = 2.0/((double)(mr2*MRATIO));
for (i=0; i<mr2; i++)
for (j=0; j<mr2; j++)
if ( (((i%MRATIO!=MRATIO/2))||(j%MRATIO!=MRATIO/2)) && (((i/MRATIO!=MRATIO/2))||(j/MRATIO!=MRATIO/2)) )
{
x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)mr2);
y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)mr2);
draw_rectangle(x, y, x+l, y+l);
}
}
break;
}
case (D_JULIA_INT):
{
/* Do nothing */
break;
}
case (D_ANNULUS_HEATED):
{
glBegin(GL_LINE_LOOP);
for (i=0; i<=NSEG; i++)
{
phi = (double)i*DPI/(double)NSEG;
x = MU + 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_MENGER_HEATED):
{
glLineWidth(3);
draw_rectangle(-1.0, -1.0, 1.0, 1.0);
/* level 1 */
if (MDEPTH > 0)
{
glLineWidth(2);
x = 1.0/((double)MRATIO);
draw_rectangle(x, x, -x, -x);
}
/* level 2 */
if (MDEPTH > 1)
{
glLineWidth(1);
mr2 = MRATIO*MRATIO;
l = 2.0/((double)mr2);
for (i=0; i<MRATIO; i++)
for (j=0; j<MRATIO; j++)
if ((i!=MRATIO/2)||(j!=MRATIO/2))
{
x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)MRATIO);
y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)MRATIO);
draw_rectangle(x, y, x+l, y+l);
}
}
/* level 3 */
if (MDEPTH > 2)
{
glLineWidth(1);
l = 2.0/((double)(mr2*MRATIO));
for (i=0; i<mr2; i++)
for (j=0; j<mr2; j++)
if ( (((i%MRATIO!=MRATIO/2))||(j%MRATIO!=MRATIO/2)) && (((i/MRATIO!=MRATIO/2))||(j/MRATIO!=MRATIO/2)) )
{
x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)mr2);
y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)mr2);
draw_rectangle(x, y, x+l, y+l);
}
}
break;
}
case (D_MENGER_H_OPEN):
{
glLineWidth(3);
// draw_rectangle(XMIN, -1.0, XMAX, 1.0);
/* level 1 */
if (MDEPTH > 0)
{
glLineWidth(2);
x = 1.0/((double)MRATIO);
draw_rectangle(x, x, -x, -x);
}
/* level 2 */
if (MDEPTH > 1)
{
glLineWidth(1);
mr2 = MRATIO*MRATIO;
l = 2.0/((double)mr2);
for (i=0; i<MRATIO; i++)
for (j=0; j<MRATIO; j++)
if ((i!=MRATIO/2)||(j!=MRATIO/2))
{
x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)MRATIO);
y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)MRATIO);
draw_rectangle(x, y, x+l, y+l);
}
}
/* level 3 */
if (MDEPTH > 2)
{
glLineWidth(1);
l = 2.0/((double)(mr2*MRATIO));
for (i=0; i<mr2; i++)
for (j=0; j<mr2; j++)
if ( (((i%MRATIO!=MRATIO/2))||(j%MRATIO!=MRATIO/2)) && (((i/MRATIO!=MRATIO/2))||(j/MRATIO!=MRATIO/2)) )
{
x = -1.0 - 0.5*l + 2.0*((double)i + 0.5)/((double)mr2);
y = -1.0 - 0.5*l + 2.0*((double)j + 0.5)/((double)mr2);
draw_rectangle(x, y, x+l, y+l);
}
}
break;
}
case (D_MANDELBROT):
{
/* Do nothing */
break;
}
case (D_MANDELBROT_CIRCLE):
{
/* Do nothing */
break;
}
case (D_JULIA):
{
/* Do nothing */
break;
}
default:
{
printf("Function draw_billiard not defined for this billiard \n");
}

204
wave_billiard.c
View File

@ -46,8 +46,10 @@
#define WINWIDTH 1280 /* window width */
#define WINHEIGHT 720 /* window height */
#define NX 640 /* number of grid points on x axis */
#define NY 360 /* number of grid points on y axis */
#define NX 1280 /* number of grid points on x axis */
#define NY 720 /* number of grid points on y axis */
// #define NX 640 /* number of grid points on x axis */
// #define NY 360 /* number of grid points on y axis */
/* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */
/* but will multiply run time by 4 */
@ -57,9 +59,11 @@
#define YMIN -1.125
#define YMAX 1.125 /* y interval for 9/16 aspect ratio */
#define JULIA_SCALE 1.1 /* scaling for Julia sets */
/* Choice of the billiard table */
#define B_DOMAIN 8 /* choice of domain shape */
#define B_DOMAIN 16 /* choice of domain shape */
#define D_RECTANGLE 0 /* rectangular domain */
#define D_ELLIPSE 1 /* elliptical domain */
@ -74,21 +78,39 @@
#define D_GRATING 10 /* diffraction grating */
#define D_EHRENFEST 11 /* Ehrenfest urn type geometry */
#define D_MENGER 15 /* Menger-Sierpinski carpet */
#define D_JULIA_INT 16 /* interior of Julia set */
/* Billiard tables for heat equation */
#define D_ANNULUS_HEATED 21 /* annulus with different temperatures */
#define D_MENGER_HEATED 22 /* Menger gasket with different temperatures */
#define D_MENGER_H_OPEN 23 /* Menger gasket with different temperatures and larger domain */
#define D_MANDELBROT 24 /* Mandelbrot set */
#define D_JULIA 25 /* Julia set */
#define D_MANDELBROT_CIRCLE 26 /* Mandelbrot set with circular conductor */
#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */
#define MU 0.05 /* parameter controlling the dimensions of domain */
#define NPOLY 8 /* number of sides of polygon */
#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */
#define MDEPTH 4 /* depth of computation of Menger gasket */
#define MRATIO 3 /* ratio defining Menger gasket */
#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */
#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */
#define FOCI 1 /* set to 1 to draw focal points of ellipse */
/* You can add more billiard tables by adapting the functions */
/* xy_in_billiard and draw_billiard below */
/* Physical patameters of wave equation */
/* Physical parameters of wave equation */
#define OMEGA 0.9 /* frequency of periodic excitation */
#define COURANT 0.01 /* Courant number */
#define GAMMA 0.0 /* damping factor in wave equation */
// #define GAMMA 5.0e-10 /* damping factor in wave equation */
#define KAPPA 5.0e-6 /* "elasticity" term enforcing oscillations */
#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */
// #define KAPPA 5.0e-6 /* "elasticity" term enforcing oscillations */
// #define KAPPA 5.0e-9 /* "elasticity" term enforcing oscillations */
// #define KAPPA 5.0e-8 /* "elasticity" term enforcing oscillations */
/* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
@ -96,13 +118,22 @@
/* Increasing COURANT speeds up the simulation, but decreases accuracy */
/* For similar wave forms, COURANT^2*GAMMA should be kept constant */
/* Boundary conditions */
#define B_COND 2
#define BC_DIRICHLET 0 /* Dirichlet boundary conditions */
#define BC_PERIODIC 1 /* periodic boundary conditions */
#define BC_ABSORBING 2 /* absorbing boundary conditions (beta version) */
/* For debugging purposes only */
#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */
#define VMAX 10.0 /* max value of wave amplitude */
/* Parameters for length and speed of simulation */
#define NSTEPS 5000 /* number of frames of movie */
#define NSTEPS 4000 /* number of frames of movie */
#define NVID 25 /* number of iterations between images displayed on screen */
#define NSEG 100 /* number of segments of boundary */
@ -120,7 +151,7 @@
#define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */
#define C_HUE 1 /* color scheme modifies hue */
#define SCALE 1 /* set to 1 to adjust color scheme to variance of field */
#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */
#define SLOPE 1.0 /* sensitivity of color on wave amplitude */
#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */
@ -128,8 +159,8 @@
#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 100.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP 80.0 /* amplitude of variation of hue for color scheme C_HUE */
#define HUEMEAN 230.0 /* mean value of hue for color scheme C_HUE */
#define HUEAMP 50.0 /* amplitude of variation of hue for color scheme C_HUE */
// #define HUEMEAN 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 */
@ -139,6 +170,10 @@
#define DPI 6.283185307
#define PID 1.570796327
double julia_x = -0.5, julia_y = 0.5; /* parameters for Julia sets */
// double julia_x = 0.33267, julia_y = 0.06395; /* parameters for Julia sets */
// double julia_x = 0.37468, julia_y = 0.21115; /* parameters for Julia sets */
#include "sub_wave.c"
double courant2; /* Courant parameter squared */
@ -162,6 +197,24 @@ void init_wave(x, y, phi, psi, xy_in)
}
}
void init_wave_flat(phi, psi, xy_in)
/* initialise flat field - phi is wave height, psi is phi at time t-1 */
double *phi[NX], *psi[NX]; short int * xy_in[NX];
{
int i, j;
double xy[2], dist2;
for (i=0; i<NX; i++)
for (j=0; j<NY; j++)
{
ij_to_xy(i, j, xy);
xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
phi[i][j] = 0.0;
psi[i][j] = 0.0;
}
}
void add_drop_to_wave(factor, x, y, phi, psi)
/* add drop at (x,y) to the field with given prefactor */
double factor, x, y, *phi[NX], *psi[NX];
@ -178,6 +231,33 @@ double factor, x, y, *phi[NX], *psi[NX];
}
}
void oscillate_linear_wave(amplitude, t, x1, y1, x2, y2, phi, psi)
/* oscillating boundary condition at (x,y) */
double amplitude, t, x1, y1, x2, y2, *phi[NX], *psi[NX];
{
int i, j, ij1[2], ij2[2], imin, imax, jmin, jmax, d = 5;
double xy[2], dist2;
xy_to_ij(x1, y1, ij1);
xy_to_ij(x2, y2, ij2);
imin = ij1[0] - d; if (imin < 0) imin = 0;
imax = ij2[0] + d; if (imax >= NX) imax = NX-1;
jmin = ij1[1] - d; if (jmin < 0) jmin = 0;
jmax = ij2[1] + d; if (jmax >= NY) jmax = NY-1;
for (i = imin; i < imax; i++)
for (j = jmin; j < jmax; j++)
{
ij_to_xy(i, j, xy);
dist2 = (xy[0]-x1)*(xy[0]-x1); /* to be improved */
// dist2 = (xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
// if (dist2 < 0.01)
if (dist2 < 0.001)
phi[i][j] = amplitude*exp(-dist2/0.001)*cos(-sqrt(dist2)/0.01)*cos(t*OMEGA);
// phi[i][j] += 0.2*exp(-dist2/0.001)*cos(-sqrt(dist2)/0.01)*cos(t*OMEGA);
}
}
@ -215,7 +295,92 @@ int time;
glEnd ();
}
void evolve_wave(phi, psi, xy_in)
void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in)
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
double *phi_in[NX], *psi_in[NX], *phi_out[NX], *psi_out[NX]; short int *xy_in[NX];
{
int i, j, iplus, iminus, jplus, jminus;
double delta, x, y, c, cc;
c = COURANT;
cc = courant2;
#pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y)
for (i=0; i<NX; i++){
for (j=0; j<NY; j++){
if (xy_in[i][j]){
/* discretized Laplacian */
if ((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;
}
delta = phi_in[iplus][j] + phi_in[iminus][j] + phi_in[i][jplus] + phi_in[i][jminus] - 4.0*phi_in[i][j];
x = phi_in[i][j];
y = psi_in[i][j];
/* evolve phi */
if (B_COND != BC_ABSORBING) phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
else
{
if ((i>0)&&(i<NX-1)&&(j>0)&&(j<NY-1))
phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y);
/* upper border */
else if (j==NY-1)
phi_out[i][j] = x - c*(x - phi_in[i][NY-2]) - KAPPA*x - GAMMA*(x-y);
/* lower border */
else if (j==0)
phi_out[i][j] = x - c*(x - phi_in[i][1]) - KAPPA*x - GAMMA*(x-y);
/* right border */
if (i==NX-1)
phi_out[i][j] = x - c*(x - phi_in[NX-2][j]) - KAPPA*x - GAMMA*(x-y);
/* left border */
else if (i==0)
phi_out[i][j] = x - c*(x - phi_in[1][j]) - KAPPA*x - GAMMA*(x-y);
}
psi_out[i][j] = x;
if (FLOOR)
{
if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX;
if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX;
if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX;
if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX;
}
}
}
}
// printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]);
}
void evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in)
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX]; short int *xy_in[NX];
{
evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in);
evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in);
}
void old_evolve_wave(phi, psi, xy_in)
/* time step of field evolution */
/* phi is value of field at time t, psi at time t-1 */
double *phi[NX], *psi[NX]; short int *xy_in[NX];
@ -286,7 +451,7 @@ double compute_variance(phi, psi, xy_in)
void animation()
{
double time, scale;
double *phi[NX], *psi[NX];
double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX];
short int *xy_in[NX];
int i, j, s;
@ -295,12 +460,15 @@ void animation()
{
phi[i] = (double *)malloc(NY*sizeof(double));
psi[i] = (double *)malloc(NY*sizeof(double));
phi_tmp[i] = (double *)malloc(NY*sizeof(double));
psi_tmp[i] = (double *)malloc(NY*sizeof(double));
xy_in[i] = (short int *)malloc(NY*sizeof(short int));
}
courant2 = COURANT*COURANT;
/* initialize wave with a drop at one point, zero elsewhere */
// init_wave_flat(phi, psi, xy_in);
init_wave(0.0, 0.0, phi, psi, xy_in);
/* add a drop at another point */
@ -329,12 +497,19 @@ void animation()
scale = sqrt(1.0 + compute_variance(phi,psi, xy_in));
// printf("Scaling factor: %5lg\n", scale);
}
else scale = 1.0;
// /* TO BE ADAPTED */
// scale = 1.0;
draw_wave(phi, psi, xy_in, scale, i);
for (j=0; j<NVID; j++) evolve_wave(phi, psi, xy_in);
for (j=0; j<NVID; j++)
{
evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
// if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
}
// for (j=0; j<NVID; j++) evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
draw_billiard();
@ -365,6 +540,9 @@ void animation()
{
free(phi[i]);
free(psi[i]);
free(phi_tmp[i]);
free(psi_tmp[i]);
free(xy_in[i]);
}
}