diff --git a/drop_billiard.c b/drop_billiard.c index 9a13700..116dca7 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -50,6 +50,8 @@ #define CIRCLE_PATTERN 0 /* pattern of circles */ +#define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ + #define NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */ // #define NCX 10 /* number of circles in x direction */ // #define NCY 15 /* number of circles in y direction */ @@ -126,10 +128,8 @@ /* animation part */ /*********************/ -void init_boundary_config(smin, smax, anglemin, anglemax, configs) +void init_boundary_config(double smin, double smax, double anglemin, double anglemax, double *configs[NPARTMAX]) /* initialize configuration: drop on the boundary, beta version */ -double smin, smax, anglemin, anglemax; -double *configs[NPARTMAX]; { int i; double ds, da, s, angle, alpha, pos[2], conf[2]; @@ -151,9 +151,8 @@ double *configs[NPARTMAX]; } } -void init_drop_config(x0, y0, angle1, angle2, configs) /* initialize configuration: drop at (x0,y0) */ -double x0, y0, angle1, angle2; -double *configs[NPARTMAX]; +void init_drop_config(double x0, double y0, double angle1, double angle2, double *configs[NPARTMAX]) +/* initialize configuration: drop at (x0,y0) */ { int i; double dalpha, alpha, pos[2]; @@ -172,9 +171,8 @@ double *configs[NPARTMAX]; } -int resample(color, configs) /* add particles where the front is stretched too thin */ -int color[NPARTMAX]; -double *configs[NPARTMAX]; +int resample(int color[NPARTMAX], double *configs[NPARTMAX]) +/* add particles where the front is stretched too thin */ { int len, i, j, k, iplus, newnparticles=nparticles, *newcolor; double dx, dy, pos[2], s1, s2, s, x, y, x1, y1, theta, alpha, beta, length2; @@ -261,10 +259,9 @@ double *configs[NPARTMAX]; else return(1); } -void draw_config(color, configs) -/* draw the wave front by ordering colors */ -int color[NPARTMAX]; -double *configs[NPARTMAX]; + +void draw_config(int color[NPARTMAX], double *configs[NPARTMAX]) +/* draw the particles */ { int i; double x1, y1, x2, y2, cosphi, sinphi, rgb[3], dist, dmax; @@ -315,10 +312,8 @@ double *configs[NPARTMAX]; } -void draw_ordered_config(color, configs) +void draw_ordered_config(int color[NPARTMAX], double *configs[NPARTMAX]) /* 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; @@ -376,10 +371,8 @@ double *configs[NPARTMAX]; } -void graph_movie(time, color, configs) +void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX]) /* compute next movie frame */ -int time, color[NPARTMAX]; -double *configs[NPARTMAX]; { int i, j; @@ -407,10 +400,8 @@ double *configs[NPARTMAX]; } } -void graph_no_movie(time, color, configs) +void graph_no_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX]) /* plot next image without making a movie */ -int time, color[NPARTMAX]; -double *configs[NPARTMAX]; { int i, j; diff --git a/global_particles.c b/global_particles.c index be7bb05..387fe5b 100644 --- a/global_particles.c +++ b/global_particles.c @@ -2,6 +2,8 @@ double circlex[NMAXCIRCLES], circley[NMAXCIRCLES], circlerad[NMAXCIRCLES]; short int circleactive[NMAXCIRCLES]; /* tells which circular scatters are active */ int ncircles = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */ +double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; +/* shooter and target positions for "laser in room of mirrors" simulations, with default values for square domain */ /* some basic math */ @@ -28,8 +30,12 @@ int ncircles = NMAXCIRCLES; /* actual number of circles, can be decre #define D_CIRCLES 20 /* several circles */ #define D_CIRCLES_IN_RECT 21 /* several circles inside a rectangle */ +#define D_CIRCLES_IN_GENUSN 22 /* several circles in polygon with identified opposite sides */ #define C_FOUR_CIRCLES 0 /* four circles almost touching each other */ #define C_SQUARE 1 /* square grid of circles */ #define C_HEX 2 /* hexagonal/triangular grid of circles */ -#define C_LASER 3 /* laser fight in a room of mirrors */ \ No newline at end of file +#define C_GOLDEN_MEAN 3 /* golden mean grid */ + +#define C_LASER 11 /* laser fight in a room of mirrors */ +#define C_LASER_GENUSN 12 /* laser fight in a translation surface */ \ No newline at end of file diff --git a/global_pdes.c b/global_pdes.c index 02eb333..916a2dc 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -42,6 +42,12 @@ #define C_CLOAK 5 /* invisibility cloak */ #define C_CLOAK_A 6 /* first optimized invisibility cloak */ +#define C_GOLDEN_MEAN 10 /* pattern based on vertical shifts by golden mean */ +#define C_GOLDEN_SPIRAL 11 /* spiral pattern based on golden mean */ +#define C_SQUARE_HEX 12 /* alternating between square and hexagonal/triangular */ + +#define C_ONE 97 /* one single circle, as for Sinai */ +#define C_TWO 98 /* two concentric circles of different type */ #define C_NOTHING 99 /* no circle at all, for comparisons */ @@ -60,10 +66,12 @@ #define BC_PERIODIC 1 /* periodic boundary conditions */ #define BC_ABSORBING 2 /* absorbing boundary conditions (beta version) */ #define BC_VPER_HABS 3 /* vertically periodic and horizontally absorbing boundary conditions */ +#define BC_ABS_REFLECT 4 /* absorbing boundary conditions, except reflecting at y=0, for comparisons */ +// #define BC_OSCILL_ABSORB 5 /* oscillating boundary condition on the left, absorbing on other walls */ /* For debugging purposes only */ -#define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ -#define VMAX 10.0 /* max value of wave amplitude */ +// #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ +// #define VMAX 10.0 /* max value of wave amplitude */ /* Plot types */ diff --git a/heat.c b/heat.c index c50a33f..4408d69 100644 --- a/heat.c +++ b/heat.c @@ -107,6 +107,7 @@ #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 BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PAUSE 100 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ @@ -160,11 +161,9 @@ double intstep1; /* integration step used in absorbing boundary conditions */ -void init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in) +void init_gaussian(double x, double y, double mean, double amplitude, double scalex, + double *phi[NX], short int * xy_in[NX]) /* 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; @@ -193,11 +192,8 @@ void init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in) } } -void init_julia_set(phi, xy_in) +void init_julia_set(double *phi[NX], short int * xy_in[NX]) /* change Julia set boundary condition */ - double *phi[NX]; - short int * xy_in[NX]; - { int i, j, in; double xy[2], dist2, module, phase, scale2; @@ -220,9 +216,8 @@ void init_julia_set(phi, xy_in) /*********************/ -void compute_gradient(phi, nablax, nablay) +void compute_gradient(double *phi[NX], double *nablax[NX], double *nablay[NX]) /* compute the gradient of the field */ -double *phi[NX], *nablax[NX], *nablay[NX]; { int i, j, iplus, iminus, jplus, jminus; double dx; @@ -241,12 +236,9 @@ double *phi[NX], *nablax[NX], *nablay[NX]; } } -void draw_field_line(x, y, xy_in, nablax, nablay, delta, nsteps) -// void draw_field_line(x, y, nablax, nablay, delta, nsteps) +void draw_field_line(double x, double y, short int *xy_in[NX], double *nablax[NX], + double *nablay[NX], double delta, int 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; @@ -303,11 +295,8 @@ short int *xy_in[NX]; glEnd(); } -void draw_wave(phi, xy_in, scale, time) +void draw_wave(double *phi[NX], short int *xy_in[NX], double scale, int 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; @@ -413,9 +402,8 @@ int time; -void evolve_wave_half(phi_in, phi_out, xy_in) +void evolve_wave_half(double *phi_in[NX], double *phi_out[NX], short int *xy_in[NX]) /* time step of field evolution */ - double *phi_in[NX], *phi_out[NX]; short int *xy_in[NX]; { int i, j, iplus, iminus, jplus, jminus; double delta1, delta2, x, y; @@ -493,18 +481,16 @@ void evolve_wave_half(phi_in, phi_out, xy_in) // printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]); } -void evolve_wave(phi, phi_tmp, xy_in) +void evolve_wave(double *phi[NX], double *phi_tmp[NX], short int *xy_in[NX]) /* time step of field evolution */ - double *phi[NX], *phi_tmp[NX]; short int *xy_in[NX]; { evolve_wave_half(phi, phi_tmp, xy_in); evolve_wave_half(phi_tmp, phi, xy_in); } -void old_evolve_wave(phi, xy_in) +void old_evolve_wave(double *phi[NX], short int *xy_in[NX]) /* 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]; @@ -595,9 +581,8 @@ void old_evolve_wave(phi, xy_in) // printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]); } -double compute_variance(phi, xy_in) +double compute_variance(double *phi[NX], short int * xy_in[NX]) /* 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; @@ -615,10 +600,8 @@ double *phi[NX]; short int * xy_in[NX]; return(variance/(double)n); } -void renormalise_field(phi, xy_in, variance) +void renormalise_field(double *phi[NX], short int * xy_in[NX], double variance) /* renormalise variance of field */ -double *phi[NX], variance; -short int * xy_in[NX]; { int i, j; double stdv; @@ -635,8 +618,7 @@ short int * xy_in[NX]; } } -void print_level(level) -int level; +void print_level(int level) { double pos[2]; char message[50]; @@ -660,10 +642,7 @@ void print_Julia_parameters() 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]; +void set_Julia_parameters(int time, double *phi[NX], short int *xy_in[NX]) { double jangle, cosj, sinj, radius = 0.15; @@ -680,10 +659,7 @@ short int *xy_in[NX]; 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]; +void set_Julia_parameters_cardioid(int time, double *phi[NX], short int *xy_in[NX]) { double jangle, cosj, sinj, yshift; diff --git a/particle_billiard.c b/particle_billiard.c index e0c9b17..e50e9d2 100644 --- a/particle_billiard.c +++ b/particle_billiard.c @@ -42,32 +42,19 @@ #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 -// #define YMAX 1.115 /* y interval for 9/16 aspect ratio */ - /* Choice of the billiard table, see global_particles.c */ -#define B_DOMAIN 14 /* choice of domain shape */ +#define B_DOMAIN 20 /* choice of domain shape */ -#define CIRCLE_PATTERN 0 /* pattern of circles */ +#define CIRCLE_PATTERN 2 /* pattern of circles */ + +#define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ #define NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */ -// #define NCX 10 /* number of circles in x direction */ -// #define NCY 15 /* number of circles in y direction */ #define NCX 15 /* number of circles in x direction */ #define NCY 20 /* number of circles in y direction */ #define LAMBDA 0.75 /* parameter controlling shape of billiard */ -// #define LAMBDA -3.346065215 /* sin(60°)/sin(15°) for Reuleaux-type triangle with 90° angles */ -// #define LAMBDA 3.0 /* parameter controlling shape of billiard */ -// #define LAMBDA 0.6 /* parameter controlling shape of billiard */ -// #define LAMBDA 0.4175295 /* sin(20°)/sin(55°) for 9-star shape with 30° angles */ -// #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.035 /* second parameter controlling shape of billiard */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NPOLY 8 /* number of sides of polygon */ @@ -81,16 +68,15 @@ /* Simulation parameters */ -#define NPART 5000 /* number of particles */ +#define NPART 30000 /* 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 SHOWTRAILS 0 /* set to 1 to keep trails of the particles */ #define NSTEPS 3000 /* number of frames of movie */ #define TIME 1000 /* time between movie frames, for fluidity of real-time simulation */ -// #define DPHI 0.000002 /* integration step */ -// #define DPHI 0.00002 /* integration step */ #define DPHI 0.000005 /* integration step */ #define NVID 150 /* number of iterations between images displayed on screen */ @@ -102,9 +88,9 @@ /* Colors and other graphical parameters */ -#define NCOLORS 16 /* number of colors */ +#define NCOLORS 32 /* number of colors */ #define COLORSHIFT 0 /* hue of initial color */ -#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */ +#define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ #define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */ #define NSEG 100 /* number of segments of boundary */ #define LENGTH 0.02 /* length of velocity vectors */ @@ -132,11 +118,9 @@ /* animation part */ /*********************/ -void init_boundary_config(smin, smax, anglemin, anglemax, configs) +void init_boundary_config(double smin, double smax, double anglemin, double anglemax, double *configs[NPARTMAX]) /* initialize configuration: drop on the boundary, beta version */ /* WORKS FOR ELLIPSE, HAS TO BE ADAPTED TO GENERAL BILLIARD */ -double smin, smax, anglemin, anglemax; -double *configs[NPARTMAX]; { int i; double ds, da, s, angle, theta, alpha, pos[2]; @@ -158,9 +142,8 @@ double *configs[NPARTMAX]; } } -void init_drop_config(x0, y0, angle1, angle2, configs) /* initialize configuration: drop at (x0,y0) */ -double x0, y0, angle1, angle2; -double *configs[NPARTMAX]; +void init_drop_config(double x0, double y0, double angle1, double angle2, double *configs[NPARTMAX]) +/* initialize configuration: drop at (x0,y0) */ { int i; double dalpha, alpha; @@ -181,10 +164,8 @@ double *configs[NPARTMAX]; } } -void init_sym_drop_config(x0, y0, angle1, angle2, configs) +void init_sym_drop_config(double x0, double y0, double angle1, double angle2, double *configs[NPARTMAX]) /* initialize configuration with two symmetric partial drops */ -double x0, y0, angle1, angle2; -double *configs[NPARTMAX]; { int i; double dalpha, alpha, meanangle; @@ -210,9 +191,8 @@ double *configs[NPARTMAX]; } -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 *configs[NPARTMAX]; +void init_line_config(double x0, double y0, double x1, double y1, double angle, double *configs[NPARTMAX]) +/* initialize configuration: line (x0,y0)-(x1,y1) in direction alpha */ { int i; double dx, dy; @@ -231,16 +211,14 @@ double *configs[NPARTMAX]; } -void draw_config(color, configs, active) +void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTMAX]) /* draw the particles */ -int color[NPARTMAX], active[NPARTMAX]; -double *configs[NPARTMAX]; { int i; double x1, y1, x2, y2, cosphi, sinphi, rgb[3]; glutSwapBuffers(); - blank(); + if (!SHOWTRAILS) blank(); if (PAINT_INT) paint_billiard_interior(); glLineWidth(PARTICLE_WIDTH); @@ -335,10 +313,8 @@ double *configs[NPARTMAX]; } -void graph_movie(time, color, configs, active) +void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTMAX]) /* compute next movie frame */ -int time, color[NPARTMAX], active[NPARTMAX]; -double *configs[NPARTMAX]; { int i, j, c; @@ -370,77 +346,6 @@ double *configs[NPARTMAX]; // draw_config(color, configs); } - -void init_circle_config() -{ - int i, j, n; - double dx, dy; - - switch (CIRCLE_PATTERN) { - case (C_FOUR_CIRCLES): - { - ncircles = 4; - - circlex[0] = 1.0; - circley[0] = 0.0; - circlerad[0] = 0.8; - - circlex[1] = -1.0; - circley[1] = 0.0; - circlerad[1] = 0.8; - - circlex[2] = 0.0; - circley[2] = 0.8; - circlerad[2] = 0.4; - - circlex[3] = 0.0; - circley[3] = -0.8; - circlerad[3] = 0.4; - - for (i=0; i<4; i++) circleactive[i] = 1; - - break; - } - case (C_SQUARE): - { - ncircles = NCX*NCY; - dy = (YMAX - YMIN)/((double)NCY); - for (i = 0; i < NCX; i++) - for (j = 0; j < NCY; j++) - { - n = NCY*i + j; - circlex[n] = ((double)(i-NCX/2) + 0.5)*dy; - circley[n] = YMIN + ((double)j + 0.5)*dy; - circlerad[n] = MU; - circleactive[n] = 1; - } - break; - } - case (C_HEX): - { - ncircles = NCX*(NCY+1); - dy = (YMAX - YMIN)/((double)NCY); - dx = dy*0.5*sqrt(3.0); - for (i = 0; i < NCX; i++) - for (j = 0; j < NCY+1; j++) - { - n = (NCY+1)*i + j; - circlex[n] = ((double)(i-NCX/2) + 0.5)*dy; - circley[n] = YMIN + ((double)j - 0.5)*dy; - if ((i+NCX)%2 == 1) circley[n] += 0.5*dy; - circlerad[n] = MU; - circleactive[n] = 1; - } - break; - } - default: - { - printf("Function init_circle_config not defined for this pattern \n"); - } - } -} - - void animation() { double time, dt, alpha, r; @@ -456,24 +361,27 @@ void animation() configs[i] = (double *)malloc(8*sizeof(double)); /* init circle configuration if the domain is D_CIRCLES */ - if (B_DOMAIN == D_CIRCLES) init_circle_config(); + if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)||(B_DOMAIN == D_CIRCLES_IN_GENUSN)) init_circle_config(); /* initialize system by putting particles in a given point with a range of velocities */ r = cos(PI/(double)NPOLY)/cos(DPI/(double)NPOLY); -// init_drop_config(0.0, 0.0, 0.0, PI, configs); +// init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0, configs); +// init_drop_config(-0.75, 0.0, -0.1, 0.1, configs); // init_drop_config(0.5, 0.5, -1.0, 1.0, configs); // init_sym_drop_config(-1.0, 0.5, -PID, PID, configs); // init_drop_config(-0.999, 0.0, -alpha, alpha, configs); // other possible initial conditions : -// init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0, configs); - init_line_config(0.0, -1.0, -1.0, 1.0, 0.25*PID, configs); + init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0, configs); +// init_line_config(0.0, -0.5, 0.0, 0.5, 0.0, configs); +// init_line_config(-1.25, -0.5, -1.25, 0.5, 0.0*PID, configs); +// init_line_config(-1.0, -0.3, -1.0, 0.3, 0.0, configs); // init_line_config(-0.7, -0.45, -0.7, 0.45, 0.0, configs); // init_line_config(-1.5, 0.1, -0.1, 1.0, -0.5*PID, configs); - blank(); + if (!SHOWTRAILS) blank(); glColor3f(0.0, 0.0, 0.0); if (DRAW_BILLIARD) draw_billiard(); @@ -556,9 +464,13 @@ void display(void) glPushMatrix(); blank(); - glutSwapBuffers(); - blank(); - glutSwapBuffers(); + + if (!SHOWTRAILS) + { + glutSwapBuffers(); + blank(); + glutSwapBuffers(); + } animation(); @@ -571,7 +483,9 @@ void display(void) int main(int argc, char** argv) { glutInit(&argc, argv); - glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH); + if (SHOWTRAILS) glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE | GLUT_DEPTH); + else glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH); +// glutInitDisplayMode(GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH); glutInitWindowSize(WINWIDTH,WINHEIGHT); glutCreateWindow("Billiard animation"); diff --git a/schrodinger.c b/schrodinger.c index 77430c0..fbf49b1 100644 --- a/schrodinger.c +++ b/schrodinger.c @@ -97,6 +97,7 @@ #define NSTEPS 200 /* number of frames of movie */ #define NVID 1200 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ +#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ @@ -140,12 +141,10 @@ double intstep1; /* integration step used in absorbing boundary conditions */ -void init_coherent_state(x, y, px, py, scalex, phi, psi, xy_in) +void init_coherent_state(double x, double y, double px, double py, double scalex, double *phi[NX], + double *psi[NX], short int *xy_in[NX]) /* 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; @@ -181,9 +180,9 @@ void init_coherent_state(x, y, px, py, scalex, phi, psi, xy_in) /* animation part */ /*********************/ -void schrodinger_color_scheme(phi, psi, scale, time, rgb) -double phi, psi, scale, rgb[3]; -int time; +void schrodinger_color_scheme(double phi, double psi, double scale, int time, double rgb[3]) +// double phi, psi, scale, rgb[3]; +// int time; { double phase, amp, lum; @@ -204,11 +203,8 @@ int time; } -void draw_wave(phi, psi, xy_in, scale, time) +void draw_wave(double *phi[NX], double *psi[NX], short int *xy_in[NX], double scale, int 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; @@ -234,10 +230,12 @@ int time; glEnd (); } -void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in) -/* time step of field evolution */ -/* phi is real part, psi is imaginary part */ - double *phi_in[NX], *psi_in[NX], *phi_out[NX], *psi_out[NX]; short int *xy_in[NX]; +void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], + short int *xy_in[NX]) +// void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in) +// /* time step of field evolution */ +// /* phi is real part, psi is imaginary part */ +// double *phi_in[NX], *psi_in[NX], *phi_out[NX], *psi_out[NX]; short int *xy_in[NX]; { int i, j, iplus, iminus, jplus, jminus; double delta1, delta2, x, y; @@ -325,19 +323,19 @@ void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in) // 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) +void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *psi_tmp[NX], short int *xy_in[NX]) /* 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]; +/* phi is real part, psi is imaginary part */ { evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in); evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in); } -double compute_variance(phi, psi, xy_in) +double compute_variance(double *phi[NX], double *psi[NX], short int *xy_in[NX]) +// 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]; +// double *phi[NX], *psi[NX]; short int * xy_in[NX]; { int i, j, n = 0; double variance = 0.0; @@ -355,10 +353,8 @@ double *phi[NX], *psi[NX]; short int * xy_in[NX]; return(variance/(double)n); } -void renormalise_field(phi, psi, xy_in, variance) +void renormalise_field(double *phi[NX], double *psi[NX], short int *xy_in[NX], double variance) /* renormalise variance of field */ -double *phi[NX], *psi[NX], variance; -short int * xy_in[NX]; { int i, j; double stdv; diff --git a/sub_part_billiard.c b/sub_part_billiard.c index 07658d8..9f29bfa 100644 --- a/sub_part_billiard.c +++ b/sub_part_billiard.c @@ -1,13 +1,14 @@ +#define DUMMY_ABSORBING -1000.0 /* dummy value of config[0] for absorbing circles */ +#define BOUNDARY_SHIFT 100.0 /* shift of boundary parametrisation for circles in domain */ + long int global_time = 0; /* counter to keep track of global time of simulation */ int nparticles=NPART; - /*********************/ /* some basic math */ /*********************/ - double vabs(x) /* absolute value */ - double x; + double vabs(double x) /* absolute value */ { double res; @@ -16,9 +17,7 @@ int nparticles=NPART; return(res); } - double module2(x, y) /* Euclidean norm */ - double x, y; - + double module2(double x, double y) /* Euclidean norm */ { double m; @@ -26,9 +25,7 @@ int nparticles=NPART; return(m); } - double argument(x, y) - double x, y; - + double argument(double x, double y) { double alph; @@ -48,8 +45,7 @@ int nparticles=NPART; return(alph); } - int polynome(a, b, c, r) - double a, b, c, r[2]; + int polynome(double a, double b, double c, double r[2]) { double delta, rdelta; int im = 1; @@ -137,9 +133,8 @@ void init() /* initialisation of window */ } -void hsl_to_rgb(h, s, l, rgb) /* color conversion from HSL to RGB */ +void hsl_to_rgb(double h, double s, double l, double rgb[3]) /* 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; @@ -177,9 +172,7 @@ double h, s, l, rgb[3]; } } -void rgb_color_scheme(i, rgb) /* color scheme */ -int i; -double rgb[3]; +void rgb_color_scheme(int i, double rgb[3]) /* color scheme */ { double hue, y, r; @@ -226,6 +219,20 @@ void save_frame() } +void write_text_fixedwidth( 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 + glutBitmapCharacter(GLUT_BITMAP_9_BY_15, st[i]); // Print a character on the screen + } +} + void write_text( double x, double y, char *st) { int l,i; @@ -239,10 +246,103 @@ void write_text( double x, double y, char *st) } } +void erase_area(double x, double y, double dx, double dy, double rgb[3]) +{ + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_QUADS); + glVertex2d(x - dx, y - dy); + glVertex2d(x + dx, y - dy); + glVertex2d(x + dx, y + dy); + glVertex2d(x - dx, y + dy); + glEnd(); +} - void compute_flower_parameters(omega, co, so, axis1, axis2, phimax) +void erase_rectangle_outside(double h, double s, double l) +{ + double rgb[3], dx; + int k; + + dx = 0.5*(XMAX - LAMBDA); + hsl_to_rgb(h, s, l, rgb); + erase_area(0.0, 1.1, 2.0, 0.1, rgb); + erase_area(0.0, -1.1, 2.0, 0.1, rgb); + erase_area(LAMBDA + dx, 0.0, dx, 2.0, rgb); + erase_area(-LAMBDA - dx, 0.0, dx, 2.0, rgb); +} + +void draw_circle(double x, double y, double r, int nseg) +{ + int i; + double phi, dphi, x1, y1; + + dphi = DPI/(double)nseg; + + glEnable(GL_LINE_SMOOTH); + glBegin(GL_LINE_LOOP); + for (i=0; i<=nseg; i++) + { + phi = (double)i*dphi; + x1 = x + r*cos(phi); + y1 = y + r*sin(phi); + glVertex2d(x1, y1); + } + glEnd (); +} + +void draw_colored_circle(double x, double y, double r, int nseg, double rgb[3]) +{ + int i, ij[2]; + double pos[2], alpha, dalpha; + + dalpha = DPI/(double)nseg; + +// glLineWidth(2); + + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_TRIANGLE_FAN); + glVertex2d(x,y); + for (i=0; i<=nseg; i++) + { + alpha = (double)i*dalpha; + glVertex2d(x + r*cos(alpha), y + r*sin(alpha)); + } + + glEnd(); +} + + + +int in_circle(double x, double y, double r) +/* test whether (x,y) is in circle of radius r */ +{ + return(x*x + y*y < r*r); +} + +int out_circle(double x, double y, double r) +/* test whether (x,y) is in circle of radius r */ +{ + return(x*x + y*y > r*r); +} + +int in_polygon(double x, double y, double r, int npoly, double apoly) +/* test whether (x,y) is in regular polygon of npoly sides inscribed in circle of radious r, turned by apoly Pi/2 */ +{ + int condition = 1, k; + double omega, cw, angle; + + omega = DPI/((double)npoly); + cw = cos(omega*0.5); + for (k=0; k= 1) { if (circleactive[k] == 2) { - hsl_to_rgb(220.0, 0.9, 0.5, rgb); + hsl_to_rgb(150.0, 0.9, 0.4, rgb); glColor3f(rgb[0], rgb[1], rgb[2]); } - glBegin(GL_LINE_LOOP); - for (i=0; i<=NSEG; i++) - { - phi = (double)i*DPI/(double)NSEG; - x = circlex[k] + circlerad[k]*cos(phi); - y = circley[k] + circlerad[k]*sin(phi); - glVertex2d(x, y); - } - glEnd (); - + draw_circle(circlex[k], circley[k], circlerad[k], NSEG); init_billiard_color(); } -// /* draw shooter position for laser pattern */ -// if (CIRCLE_PATTERN == C_LASER) -// { -// hsl_to_rgb(0.0, 0.9, 0.5, rgb); -// glColor3f(rgb[0], rgb[1], rgb[2]); -// -// glBegin(GL_LINE_LOOP); -// for (i=0; i<=NSEG; i++) -// { -// phi = (double)i*DPI/(double)NSEG; -// x = X_SHOOTER + circlerad[ncircles-1]*cos(phi); -// y = Y_SHOOTER + circlerad[ncircles-1]*sin(phi); -// glVertex2d(x, y); -// } -// glEnd (); -// } + /* draw shooter position for laser pattern */ + if (CIRCLE_PATTERN == C_LASER) + { + hsl_to_rgb(0.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + draw_circle(x_shooter, y_shooter, circlerad[ncircles-1], NSEG); + } init_billiard_color(); @@ -903,6 +938,44 @@ void draw_billiard() /* draws the billiard boundary */ glEnd(); break; } + case (D_CIRCLES_IN_GENUSN): + { +// for (k=0; k= 1) + for (k=0; k= 1)&&(in_polygon(circlex[k], circley[k], 1.0, NPOLY, APOLY))) + { + if (circleactive[k] == 2) + { + hsl_to_rgb(150.0, 0.9, 0.4, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + } + draw_circle(circlex[k], circley[k], circlerad[k], NSEG); + init_billiard_color(); + } + + /* draw shooter position for laser pattern */ + if ((CIRCLE_PATTERN == C_LASER)||(CIRCLE_PATTERN == C_LASER_GENUSN)) + { + hsl_to_rgb(0.0, 0.9, 0.5, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + draw_circle(x_shooter, y_shooter, circlerad[ncircles-1], NSEG); + } + + /* draw polygon */ + init_billiard_color(); + + 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(SCALING_FACTOR*x, SCALING_FACTOR*y); + } + glEnd (); + break; + } default: { printf("Function draw_billiard not defined for this billiard \n"); @@ -928,20 +1001,17 @@ void draw_billiard() /* draws the billiard boundary */ */ -void print_config(conf) /* for debugging purposes */ -double conf[8]; +void print_config(double conf[8]) /* for debugging purposes */ { 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]; +void print_config_23(double conf[8]) /* for debugging purposes */ { printf("t = %.8f, L = %.8f\n", conf[2], conf[3]); } -void print_colors(color) /* for debugging purposes */ -int color[NPARTMAX]; +void print_colors(int color[NPARTMAX]) /* for debugging purposes */ { int i; @@ -950,17 +1020,16 @@ int color[NPARTMAX]; } + /****************************************************************************************/ /* rectangle billiard */ /****************************************************************************************/ - int pos_rectangle(conf, pos, alpha) + int pos_rectangle(double conf[2], double pos[2], double *alpha) /* determine position on boundary of rectangle */ /* the corners of the rectangle are (-LAMBDA,-1), ..., (LAMBDA,1) */ /* returns the number of the side hit, or 0 if hitting a corner */ - double conf[2], pos[2], *alpha; - { double s, theta; @@ -1045,9 +1114,8 @@ int color[NPARTMAX]; } - int vrectangle_xy(config, alpha, pos) + int vrectangle_xy(double config[8], double alpha, double pos[2]) /* determine initial configuration for start at point pos = (x,y) */ - double config[8], alpha, pos[2]; { double l, s0, c0, x1, y1, margin = 1e-12; int c, intb=1; @@ -1126,9 +1194,7 @@ int color[NPARTMAX]; return(c); } - int vrectangle(config) - double config[8]; - + int vrectangle(double config[8]) { double pos[2], alpha; int c; @@ -1147,10 +1213,8 @@ int color[NPARTMAX]; /* elliptic billiard */ /****************************************************************************************/ - int pos_ellipse(conf, pos, alpha) + int pos_ellipse(double conf[2], double pos[2], double *alpha) /* determine position on boundary of ellipse */ - double conf[2], pos[2], *alpha; - { double theta; @@ -1164,10 +1228,8 @@ int color[NPARTMAX]; } - int vellipse_xy(config, alpha, pos) + int vellipse_xy(double config[8], double alpha, double pos[2]) /* determine initial configuration for start at point pos = (x,y) */ - double config[8], alpha, pos[2]; - { double c0, s0, lam2, a, b, c, x1, y1, t, theta; int i; @@ -1206,10 +1268,8 @@ int color[NPARTMAX]; return(1); } - int vellipse(config) + int vellipse(double config[8]) /* determine initial configuration when starting from boundary */ - double config[8]; - { double pos[2], theta, alpha; int i; @@ -1229,10 +1289,8 @@ int color[NPARTMAX]; /* stadium billiard */ /****************************************************************************************/ - int pos_stade(conf, pos, alpha) + int pos_stade(double conf[2], double pos[2], double *alpha) /* determine position on boundary of stadium */ - double conf[2], pos[2], *alpha; - { double s, theta, l, psi0, psi; @@ -1294,10 +1352,8 @@ int color[NPARTMAX]; } - int vstade_xy(config, alpha, pos) + int vstade_xy(double config[8], double alpha, double pos[2]) /* determine initial configuration for start at point pos = (x,y) */ - double config[8], alpha, pos[2]; - { double l, s0, c0, t, x, y, x1, y1, a, b, res[2]; double smin, psi, margin = 1e-12; @@ -1387,8 +1443,7 @@ int color[NPARTMAX]; return(c); } - int vstade(config) - double config[8]; + int vstade(double config[8]) { double alpha, pos[2]; int c; @@ -1404,11 +1459,9 @@ int color[NPARTMAX]; /* Sinai billiard */ /****************************************************************************************/ - int pos_sinai(conf, pos, alpha) + int pos_sinai(double conf[2], double pos[2], double *alpha) /* determine position on boundary of Sinai billiard */ /* s in [0,2 Pi) is on the circle, other s are on boundary of window */ - double conf[2], pos[2], *alpha; - { double s, theta, psi0, psi, s1, s2, s3, s4; @@ -1460,10 +1513,8 @@ int color[NPARTMAX]; } - int vsinai_xy(config, alpha, pos) + int vsinai_xy(double config[8], double alpha, double pos[2]) /* 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], s1, s2, s3, s4, s, r; double psi, lam2, margin = 1e-12; @@ -1582,9 +1633,7 @@ int color[NPARTMAX]; /* return a negative value if the disc is not hit, for color scheme */ } - int vsinai(config) - double config[8]; - + int vsinai(double config[8]) { double alpha, pos[2]; int c; @@ -1604,12 +1653,10 @@ int color[NPARTMAX]; /****************************************************************************************/ - int pos_triangle(conf, pos, alpha) + int pos_triangle(double conf[2], double pos[2], double *alpha) /* determine position on boundary of triangle */ /* the corners of the triangle are (-LAMBDA,-1), (LAMBDA,-1), (-LAMBDA,1) */ /* we use arclength for horizontal and vertical side, x for diagonal */ - double conf[2], pos[2], *alpha; - { double s, theta; @@ -1641,11 +1688,9 @@ int color[NPARTMAX]; } - int vtriangle_xy(config, alpha, pos) + int vtriangle_xy(double config[8], double alpha, double pos[2]) /* determine initial configuration for start at point pos = (x,y) */ /* Warning: reflection in the corners is not yet implemented correctly */ - double config[8], alpha, pos[2]; - { double s0, c0, t, x, y, x1, y1, psi; int c, intb=1, intc, i; @@ -1706,9 +1751,7 @@ int color[NPARTMAX]; return(c); } - int vtriangle(config) - double config[8]; - + int vtriangle(double config[8]) { double alpha, pos[2]; int c; @@ -1727,10 +1770,8 @@ int color[NPARTMAX]; /* annulus billiard */ /****************************************************************************************/ - int pos_annulus(conf, pos, alpha) + int pos_annulus(double conf[2], double pos[2], double *alpha) /* determine position on boundary of annulus */ - double conf[2], pos[2], *alpha; - { double s, theta, psi0, psi, s1, s2, s3, s4; @@ -1758,10 +1799,8 @@ int color[NPARTMAX]; } } - int vannulus_xy(config, alpha, pos) + int vannulus_xy(double config[8], double alpha, double pos[2]) /* 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; @@ -1825,10 +1864,8 @@ int color[NPARTMAX]; return(c); } - int vannulus(config) + int vannulus(double config[8]) /* determine initial configuration when starting from boundary */ - double config[8]; - { double pos[2], alpha; int c; @@ -1844,11 +1881,9 @@ int color[NPARTMAX]; /* polygonal billiard */ /****************************************************************************************/ - int pos_polygon(conf, pos, alpha) + int pos_polygon(double conf[2], double pos[2], double *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; @@ -1876,10 +1911,8 @@ int color[NPARTMAX]; return(c); } - int vpolygon_xy(config, alpha, pos) + int vpolygon_xy(double config[8], double alpha, double pos[2]) /* 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; @@ -1932,10 +1965,8 @@ int color[NPARTMAX]; return(c); } - int vpolygon(config) + int vpolygon(double config[8]) /* determine initial configuration when starting from boundary */ - double config[8]; - { double pos[2], alpha; int c; @@ -1951,11 +1982,9 @@ int color[NPARTMAX]; /* Reuleaux-type and star-shaped billiard */ /****************************************************************************************/ - int pos_reuleaux(conf, pos, alpha) + int pos_reuleaux(double conf[2], double pos[2], double *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; @@ -1989,10 +2018,8 @@ int color[NPARTMAX]; return(c); } - int vreuleaux_xy(config, alpha, pos) + int vreuleaux_xy(double config[8], double alpha, double pos[2]) /* 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]; @@ -2075,10 +2102,8 @@ int color[NPARTMAX]; return(c); } - int vreuleaux(config) + int vreuleaux(double config[8]) /* determine initial configuration when starting from boundary */ - double config[8]; - { double pos[2], alpha; int c; @@ -2094,11 +2119,9 @@ int color[NPARTMAX]; /* Bunimovich flower billiard */ /****************************************************************************************/ - int pos_flower(conf, pos, alpha) + int pos_flower(double conf[2], double pos[2], double *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; @@ -2128,10 +2151,9 @@ int color[NPARTMAX]; return(c); } - int vflower_xy(config, alpha, pos) + int vflower_xy(double config[8], double alpha, double pos[2]) /* determine initial configuration for start at point pos = (x,y) */ - double config[8], alpha, pos[2]; - +// 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; @@ -2201,105 +2223,9 @@ int color[NPARTMAX]; 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 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= 0) + { + ncirc = (int)(conf[0]/DPI); + if (ncirc >= ncircles) ncirc = ncircles - 1; + + angle = conf[0] - (double)ncirc*DPI; + + pos[0] = circlex[ncirc] + circlerad[ncirc]*cos(angle); + pos[1] = circley[ncirc] + circlerad[ncirc]*sin(angle); + + *alpha = angle + PID + conf[1]; + + return(ncirc); + } + else /* particle starts on boundary */ + { + omega = DPI/((double)NPOLY); + length = 2.0*sin(0.5*omega); + +// conf[0] += 2.0*length*(double)NPOLY; + conf[0] += BOUNDARY_SHIFT; + c = pos_genusn(conf, pos, alpha); + +// conf[0] -= 2.0*length*(double)NPOLY; + conf[0] -= BOUNDARY_SHIFT; + + return(-c); + } + } + + + int vcircles_in_genusn_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double c0, s0, b, c, t, theta, delta, margin = 1.0e-12, tmin, rlarge = 1.0e10, omega, length, angle, cw; + double tval[ncircles], xint[ncircles], yint[ncircles], phiint[ncircles]; + int i, k, nt = 0, nscat[ncircles], ntmin, side, condition; + + c0 = cos(alpha); + s0 = sin(alpha); + omega = DPI/((double)NPOLY); + length = 2.0*sin(0.5*omega); + cw = cos(omega*0.5); + + for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nscat[nt] = i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + phiint[nt] = argument(xint[nt] - circlex[i], yint[nt] - circley[i]); + + /* test wether intersection is in polygon */ + if (in_polygon(xint[nt], yint[nt], 1.0, NPOLY, APOLY)) nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= DPI) phiint[ntmin] -= DPI; + + config[0] = (double)nscat[ntmin]*DPI + phiint[ntmin]; + config[1] = PID - alpha + phiint[ntmin]; /* CHECK */ + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + + config[2] = 0.0; /* running time */ + config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = xint[ntmin]; /* position of collision */ + config[7] = yint[ntmin]; + + + /* set dummy coordinates if circles are absorbing */ + if (ABSORBING_CIRCLES) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + + return(nscat[ntmin]); + } + else /* there is no intersection with the circles - compute intersection with boundary */ + { + side = vgenusn_xy(config, alpha, pos); + +// config[0] -= 2.0*length*(double)NPOLY; + config[0] -= BOUNDARY_SHIFT; + + return(side); + } + } + + int vcircles_in_genusn(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_circles_in_genusn(config, pos, &alpha); + + vcircles_in_genusn_xy(config, alpha, pos); + + return(c); + } + + /****************************************************************************************/ /* general billiard */ /****************************************************************************************/ - int pos_billiard(conf, pos, alpha) + int pos_billiard(double conf[8], double pos[2], double *alpha) /* determine initial configuration for start at point pos = (x,y) */ - double conf[8], pos[2], *alpha; { switch (B_DOMAIN) { case (D_RECTANGLE): { - return(pos_rectangle(conf, pos, &alpha)); + return(pos_rectangle(conf, pos, alpha)); break; } case (D_ELLIPSE): { - return(pos_ellipse(conf, pos, &alpha)); + return(pos_ellipse(conf, pos, alpha)); break; } case (D_STADIUM): { - return(pos_stade(conf, pos, &alpha)); + return(pos_stade(conf, pos, alpha)); break; } case (D_SINAI): { - return(pos_sinai(conf, pos, &alpha)); + return(pos_sinai(conf, pos, alpha)); break; } case (D_TRIANGLE): { - return(pos_triangle(conf, pos, &alpha)); + return(pos_triangle(conf, pos, alpha)); break; } case (D_ANNULUS): { - return(pos_annulus(conf, pos, &alpha)); + return(pos_annulus(conf, pos, alpha)); break; } case (D_POLYGON): { - return(pos_polygon(conf, pos, &alpha)); + return(pos_polygon(conf, pos, alpha)); break; } case (D_REULEAUX): { - return(pos_reuleaux(conf, pos, &alpha)); + return(pos_reuleaux(conf, pos, alpha)); break; } case (D_FLOWER): { - return(pos_flower(conf, pos, &alpha)); + return(pos_flower(conf, pos, alpha)); break; } case (D_ALT_REU): { - return(pos_alt_reuleaux(conf, pos, &alpha)); + return(pos_alt_reuleaux(conf, pos, alpha)); break; } case (D_ANGLE): { - return(pos_angle(conf, pos, &alpha)); + return(pos_angle(conf, pos, alpha)); break; } case (D_LSHAPE): { - return(pos_lshape(conf, pos, &alpha)); + return(pos_lshape(conf, pos, alpha)); break; } case (D_GENUSN): { - return(pos_genusn(conf, pos, &alpha)); + return(pos_genusn(conf, pos, alpha)); break; } case (D_CIRCLES): { - return(pos_circles(conf, pos, &alpha)); + return(pos_circles(conf, pos, alpha)); break; } case (D_CIRCLES_IN_RECT): { - return(pos_circles_in_rect(conf, pos, &alpha)); + return(pos_circles_in_rect(conf, pos, alpha)); + break; + } + case (D_CIRCLES_IN_GENUSN): + { + return(pos_circles_in_genusn(conf, pos, alpha)); break; } default: @@ -3299,9 +3348,8 @@ int color[NPARTMAX]; - int vbilliard_xy(config, alpha, pos) + int vbilliard_xy(double config[8], double alpha, double pos[2]) /* determine initial configuration for start at point pos = (x,y) */ - double config[8], alpha, pos[2]; { switch (B_DOMAIN) { case (D_RECTANGLE): @@ -3379,6 +3427,11 @@ int color[NPARTMAX]; return(vcircles_in_rect_xy(config, alpha, pos)); break; } + case (D_CIRCLES_IN_GENUSN): + { + return(vcircles_in_genusn_xy(config, alpha, pos)); + break; + } default: { printf("Function vbilliard_xy not defined for this billiard \n"); @@ -3388,9 +3441,8 @@ int color[NPARTMAX]; /* TO DO: fix returned value */ - int vbilliard(config) + int vbilliard(double config[8]) /* determine initial configuration when starting from boundary */ - double config[8]; { double pos[2], theta, alpha; int c; @@ -3400,105 +3452,112 @@ int color[NPARTMAX]; { c = pos_rectangle(config, pos, &alpha); - return(vrectangle(config, alpha, pos)); + return(vrectangle(config)); break; } case (D_ELLIPSE): { c = pos_ellipse(config, pos, &alpha); - return(vellipse(config, alpha, pos)); + return(vellipse(config)); break; } case (D_STADIUM): { c = pos_stade(config, pos, &alpha); - return(vstade(config, alpha, pos)); + return(vstade(config)); break; } case (D_SINAI): { c = pos_sinai(config, pos, &alpha); - return(vsinai(config, alpha, pos)); + return(vsinai(config)); break; } case (D_TRIANGLE): { c = pos_triangle(config, pos, &alpha); - return(vtriangle(config, alpha, pos)); + return(vtriangle(config)); break; } case (D_ANNULUS): { c = pos_annulus(config, pos, &alpha); - return(vannulus(config, alpha, pos)); + return(vannulus(config)); break; } case (D_POLYGON): { c = pos_polygon(config, pos, &alpha); - return(vpolygon(config, alpha, pos)); + return(vpolygon(config)); break; } case (D_REULEAUX): { c = pos_reuleaux(config, pos, &alpha); - return(vreuleaux(config, alpha, pos)); + return(vreuleaux(config)); break; } case (D_FLOWER): { c = pos_flower(config, pos, &alpha); - return(vflower(config, alpha, pos)); + return(vflower(config)); break; } case (D_ALT_REU): { c = pos_alt_reuleaux(config, pos, &alpha); - return(valt_reuleaux(config, alpha, pos)); + return(valt_reuleaux(config)); break; } case (D_ANGLE): { c = pos_angle(config, pos, &alpha); - return(vangle(config, alpha, pos)); + return(vangle(config)); break; } case (D_LSHAPE): { c = pos_lshape(config, pos, &alpha); - return(vlshape(config, alpha, pos)); + return(vlshape(config)); break; } case (D_GENUSN): { c = pos_genusn(config, pos, &alpha); - return(vgenusn(config, alpha, pos)); + return(vgenusn(config)); break; } case (D_CIRCLES): { c = pos_circles(config, pos, &alpha); - return(vcircles(config, alpha, pos)); + return(vcircles(config)); break; } case (D_CIRCLES_IN_RECT): { c = pos_circles_in_rect(config, pos, &alpha); - return(vcircles_in_rect(config, alpha, pos)); + return(vcircles_in_rect(config)); + break; + } + case (D_CIRCLES_IN_GENUSN): + { + c = pos_circles_in_genusn(config, pos, &alpha); + + return(vcircles_in_genusn(config)); break; } default: @@ -3508,9 +3567,8 @@ int color[NPARTMAX]; } } - int xy_in_billiard(x, y) + int xy_in_billiard(double x, double y) /* returns 1 if (x,y) represents a point in the billiard */ - double x, y; { double l2, r1, r2, omega, omega2, c, angle, x1, y1, x2, co, so, x2plus, x2minus; int condition, k; @@ -3570,15 +3628,7 @@ int color[NPARTMAX]; } case D_POLYGON: { - condition = 1; - omega = DPI/((double)NPOLY); - c = cos(omega*0.5); - for (k=0; k circlerad[k]*circlerad[k])*circleactive[k]; - } + if (circleactive[k]) condition = condition*out_circle(x-circlex[k], y-circley[k], circlerad[k]); return(condition); break; } case D_CIRCLES_IN_RECT: { -// if ((vabs(x) = LAMBDA)||(vabs(y) >= 1.0)) return(0); else { condition = 1; for (k=0; k circlerad[k]*circlerad[k])*circleactive[k]; - } + if (circleactive[k]) condition = condition*out_circle(x-circlex[k], y-circley[k], circlerad[k]); + return(condition); + } + break; + } + case D_CIRCLES_IN_GENUSN: + { + condition = in_polygon(x, y, 1.0, NPOLY, APOLY); + if (condition == 0) return(0); + else /* test whether (x,y) outside all circles */ + { + condition = 1; + for (k=0; k YMAX) y -= height; + circlerad[n] = MU; + circleactive[n] = 1; + } + + break; + } +// case (C_LASER): +// { +// ncircles = 17; +// +// xx[0] = 0.5*(X_SHOOTER + X_TARGET); +// xx[1] = LAMBDA - 0.5*(X_TARGET - X_SHOOTER); +// xx[2] = -xx[0]; +// xx[3] = -xx[1]; +// +// yy[0] = 0.5*(Y_SHOOTER + Y_TARGET); +// yy[1] = 1.0 - 0.5*(Y_TARGET - Y_SHOOTER); +// yy[2] = -yy[0]; +// yy[3] = -yy[1]; +// +// for (i = 0; i < 4; i++) +// for (j = 0; j < 4; j++) +// { +// circlex[4*i + j] = xx[i]; +// circley[4*i + j] = yy[j]; +// +// } +// +// circlex[ncircles - 1] = X_TARGET; +// circley[ncircles - 1] = Y_TARGET; +// +// for (i=0; i YMIN - MU)) circleactive[n] = 1; + else circleactive[n] = 0; } break; } @@ -388,9 +444,9 @@ void init_circle_config() for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; - circlex[n] = ((double)(i-NGRIDX/2) + 0.5 + 0.5*(double)rand()/RAND_MAX)*dy; - circley[n] = YMIN + ((double)j + 0.5 + 0.5*(double)rand()/RAND_MAX)*dy; - circlerad[n] = MU*(1.0 + 0.35*(double)rand()/RAND_MAX); + circlex[n] = ((double)(i-NGRIDX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circley[n] = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circlerad[n] = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); circleactive[n] = 1; } break; @@ -461,6 +517,129 @@ void init_circle_config() } break; } + case (C_GOLDEN_MEAN): + { + ncircles = 300; + gamma = (sqrt(5.0) - 1.0)*0.5; /* golden mean */ + height = YMAX - YMIN; + dx = 2.0*LAMBDA/((double)ncircles); + for (n = 0; n < ncircles; n++) + { + circlex[n] = -LAMBDA + n*dx; + circley[n] = y; + y += height*gamma; + if (y > YMAX) y -= height; + circlerad[n] = MU; + circleactive[n] = 1; + } + + /* test for circles that overlap top or bottom boundary */ + ncirc0 = ncircles; + for (n=0; n < ncirc0; n++) + { + if (circley[n] + circlerad[n] > YMAX) + { + circlex[ncircles] = circlex[n]; + circley[ncircles] = circley[n] - height; + circlerad[ncircles] = MU; + circleactive[ncircles] = 1; + ncircles ++; + } + else if (circley[n] - circlerad[n] < YMIN) + { + circlex[ncircles] = circlex[n]; + circley[ncircles] = circley[n] + height; + circlerad[ncircles] = MU; + circleactive[ncircles] = 1; + ncircles ++; + } + } + break; + } + case (C_GOLDEN_SPIRAL): + { + ncircles = 1; + circlex[0] = 0.0; + circley[0] = 0.0; + + gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ + phi = 0.0; + r0 = 2.0*MU; + r = r0 + MU; + + for (i=0; i<1000; i++) + { + x = r*cos(phi); + y = r*sin(phi); + + phi += gamma; + r += MU*r0/r; + + if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU)) + { + circlex[ncircles] = x; + circley[ncircles] = y; + ncircles++; + } + } + + for (i=0; i YMIN - MU)) circleactive[i] = 1; +// printf("i = %i, circlex = %.3lg, circley = %.3lg\n", i, circlex[i], circley[i]); + } + break; + } + case (C_SQUARE_HEX): + { + ncircles = NGRIDX*(NGRIDY+1); + dy = (YMAX - YMIN)/((double)NGRIDY); + dx = dy*0.5*sqrt(3.0); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY+1; j++) + { + n = (NGRIDY+1)*i + j; + circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy; /* is +0.5 needed? */ + circley[n] = YMIN + ((double)j - 0.5)*dy; + if (((i+NGRIDX)%4 == 2)||((i+NGRIDX)%4 == 3)) circley[n] += 0.5*dy; + circlerad[n] = MU; + /* activate only circles that intersect the domain */ + if ((circley[n] < YMAX + MU)&&(circley[n] > YMIN - MU)) circleactive[n] = 1; + else circleactive[n] = 0; + } + break; + } + case (C_ONE): + { + circlex[ncircles] = 0.0; + circley[ncircles] = 0.0; + circlerad[ncircles] = MU; + circleactive[ncircles] = 1; + ncircles += 1; + break; + } + case (C_TWO): /* used for comparison with cloak */ + { + circlex[ncircles] = 0.0; + circley[ncircles] = 0.0; + circlerad[ncircles] = MU; + circleactive[ncircles] = 2; + ncircles += 1; + + circlex[ncircles] = 0.0; + circley[ncircles] = 0.0; + circlerad[ncircles] = 2.0*MU; + circleactive[ncircles] = 1; + ncircles += 1; + break; + } + case (C_NOTHING): + { + ncircles += 0; + break; + } default: { printf("Function init_circle_config not defined for this pattern \n"); @@ -477,19 +656,19 @@ double x, y; int i, j, k, k1, k2, condition; switch (B_DOMAIN) { - case D_RECTANGLE: + case (D_RECTANGLE): { if ((vabs(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); @@ -497,13 +676,13 @@ double x, y; else return(0); break; } - case D_SINAI: + case (D_SINAI): { if (x*x + y*y > LAMBDA*LAMBDA) return(1); else return(0); break; } - case D_DIAMOND: + case (D_DIAMOND): { l2 = LAMBDA*LAMBDA; r2 = l2 + (LAMBDA-1.0)*(LAMBDA-1.0); @@ -514,26 +693,26 @@ double x, y; else return(0); break; } - case D_TRIANGLE: + case (D_TRIANGLE): { if ((x>-LAMBDA)&&(y>-1.0)&&(LAMBDA*y+x<0.0)) return(1); else return(0); break; } - case D_FLAT: + case (D_FLAT): { if (y > -LAMBDA) return(1); else return(0); break; } - case D_ANNULUS: + 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: + case (D_POLYGON): { condition = 1; omega = DPI/((double)NPOLY); @@ -546,13 +725,13 @@ double x, y; // for (k=0; k MU)) return(1); else if ((vabs(y-LAMBDA) < MU)||(vabs(y+LAMBDA) < MU)) return (1); else return(0); } - case D_GRATING: + case (D_GRATING): { k1 = -(int)((-YMIN)/LAMBDA); k2 = (int)(YMAX/LAMBDA); @@ -565,11 +744,11 @@ double x, y; // printf("x = %.3lg, y = %.3lg, k1 = %i, k2 = %i, condition = %i\n", x, y, k1, k2, condition); return(condition); } - case D_EHRENFEST: + 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))); } - case D_DISK_GRID: + case (D_DISK_GRID): { dy = (YMAX - YMIN)/((double)NGRIDY); for (i = -NGRIDX/2; i < NGRIDX/2; i++) @@ -581,7 +760,7 @@ double x, y; } return(1); } - case D_DISK_HEX: + case (D_DISK_HEX): { dy = (YMAX - YMIN)/((double)NGRIDY); dx = dy*0.5*sqrt(3.0); @@ -595,7 +774,7 @@ double x, y; } return(1); } - case D_CIRCLES: + case (D_CIRCLES): { for (i = 0; i < ncircles; i++) if (circleactive[i]) @@ -607,7 +786,7 @@ double x, y; } return(1); } - case D_MENGER: + case (D_MENGER): { x1 = 0.5*(x+1.0); y1 = 0.5*(y+1.0); @@ -619,7 +798,7 @@ double x, y; } return(1); } - case D_JULIA_INT: + case (D_JULIA_INT): { u = x/JULIA_SCALE; v = y/JULIA_SCALE; @@ -634,7 +813,7 @@ double x, y; if (u*u + v*v < MANDELLIMIT) return(1); else return(0); } - case D_MENGER_ROTATED: + case (D_MENGER_ROTATED): { x2 = 1.0*(x + y); y2 = 1.0*(x - y); @@ -651,7 +830,7 @@ double x, y; } return(1); } - case D_ANNULUS_HEATED: /* returns 2 if in inner circle */ + case (D_ANNULUS_HEATED): /* returns 2 if in inner circle */ { l2 = LAMBDA*LAMBDA; r2 = x*x + y*y; @@ -660,7 +839,7 @@ double x, y; else if (r2mu <= l2) return(2); else return (0); } - case D_MENGER_HEATED: + case (D_MENGER_HEATED): { if ((vabs(x) >= 1.0)||(vabs(y) >= 1.0)) return(0); else @@ -676,7 +855,7 @@ double x, y; return(1); } } - case D_MENGER_H_OPEN: /* returns 2 if in inner circle */ + case (D_MENGER_H_OPEN): /* returns 2 if in inner circle */ { x1 = 0.5*(x+1.0); y1 = 0.5*(y+1.0); @@ -688,7 +867,7 @@ double x, y; } return(1); } - case D_MANDELBROT: + case (D_MANDELBROT): { u = 0.0; v = 0.0; @@ -707,7 +886,7 @@ double x, y; 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: + case (D_MANDELBROT_CIRCLE): { u = 0.0; v = 0.0; @@ -723,7 +902,7 @@ double x, y; else if ((x-LAMBDA)*(x-LAMBDA) + (y-0.5)*(y-0.5) < MU*MU) return(2); else return(1); } - case D_JULIA: + case (D_JULIA): { u = x/JULIA_SCALE; v = y/JULIA_SCALE; @@ -750,9 +929,8 @@ double x, y; } } -int ij_in_billiard(i, j) +int ij_in_billiard(int i, int j) /* returns 1 if (i,j) represents a point in the billiard */ -int i, j; { double xy[2]; @@ -811,31 +989,13 @@ void draw_billiard() /* draws the billiard boundary */ 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(); + + draw_circle(x0, 0.0, r, NSEG); + draw_circle(-x0, 0.0, r, NSEG); } break; } - case D_STADIUM: + case (D_STADIUM): { glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) @@ -857,21 +1017,12 @@ void draw_billiard() /* draws the billiard boundary */ glEnd(); break; } - case D_SINAI: + 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(); + draw_circle(0.0, 0.0, LAMBDA, NSEG); break; } - case D_DIAMOND: + case (D_DIAMOND): { alpha = atan(1.0 - 1.0/LAMBDA); dphi = (PID - 2.0*alpha)/(double)NSEG; @@ -936,26 +1087,8 @@ void draw_billiard() /* draws the billiard boundary */ } 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 (); + draw_circle(0.0, 0.0, LAMBDA, NSEG); + draw_circle(0.0, 0.0, 1.0, NSEG); break; } case (D_POLYGON): @@ -1008,27 +1141,18 @@ void draw_billiard() /* draws the billiard boundary */ glEnd(); break; } - case D_GRATING: + 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 (); + draw_circle(0.0, z, MU, NSEG); } break; } - case D_EHRENFEST: + case (D_EHRENFEST): { alpha = asin(MU/LAMBDA); x0 = 1.0 - sqrt(LAMBDA*LAMBDA - MU*MU); @@ -1053,7 +1177,7 @@ void draw_billiard() /* draws the billiard boundary */ glEnd (); break; } - case D_DISK_GRID: + case (D_DISK_GRID): { glLineWidth(2); for (i = -NGRIDX/2; i < NGRIDX/2; i++) @@ -1063,20 +1187,11 @@ void draw_billiard() /* draws the billiard boundary */ dx = dy*0.5*sqrt(3.0); x1 = ((double)i + 0.5)*dy; y1 = YMIN + ((double)j + 0.5)*dy; - glBegin(GL_LINE_LOOP); - for (k=0; k<=NSEG; k++) - { - phi = (double)k*DPI/(double)NSEG; - x = x1 + MU*cos(phi); - y = y1 + MU*sin(phi); - xy_to_pos(x, y, pos); - glVertex2d(pos[0], pos[1]); - } - glEnd (); + draw_circle(x1, y1, MU, NSEG); } break; } - case D_DISK_HEX: + case (D_DISK_HEX): { glLineWidth(2); for (i = -NGRIDX/2; i < NGRIDX/2; i++) @@ -1086,16 +1201,7 @@ void draw_billiard() /* draws the billiard boundary */ x1 = ((double)i + 0.5)*dy; y1 = YMIN + ((double)j + 0.5)*dy; if ((i+NGRIDX)%2 == 1) y1 += 0.5*dy; - glBegin(GL_LINE_LOOP); - for (k=0; k<=NSEG; k++) - { - phi = (double)k*DPI/(double)NSEG; - x = x1 + MU*cos(phi); - y = y1 + MU*sin(phi); - xy_to_pos(x, y, pos); - glVertex2d(pos[0], pos[1]); - } - glEnd (); + draw_circle(x1, y1, MU, NSEG); } break; } @@ -1103,19 +1209,7 @@ void draw_billiard() /* draws the billiard boundary */ { glLineWidth(BOUNDARY_WIDTH); for (i = 0; i < ncircles; i++) - if (circleactive[i]) - { - glBegin(GL_LINE_LOOP); - for (k=0; k<=NSEG; k++) - { - phi = (double)k*DPI/(double)NSEG; - x = circlex[i] + circlerad[i]*cos(phi); - y = circley[i] + circlerad[i]*sin(phi); - xy_to_pos(x, y, pos); - glVertex2d(pos[0], pos[1]); - } - glEnd (); - } + if (circleactive[i]) draw_circle(circlex[i], circley[i], circlerad[i], NSEG); break; } case (D_MENGER): @@ -1221,26 +1315,8 @@ void draw_billiard() /* draws the billiard boundary */ } 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 (); + draw_circle(MU, 0.0, LAMBDA, NSEG); + draw_circle(0.0, 0.0, 1.0, NSEG); break; } case (D_MENGER_HEATED): diff --git a/sub_wave_comp.c b/sub_wave_comp.c new file mode 100644 index 0000000..4b0cb74 --- /dev/null +++ b/sub_wave_comp.c @@ -0,0 +1,667 @@ +/* some changes in sub_wave.c required by wave_comparison.c */ + +short int circletop[NMAXCIRCLES]; /* set to 1 if circle is in top half */ + + +void init_circle_config_half(int pattern, int top) +/* initialise the arrays circlex, circley, circlerad and circleactive */ +/* for billiard shape D_CIRCLES */ +{ + int i, j, n, ncirc0; + double dx, dy, p, phi, r, ra[5], sa[5], height, y = 0.0, gamma, ymean, ytop, ybottom; + + ymean = 0.5*(YMIN + YMAX); + switch (pattern) { + case (C_SQUARE): + { + dy = (YMAX - YMIN)/((double)NGRIDY); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY/2; j++) + { + printf("i = %i, j = %i, n = %i\n", i, j, n); + n = ncircles + NGRIDY*i/2 + j; + circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy; + y = ((double)j + 0.5)*dy; + if (top) circley[n] = ymean + y; + else circley[n] = ymean - y; + circlerad[n] = MU; + circleactive[n] = 1; + circletop[n] = top; + } + ncircles += NGRIDX*NGRIDY/2; + break; + } + case (C_HEX): + { + dy = (YMAX - YMIN)/((double)NGRIDY); + dx = dy*0.5*sqrt(3.0); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY/2+1; j++) + { + n = ncircles + (NGRIDY+1)*i/2 + j; + circlex[n] = ((double)(i-NGRIDX/2))*dy; + y = ((double)j - 0.5)*dy; + if ((i+NGRIDX)%2 == 1) y += 0.5*dy; + if (top) circley[n] = ymean + 0.5*dy + y; + else circley[n] = ymean - 0.5*dy - y; + circlerad[n] = MU; + circleactive[n] = 1; + circletop[n] = top; + } + ncircles += NGRIDX*(NGRIDY+1)/2; + break; + } + case (C_RAND_DISPLACED): + { + dy = (YMAX - YMIN)/((double)NGRIDY); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY/2; j++) + { + n = ncircles + NGRIDY*i/2 + j; + circlex[n] = ((double)(i-NGRIDX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + if (NGRIDX%2 == 0) circlex[n] += 0.5*dy; + y = ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + if (top) circley[n] = ymean + y; + else circley[n] = ymean - y; + circlerad[n] = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); + circleactive[n] = 1; + circletop[n] = top; + printf("n = %i, x = %.3lg\n", n, circlex[n]); + } + ncircles += NGRIDX*NGRIDY/2; + break; + } + case (C_RAND_PERCOL): + { + dy = (YMAX - YMIN)/((double)NGRIDY); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY/2; j++) + { + n = ncircles + NGRIDY*i/2 + j; + circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy; + y = YMIN + ((double)j + 0.5)*dy; + if ( ((top)&&(y > 0.0))||((!top)&&(y <= 0.0)) ) + circley[n] = y; + circlerad[n] = MU; + p = (double)rand()/RAND_MAX; + if (p < P_PERCOL) circleactive[n] = 1; + else circleactive[n] = 0; + circletop[n] = top; + } + ncircles += NGRIDX*NGRIDY/2; + break; + } + case (C_RAND_POISSON): + { + for (i = 0; i < NPOISSON/2; i++) + { + n = ncircles + i; + circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + if (top) y = ymean + (YMAX-ymean)*(double)rand()/RAND_MAX; + else y = ymean + (YMIN-ymean)*(double)rand()/RAND_MAX; + if ( ((top)&&(y > 0.0))||((!top)&&(y <= 0.0)) ) + circley[n] = y; + circlerad[n] = MU; + circleactive[n] = 1; + circletop[n] = top; + printf("n = %i, x = %.3lg\n", n, circlex[n]); + } + ncircles += NPOISSON/2; + break; + } + case (C_CLOAK): + { + for (i = 0; i < 40; i++) + for (j = 0; j < 5; j++) + { + n = ncircles + 5*i + j; + phi = (double)i*DPI/40.0; + r = LAMBDA*0.5*(1.0 + (double)j/5.0); + circlex[n] = r*cos(phi); + y = r*sin(phi); + if ( ((top)&&(y > 0.0))||((!top)&&(y <= 0.0)) ) + circley[n] = y; + circlerad[n] = MU; + circleactive[n] = 1; + circletop[n] = top; + } + ncircles += 200; + break; + } + case (C_CLOAK_A): /* optimized model A1 by C. Jo et al */ + { + ra[0] = 0.0731; sa[0] = 1.115; + ra[1] = 0.0768; sa[1] = 1.292; + ra[2] = 0.0652; sa[2] = 1.464; + ra[3] = 0.056; sa[3] = 1.633; + ra[4] = 0.0375; sa[4] = 1.794; + for (i = 0; i < 21; i++) + for (j = 0; j < 5; j++) + { + n = ncircles + 5*i + j; + phi = (double)i*DPI/40.0; + r = LAMBDA*sa[j]; + circlex[n] = r*cos(phi); + circley[n] = r*sin(phi); + + if (top) y = r*sin(phi); + else y = -r*sin(phi); + + circley[n] = y; + + circlerad[n] = LAMBDA*ra[j]; + circleactive[n] = 1; + circletop[n] = top; + } + ncircles += 105; + + /* add circle in the center */ + circlex[ncircles] = 0.0; + circley[ncircles] = 0.0; + circlerad[ncircles] = MU; + circleactive[ncircles] = 2; + circletop[ncircles] = top; + ncircles += 1; + break; + } + case (C_GOLDEN_MEAN): + { + gamma = (sqrt(5.0) - 1.0)*0.5; /* golden mean */ + height = YMAX - ymean; + dx = 2.0*LAMBDA/150.0; + if (top) y = ymean + 0.5*height; + else y = ymean - 0.5*height; + for (n = 0; n < 150; n++) + { + circlex[ncircles + n] = -LAMBDA + n*dx; + circley[ncircles + n] = y; + + if (top) + { + y += height*gamma; + if (y > YMAX) y -= height; + } + else + { + y -= height*gamma; + if (y < YMIN) y += height; + + } + + circlerad[ncircles + n] = MU; + circleactive[ncircles + n] = 1; + circletop[ncircles] = top; + } + + /* test for circles that overlap top or bottom boundary */ + ncirc0 = ncircles; + ncircles += 150; + if (top) ytop = YMAX; + else ytop = ymean; + + if (top) ybottom = ymean; + else ybottom = YMIN; + + for (n=0; n < 150; n++) + { + if (circley[ncirc0 + n] + circlerad[ncirc0 + n] > ytop) + { + circlex[ncircles] = circlex[ncirc0 + n]; + circley[ncircles] = circley[ncirc0 + n] - height; + circlerad[ncircles] = MU; + circleactive[ncircles] = 1; + ncircles ++; + } + else if (circley[ncirc0 + n] - circlerad[ncirc0 + n] < ybottom) + { + circlex[ncircles] = circlex[ncirc0 + n]; + circley[ncircles] = circley[ncirc0 + n] + height; + circlerad[ncircles] = MU; + circleactive[ncircles] = 1; + ncircles ++; + } + } + + break; + } + case (C_ONE): + { + circlex[ncircles] = 0.0; + circley[ncircles] = 0.0; + if (top) circlerad[ncircles] = MU; + else circlerad[ncircles] = MUB; + circleactive[ncircles] = 1; + circletop[ncircles] = top; + ncircles += 1; + break; + } + case (C_TWO): + { + circlex[ncircles] = 0.0; + circley[ncircles] = 0.0; + if (top) circlerad[ncircles] = MU; + else circlerad[ncircles] = MUB; + circleactive[ncircles] = 2; + circletop[ncircles] = top; + ncircles += 1; + + circlex[ncircles] = 0.0; + circley[ncircles] = 0.0; + if (top) circlerad[ncircles] = 2.0*MU; + else circlerad[ncircles] = 2.0*MUB; + circleactive[ncircles] = 1; + circletop[ncircles] = top; + ncircles += 1; + break; + } + case (C_NOTHING): + { + ncircles += 0; + break; + } + default: + { + printf("Function init_circle_config not defined for this pattern \n"); + } + } +} + + +void init_circle_config_comp() +/* initialise the arrays circlex, circley, circlerad and circleactive */ +/* for billiard shape D_CIRCLES */ +{ + ncircles = 0; + init_circle_config_half(CIRCLE_PATTERN, 1); + init_circle_config_half(CIRCLE_PATTERN_B, 0); +} + + +int xy_in_billiard_half(double x, double y, int domain, int pattern, int top) +/* returns 1 if (x,y) represents a point in the billiard */ +{ + double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy; + int i, j, k, k1, k2, condition, type; + + switch (domain) { + case D_MENGER: + { + x1 = 0.5*(x+1.0); + y1 = 0.5*(y+1.0); + for (k=0; k 0.0)&&(vabs(x)<1.0)&&(vabs(y)<1.0)&&(((int)x1 % MRATIO)==MRATIO/2)&&(((int)y1 % MRATIO)==MRATIO/2)) return(0); + else if ((!top)&&(y < 0.0)&&(vabs(x)<1.0)&&(vabs(y)<1.0)&&(((int)x1 % MRATIO)==MRATIO/2)&&(((int)y1 % MRATIO)==MRATIO/2)) return(0); + } + return(1); + } + case D_CIRCLES: + { + for (i = 0; i < ncircles; i++) + if (circleactive[i] != 0) + { + /* choose specific type according to value of circleactive[i] */ + if (circleactive[i] == 1) type = 0; + else type = circleactive[i]; + + x1 = circlex[i]; + y1 = circley[i]; + r2 = circlerad[i]*circlerad[i]; + if ((top)&&(circletop[i])&&(y > 0.0)&&((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2)) return(type); + else if ((!top)&&(!circletop[i])&&(y < 0.0)&&((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2)) return(type); + } + return(1); + } + default: + { + printf("Function ij_in_billiard not defined for this billiard \n"); + return(0); + } + } +} + + +int xy_in_billiard_comp(double x, double y) +/* returns 1 if (x,y) represents a point in the billiard */ +{ + if (y > 0) return(xy_in_billiard_half(x, y, B_DOMAIN, CIRCLE_PATTERN, 1)); + else return(xy_in_billiard_half(x, y, B_DOMAIN_B, CIRCLE_PATTERN_B, 0)); +} + + +int ij_in_billiard_comp(int i, int j) +/* returns 1 if (i,j) represents a point in the billiard */ +{ + double xy[2]; + + ij_to_xy(i, j, xy); + + return(xy_in_billiard_comp(xy[0], xy[1])); +} + + +void draw_billiard_half(int domain, int pattern, int top) /* draws the billiard boundary */ +/* two domain version implemented for D_CIRCLES */ +{ + double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, signtop; + int i, j, k, k1, k2, mr2; + + if (BLACK) glColor3f(1.0, 1.0, 1.0); + else glColor3f(0.0, 0.0, 0.0); + glLineWidth(5); + + if (top) signtop = 1.0; + else signtop = -1.0; + + glEnable(GL_LINE_SMOOTH); + + switch (domain) { + 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, 0.0, x, signtop*x); + } + + /* level 2 */ + if (MDEPTH > 1) + { + glLineWidth(1); + mr2 = MRATIO*MRATIO; + l = 2.0/((double)mr2); + + for (i=0; i 2) + { + glLineWidth(1); + l = 2.0/((double)(mr2*MRATIO)); + + for (i=0; i 0.0)) y = 0.0; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + glEnd (); + } + 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"); + } + } +} + + +void draw_billiard_comp() /* draws the billiard boundary */ +{ + draw_billiard_half(B_DOMAIN, CIRCLE_PATTERN, 1); + draw_billiard_half(B_DOMAIN_B, CIRCLE_PATTERN_B, 0); +} + + +void int_planar_wave_comp(double x, double y, double *phi[NX], double *psi[NX], short int * xy_in[NX]) +/* initialise field with drop at (x,y) - phi is wave height, psi is phi at time t-1 */ +/* beta version, works for vertical planar wave only so far */ +{ + int i, j; + double xy[2], dist2; + + for (i=0; i NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); + else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + } + glColor3f(rgb[0], rgb[1], rgb[2]); + + glVertex2i(i, j); + glVertex2i(i+1, j); + glVertex2i(i+1, j+1); + glVertex2i(i, j+1); + } + } + + glEnd (); + + /* draw horizontal mid line */ + glColor3f(1.0, 1.0, 1.0); + glBegin(GL_LINE_STRIP); + xy_to_pos(XMIN, 0.5*(YMIN+YMAX), pos); + glVertex2d(pos[0], pos[1]); + xy_to_pos(XMAX, 0.5*(YMIN+YMAX), pos); + glVertex2d(pos[0], pos[1]); + glEnd(); +} + +void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX], double energies[6]) +/* compute total energy in top/bottom left/right boxes */ +{ + int i, j, ij[2]; + double energy = 0.0, pos, xleft = XMAX, xright = XMIN; + static short int first = 1; + static int ileft, iright, jmid = NY/2; + + if (first) /* compute box limits */ + { + /* find leftmost and rightmost circle */ + for (i=0; i xright)) xright = circlex[i] + circlerad[i]; + + xy_to_ij(xleft, 0.0, ij); + ileft = ij[0]; + xy_to_ij(xright, 0.0, ij); + iright = ij[0]; + first = 0; + + printf("xleft = %.3lg, xright = %.3lg", xleft, xright); + } + + /* top left box */ + for (i=0; i #define MOVIE 0 /* set to 1 to generate movie */ +#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ /* General geometrical parameters */ @@ -58,32 +62,28 @@ #define XMAX 2.0 /* x interval */ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -// #define XMIN -1.8 -// #define XMAX 1.8 /* x interval */ -// #define YMIN -1.0125 -// #define YMAX 1.0125 /* y interval for 9/16 aspect ratio */ #define JULIA_SCALE 1.0 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 3 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 4 /* pattern of circles, see list in global_pdes.c */ +#define CIRCLE_PATTERN 11 /* pattern of circles, see list in global_pdes.c */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define LAMBDA 0.75 /* parameter controlling the dimensions of domain */ -#define MU 0.025 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.85 /* parameter controlling the dimensions of domain */ +#define MU 0.03 /* parameter controlling the dimensions of domain */ #define NPOLY 3 /* number of sides of polygon */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ -#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ -#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ +#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ +#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 15 /* number of grid point for grid of disks */ +#define NGRIDX 15 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ /* You can add more billiard tables by adapting the functions */ @@ -91,16 +91,21 @@ /* Physical parameters of wave equation */ -#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ +#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ -#define OMEGA 0.9 /* frequency of periodic excitation */ -#define COURANT 0.01 /* Courant number */ -#define COURANTB 0.003 /* Courant number in medium B */ -#define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMA_BC 1.0e-4 /* damping factor on boundary */ -// #define GAMMA 5.0e-10 /* damping factor in wave equation */ -#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ -#define KAPPA_BC 5.0e-4 /* "elasticity" term on absorbing boundary */ +#define OMEGA 0.002 /* frequency of periodic excitation */ +#define AMPLITUDE 1.0 /* amplitude of periodic excitation */ +#define COURANT 0.02 /* Courant number */ +#define COURANTB 0.01 /* Courant number in medium B */ +#define GAMMA 0.0 /* damping factor in wave equation */ +#define GAMMAB 1.0e-6 /* damping factor in wave equation */ +#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ +#define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */ +#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ +#define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ +#define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */ /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */ /* The physical damping coefficient is given by GAMMA/(DT)^2 */ /* Increasing COURANT speeds up the simulation, but decreases accuracy */ @@ -112,10 +117,11 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 7000 /* number of frames of movie */ -#define NVID 50 /* number of iterations between images displayed on screen */ +#define NSTEPS 3000 /* number of frames of movie */ +#define NVID 20 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ -#define INITIAL_TIME 200 /* time after which to start saving frames */ +#define INITIAL_TIME 250 /* time after which to start saving frames */ +#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ @@ -124,7 +130,9 @@ /* Plot type, see list in global_pdes.c */ -#define PLOT 1 +#define PLOT 0 + +#define PLOT_B 1 /* plot type for second movie */ /* Color schemes */ @@ -133,10 +141,10 @@ #define COLOR_SCHEME 1 /* choice of color scheme, see list in global_pdes.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ -#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ -// #define SLOPE 0.5 /* sensitivity of color on wave amplitude */ +#define SLOPE 1.5 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 750.0 /* scaling factor for energy representation */ +#define E_SCALE 2000.0 /* scaling factor for energy representation */ +// #define E_SCALE 2500.0 /* scaling factor for energy representation */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ #define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ @@ -144,199 +152,34 @@ #define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ #define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -220.0 /* amplitude of variation of hue for color scheme C_HUE */ -// #define HUEMEAN 320.0 /* mean value of hue for color scheme C_HUE */ -// #define HUEAMP 100.0 /* amplitude of variation of hue for color scheme C_HUE */ + +/* 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 */ -#include "global_pdes.c" -#include "sub_wave.c" +#include "global_pdes.c" /* constants and global variables */ +#include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ +#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ + double courant2, courantb2; /* Courant parameters squared */ -void init_wave(x, y, phi, psi, xy_in) -/* initialise field with drop at (x,y) - phi is wave height, psi is phi at time t-1 */ - double x, y, *phi[NX], *psi[NX]; short int * xy_in[NX]; - -{ - int i, j; - double xy[2], dist2; - - for (i=0; i= 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); - } -} - - - - /*********************/ /* animation part */ /*********************/ -double compute_energy(phi, psi, xy_in, i, j) -double *phi[NX], *psi[NX]; -short int *xy_in[NX]; -int i, j; -{ - double velocity, energy, gradientx2, gradienty2; - int iplus, iminus, jplus, jminus; - - velocity = (phi[i][j] - psi[i][j]); - - iplus = (i+1); if (iplus == NX) iplus = NX-1; - iminus = (i-1); if (iminus == -1) iminus = 0; - jplus = (j+1); if (jplus == NY) jplus = NY-1; - jminus = (j-1); if (jminus == -1) jminus = 0; - - gradientx2 = (phi[iplus][j]-phi[i][j])*(phi[iplus][j]-phi[i][j]) - + (phi[i][j] - phi[iminus][j])*(phi[i][j] - phi[iminus][j]); - gradienty2 = (phi[i][jplus]-phi[i][j])*(phi[i][jplus]-phi[i][j]) - + (phi[i][j] - phi[i][jminus])*(phi[i][j] - phi[i][jminus]); - if (xy_in[i][j]) return(E_SCALE*E_SCALE*(velocity*velocity + 0.5*COURANT*COURANT*(gradientx2+gradienty2))); - else if (TWOSPEEDS) return(E_SCALE*E_SCALE*(velocity*velocity + 0.5*COURANTB*COURANTB*(gradientx2+gradienty2))); - else return(0); -} - -void draw_wave(phi, psi, xy_in, scale, time) -/* draw the field */ -double *phi[NX], *psi[NX], scale; -short int *xy_in[NX]; -int time; -{ - int i, j, iplus, iminus, jplus, jminus; - double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2; - static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX); - - glBegin(GL_QUADS); - -// printf("dtinverse = %.5lg\n", dtinverse); - - for (i=0; i NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); - else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); - } - glColor3f(rgb[0], rgb[1], rgb[2]); - - glVertex2i(i, j); - glVertex2i(i+1, j); - glVertex2i(i+1, j+1); - glVertex2i(i, j+1); - } - } - - glEnd (); -} - -void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in) +void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], + short int *xy_in[NX]) /* 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; + double delta, x, y, c, cc, gamma; + static long time = 0; + + time++; // c = COURANT; // cc = courant2; @@ -348,11 +191,13 @@ void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in) { c = COURANT; cc = courant2; + gamma = GAMMA; } else if (TWOSPEEDS) { c = COURANTB; - cc = courantb2; + cc = courantb2; + gamma = GAMMAB; } if ((TWOSPEEDS)||(xy_in[i][j])){ @@ -381,6 +226,14 @@ void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in) jminus = (j-1) % NY; if (jminus < 0) jminus += NY; } + + /* imposing linear wave on top and bottom by making Laplacian 1d */ + if (OSCILLATE_TOPBOT) + { + if (j == NY-1) jminus = NY-1; + else if (j == 0) jplus = 0; + } + 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]; @@ -388,42 +241,46 @@ void evolve_wave_half(phi_in, psi_in, phi_out, psi_out, xy_in) /* evolve phi */ if ((B_COND == BC_PERIODIC)||(B_COND == BC_DIRICHLET)) - phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - GAMMA*(x-y); + phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - gamma*(x-y); else if (B_COND == BC_ABSORBING) { if ((i>0)&&(i0)&&(j0)&&(i= INITIAL_TIME) save_frame(); else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); + + if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) + { + draw_wave(phi, psi, xy_in, scale, i, PLOT_B); + draw_billiard(); + glutSwapBuffers(); + save_frame_counter(NSTEPS + 21 + counter); + counter++; + } /* 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 */ @@ -561,7 +406,21 @@ void animation() if (MOVIE) { + if (DOUBLE_MOVIE) + { + draw_wave(phi, psi, xy_in, scale, i, PLOT); + draw_billiard(); + glutSwapBuffers(); + } for (i=0; i<20; i++) save_frame(); + if (DOUBLE_MOVIE) + { + draw_wave(phi, psi, xy_in, scale, i, PLOT_B); + draw_billiard(); + glutSwapBuffers(); + for (i=0; i<20; i++) save_frame_counter(NSTEPS + 21 + counter + i); + } + s = system("mv wave*.tif tif_wave/"); } for (i=0; i= 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); + } +} + + +double compute_energy(double *phi[NX], double *psi[NX], short int *xy_in[NX], int i, int j) +{ + double velocity, energy, gradientx2, gradienty2; + int iplus, iminus, jplus, jminus; + + velocity = (phi[i][j] - psi[i][j]); + + iplus = (i+1); if (iplus == NX) iplus = NX-1; + iminus = (i-1); if (iminus == -1) iminus = 0; + jplus = (j+1); if (jplus == NY) jplus = NY-1; + jminus = (j-1); if (jminus == -1) jminus = 0; + + gradientx2 = (phi[iplus][j]-phi[i][j])*(phi[iplus][j]-phi[i][j]) + + (phi[i][j] - phi[iminus][j])*(phi[i][j] - phi[iminus][j]); + gradienty2 = (phi[i][jplus]-phi[i][j])*(phi[i][jplus]-phi[i][j]) + + (phi[i][j] - phi[i][jminus])*(phi[i][j] - phi[i][jminus]); + if (xy_in[i][j]) return(E_SCALE*E_SCALE*(velocity*velocity + 0.5*COURANT*COURANT*(gradientx2+gradienty2))); + else if (TWOSPEEDS) return(E_SCALE*E_SCALE*(velocity*velocity + 0.5*COURANTB*COURANTB*(gradientx2+gradienty2))); + else return(0); +} + + +double compute_variance(double *phi[NX], double *psi[NX], short int *xy_in[NX]) +/* compute the variance of the field, to adjust color scheme */ +{ + int i, j, n = 0; + double variance = 0.0; + + for (i=1; i NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); + else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + } + glColor3f(rgb[0], rgb[1], rgb[2]); + + glVertex2i(i, j); + glVertex2i(i+1, j); + glVertex2i(i+1, j+1); + glVertex2i(i, j+1); + } + } + + glEnd (); +} + + + + + diff --git a/wave_comparison.c b/wave_comparison.c new file mode 100644 index 0000000..238014c --- /dev/null +++ b/wave_comparison.c @@ -0,0 +1,486 @@ +/*********************************************************************************/ +/* */ +/* Animation of wave equation in a planar domain */ +/* */ +/* N. Berglund, december 2012, may 2021 */ +/* */ +/* UPDATE 24/04: distinction between damping and "elasticity" parameters */ +/* UPDATE 27/04: new billiard shapes, bug in color scheme fixed */ +/* UPDATE 28/04: code made more efficient, with help of Marco Mancini */ +/* */ +/* 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 wave_billiard wave_billiard.c */ +/* -L/usr/X11R6/lib -ltiff -lm -lGL -lGLU -lX11 -lXmu -lglut -O3 -fopenmp */ +/* */ +/* OMP acceleration may be more effective after executing */ +/* export OMP_NUM_THREADS=2 in the shell before running the program */ +/* */ +/* To make a video, set MOVIE to 1 and create subfolder tif_wave */ +/* 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 +#include +#include +#include +#include +#include +#include /* Sam Leffler's libtiff library. */ +#include + +#define MOVIE 0 /* set to 1 to generate movie */ + +#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 XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.125 +#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ + +#define JULIA_SCALE 1.0 /* scaling for Julia sets */ + +/* Choice of the billiard table */ + +#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */ + +#define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_pdes.c */ +#define CIRCLE_PATTERN_B 10 /* pattern of circles, see list in global_pdes.c */ + +#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ +#define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ + +#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ +#define MU 0.03 /* parameter controlling the dimensions of domain */ +#define MUB 0.03 /* parameter controlling the dimensions of domain */ +#define NPOLY 3 /* number of sides of polygon */ +#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define MDEPTH 4 /* depth of computation of Menger gasket */ +#define MRATIO 3 /* ratio defining Menger gasket */ +#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ +#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ +#define FOCI 1 /* set to 1 to draw focal points of ellipse */ +#define NGRIDX 15 /* number of grid point for grid of disks */ +#define NGRIDY 20 /* number of grid point for grid of disks */ + +/* You can add more billiard tables by adapting the functions */ +/* xy_in_billiard and draw_billiard below */ + +/* Physical parameters of wave equation */ + +#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ +#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ + +#define OMEGA 0.0035 /* frequency of periodic excitation */ +#define AMPLITUDE 1.0 /* amplitude of periodic excitation */ +#define COURANT 0.02 /* Courant number */ +#define COURANTB 0.0075 /* Courant number in medium B */ +#define GAMMA 0.0 /* damping factor in wave equation */ +#define GAMMAB 1.0e-5 /* damping factor in wave equation */ +#define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ +#define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */ +#define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ +#define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ +#define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */ +/* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */ +/* The physical damping coefficient is given by GAMMA/(DT)^2 */ +/* Increasing COURANT speeds up the simulation, but decreases accuracy */ +/* For similar wave forms, COURANT^2*GAMMA should be kept constant */ + +/* Boundary conditions, see list in global_pdes.c */ + +#define B_COND 3 + +/* Parameters for length and speed of simulation */ + +#define NSTEPS 6000 /* number of frames of movie */ +#define NVID 25 /* number of iterations between images displayed on screen */ +#define NSEG 100 /* number of segments of boundary */ +#define INITIAL_TIME 200 /* time after which to start saving frames */ +#define COMPUTE_ENERGIES 1 /* set to 1 to compute and print energies */ +#define BOUNDARY_WIDTH 2 /* width of billiard 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 */ + +/* Plot type, see list in global_pdes.c */ + +#define PLOT 1 + +/* Color schemes */ + +#define BLACK 1 /* background */ + +#define COLOR_SCHEME 1 /* choice of color scheme, see list in global_pdes.c */ + +#define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ +#define SLOPE 1.5 /* sensitivity of color on wave amplitude */ +#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ +#define E_SCALE 2000.0 /* scaling factor for energy representation */ + +#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ +#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ +#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */ +#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ +#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */ +#define HUEAMP -220.0 /* amplitude of variation of hue for color scheme C_HUE */ + +/* 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 */ + + +#include "global_pdes.c" /* constants and global variables */ +#include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ +#include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ +#include "sub_wave_comp.c" /* some functions specific to wave_comparison */ + +double courant2, courantb2; /* Courant parameters squared */ + + +/*********************/ +/* animation part */ +/*********************/ + +void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], + short int *xy_in[NX]) +/* time step of field evolution */ +/* phi is value of field at time t, psi at time t-1 */ +{ + int i, j, iplus, iminus, jplus, jminus, jmid = NY/2; + double delta, x, y, c, cc, gamma; + static long time = 0; + + time++; + + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y) + for (i=0; i= NY) jplus -= jmid; + jminus = j-1; + if (jminus < jmid) jminus += jmid; + } + } + else if (B_COND == BC_VPER_HABS) + { + iplus = (i+1); if (iplus == NX) iplus = NX-1; + iminus = (i-1); if (iminus == -1) iminus = 0; + if (j < jmid) /* lower half */ + { + jplus = (j+1) % jmid; + jminus = (j-1) % jmid; + if (jminus < 0) jminus += jmid; + } + else /* upper half */ + { + jplus = j+1; + if (jplus >= NY) jplus -= jmid; + jminus = j-1; + if (jminus < jmid) jminus += jmid; + } + } + + /* imposing linear wave on top and bottom by making Laplacian 1d */ + if (OSCILLATE_TOPBOT) + { + if (j == NY-1) jminus = NY-1; + else if (j == 0) jplus = 0; + } + + 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_PERIODIC)||(B_COND == BC_DIRICHLET)) + phi_out[i][j] = -y + 2*x + cc*delta - KAPPA*x - gamma*(x-y); + else if ((B_COND == BC_ABSORBING)||(B_COND == BC_ABS_REFLECT)) + { + if ((i>0)&&(i0)&&(j0)&&(i 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(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double *psi_tmp[NX], short int *xy_in[NX]) +/* time step of field evolution */ +/* phi is value of field at time t, psi at time t-1 */ +{ + evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in); + evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in); +} + + + +void animation() +{ + double time, scale, energies[6], top_energy, bottom_energy; + double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[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= INITIAL_TIME) save_frame(); + else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); + + /* it seems that saving too many files too fast can cause trouble with the file system */ + /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ + if (i % PAUSE == PAUSE - 1) + { + printf("Making a short pause\n"); + sleep(PSLEEP); + s = system("mv wave*.tif tif_wave/"); + } + } + + } + + if (MOVIE) + { + for (i=0; i<20; i++) save_frame(); + s = system("mv wave*.tif tif_wave/"); + } + for (i=0; i