From 1f962fc7c86b0798aba0993ccf77fbe0580b6fab Mon Sep 17 00:00:00 2001 From: nilsberglund-orleans <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Fri, 12 Nov 2021 16:22:24 +0100 Subject: [PATCH] Add files via upload --- drop_billiard.c | 57 ++- global_particles.c | 35 +- global_pdes.c | 41 +- heat.c | 126 +----- mangrove.c | 830 +++++++++++++++++++++++----------- particle_billiard.c | 253 +++++++++-- particle_pinball.c | 56 +-- schrodinger.c | 246 ++++++++-- sub_part_billiard.c | 883 ++++++++++++++++++++++++++++++------ sub_part_pinball.c | 334 +++++++++----- sub_wave.c | 1056 +++++++++++++++++++++++++++++++++---------- sub_wave_comp.c | 220 ++++----- wave_billiard.c | 372 +++++++++++++-- wave_common.c | 61 +-- wave_comparison.c | 373 ++++++++++++++- wave_energy.c | 324 +++++++++++-- 16 files changed, 4046 insertions(+), 1221 deletions(-) diff --git a/drop_billiard.c b/drop_billiard.c index 883f7a1..cf0a7cd 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -30,14 +30,11 @@ #define MOVIE 0 /* set to 1 to generate movie */ -#define WINWIDTH 720 /* window width */ -// #define WINWIDTH 1280 /* window width */ +#define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -#define XMIN -1.125 -#define XMAX 1.125 /* x interval */ +#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 */ @@ -48,26 +45,29 @@ #define B_DOMAIN 15 /* choice of domain shape */ #define CIRCLE_PATTERN 0 /* pattern of circles */ +#define POLYLINE_PATTERN 1 /* pattern of polyline */ #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 NMAXPOLY 1000 /* total number of sides of polygonal line */ // #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 NPOISSON 500 /* number of points for Poisson C_RAND_POISSON arrangement */ #define NGOLDENSPIRAL 2000 /* max number of points for C_GOLDEN_SPIRAL arrandement */ +#define SDEPTH 1 /* Sierpinski gastket depth */ -#define LAMBDA 0.0 /* parameter controlling shape of billiard */ +#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */ // #define LAMBDA 1.124950941 /* sin(36°)/sin(31.5°) for 5-star shape with 45° angles */ // #define LAMBDA 1.445124904 /* sin(36°)/sin(24°) for 5-star shape with 60° angles */ // #define LAMBDA 3.75738973 /* sin(36°)/sin(9°) for 5-star shape with 90° angles */ // #define LAMBDA -1.73205080756888 /* -sqrt(3) for Reuleaux triangle */ // #define LAMBDA 1.73205080756888 /* sqrt(3) for triangle tiling plane */ -#define MU 1.0 /* second parameter controlling shape of billiard */ +#define MU 0.7 /* parameter controlling the dimensions of domain */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NPOLY 6 /* number of sides of polygon */ +#define NPOLY 4 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define DRAW_BILLIARD 1 /* set to 1 to draw billiard */ #define DRAW_CONSTRUCTION_LINES 1 /* set to 1 to draw additional construction lines for billiard */ @@ -75,7 +75,7 @@ #define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */ -#define NPART 10000 /* number of particles */ +#define NPART 20000 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ #define NSTEPS 5000 /* number of frames of movie */ @@ -92,7 +92,7 @@ #define LMAX 0.01 /* minimal segment length triggering resampling */ #define LPERIODIC 0.1 /* lines longer than this are not drawn (useful for Sinai billiard) */ -#define LCUT 2.0 /* controls the max size of segments not considered as being cut */ +#define LCUT 10000.0 /* controls the max size of segments not considered as being cut */ #define DMIN 0.02 /* minimal distance to boundary for triggering resampling */ #define CYCLE 0 /* set to 1 for closed curve (start in all directions) */ #define ORDER_COLORS 1 /* set to 1 if colors should be drawn in order */ @@ -101,18 +101,18 @@ #define COLOR_PALETTE 1 /* Color palette, see list in global_pdes.c */ -#define NCOLORS 6 /* number of colors */ +#define NCOLORS 14 /* number of colors */ #define COLORSHIFT 3 /* hue of initial color */ #define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ #define NSEG 100 /* number of segments of boundary */ #define BILLIARD_WIDTH 4 /* width of billiard */ #define FRONT_WIDTH 4 /* width of wave front */ -#define BLACK 1 /* set to 1 for black background */ +#define BLACK 0 /* set to 1 for black background */ #define COLOR_OUTSIDE 0 /* set to 1 for colored outside */ #define OUTER_COLOR 300.0 /* color outside billiard */ -#define PAINT_INT 1 /* set to 1 to paint interior in other color (for polygon) */ -#define PAINT_EXT 0 /* set to 1 to paint exterior of billiard */ +#define PAINT_INT 0 /* set to 1 to paint interior in other color (for polygon) */ +#define PAINT_EXT 1 /* set to 1 to paint exterior of billiard */ #define PAUSE 1000 /* number of frames after which to pause */ @@ -125,6 +125,7 @@ #define DPI 6.283185307 #define PID 1.570796327 + #include "global_particles.c" #include "sub_part_billiard.c" @@ -175,6 +176,27 @@ void init_drop_config(double x0, double y0, double angle1, double angle2, double } } +void init_partial_drop_config(int i1, int i2, double x0, double y0, double angle1, double angle2, double *configs[NPARTMAX]) +/* initialize configuration: drop at (x0,y0) */ +{ + int i; + double dalpha, alpha, pos[2]; + + while (angle2 < angle1) angle2 += DPI; + dalpha = (angle2 - angle1)/((double)(i2 - i1)); + if (i2 >= NPART) i2 = NPART; + + for (i=i1; i0)&&(i0)&&(j VMAX) phi[i][j] = VMAX; - if (newphi[i][j] < -VMAX) phi[i][j] = -VMAX; - } - } - } - } - - for (i=0; i0)&&(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_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 value of field at time t, psi at time t-1 */ +/* this version of the function has been rewritten in order to minimize the number of if-branches */ { int i, j, iplus, iminus, jplus, jminus, tb_shift; double delta, x, y, c, cc, gamma, kappa, phase, phasemin; static long time = 0; - static int init_bc = 1; - static double left_bc[NY], top_bc[NX], bot_bc[NX]; + static double tc[NX][NY], tcc[NX][NY], tgamma[NX][NY], left_bc[NY], top_bc[NX], bot_bc[NX]; + static short int first = 1, init_bc = 1; time++; @@ -271,153 +471,272 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX init_bc = 0; } - #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y,c,cc,gamma,kappa) - for (i=0; i0)&&(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; + iplus = (i+1); if (iplus == NX) iplus = NX-1; + iminus = (i-1); if (iminus == -1) iminus = 0; + + delta = phi_in[iplus][0] + phi_in[iminus][0] + phi_in[i][1] - 3.0*x; + phi_out[i][0] = -y + 2*x + tcc[i][0]*delta - KAPPA*x - tgamma[i][0]*(x-y); + break; } + case (BC_PERIODIC): + { + iplus = (i+1) % NX; + iminus = (i-1) % NX; + if (iminus < 0) iminus += NX; + + delta = phi_in[iplus][0] + phi_in[iminus][0] + phi_in[i][1] + phi_in[i][NY-1] - 4.0*x; + phi_out[i][0] = -y + 2*x + tcc[i][0]*delta - KAPPA*x - tgamma[i][0]*(x-y); + break; + } + case (BC_ABSORBING): + { + iplus = (i+1); if (iplus == NX) iplus = NX-1; + iminus = (i-1); if (iminus == -1) iminus = 0; + + delta = phi_in[iplus][0] + phi_in[iminus][0] + phi_in[i][1] - 3.0*x; + phi_out[i][0] = x - tc[i][0]*(x - phi_in[i][1]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y); + break; + } + case (BC_VPER_HABS): + { + iplus = (i+1); if (iplus == NX) iplus = NX-1; + iminus = (i-1); if (iminus == -1) iminus = 0; + + delta = phi_in[iplus][0] + phi_in[iminus][0] + phi_in[i][1] + phi_in[i][NY-1] - 4.0*x; + phi_out[i][0] = -y + 2*x + tcc[i][0]*delta - KAPPA*x - tgamma[i][0]*(x-y); + break; + } + } + psi_out[i][0] = x; + } + } + + /* add oscillating boundary condition on the left corners - NEEDED ? */ + if ((i == 0)&&(OSCILLATE_LEFT)) + { + phi_out[i][0] = AMPLITUDE*cos((double)time*OMEGA); + phi_out[i][NY-1] = AMPLITUDE*cos((double)time*OMEGA); + } + + /* for debugging purposes/if there is a risk of blow-up */ + if (FLOOR) for (i=0; 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]); } @@ -429,6 +748,7 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *phi_tmp[NX], double * evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in); } + void hash_xy_to_ij(double x, double y, int ij[2]) { static int first = 1; @@ -457,18 +777,18 @@ void hash_xy_to_ij(double x, double y, int ij[2]) } -void compute_repelling_force(int i, int j, double force[2]) +void compute_repelling_force(int i, int j, double force[2], t_mangrove* mangrove) /* compute repelling force of mangrove j on mangrove i */ { double x1, y1, x2, y2, distance, r, f; - x1 = circlex[i]; - y1 = circley[i]; - x2 = circlex[j]; - y2 = circley[j]; + x1 = mangrove[i].xc; + y1 = mangrove[i].yc; + x2 = mangrove[j].xc; + y2 = mangrove[j].yc; distance = module2(x2 - x1, y2 - y1); - r = circlerad[i] + circlerad[j]; + r = mangrove[i].radius + mangrove[j].radius; if (r <= 0.0) r = 0.001*MU; f = KREPEL/(0.001 + distance*distance); @@ -485,7 +805,7 @@ void compute_repelling_force(int i, int j, double force[2]) } -void update_hashgrid(int* mangrove_hashx, int* mangrove_hashy, int* hashgrid_number, int* hashgrid_mangroves) +void update_hashgrid(t_mangrove* mangrove, int* hashgrid_number, int* hashgrid_mangroves) { int i, j, k, n, m, ij[2], max = 0; @@ -498,7 +818,7 @@ void update_hashgrid(int* mangrove_hashx, int* mangrove_hashy, int* hashgrid_num // if (circleactive[k]) { // printf("placing circle %i\t", k); - hash_xy_to_ij(circlex[k], circley[k], ij); + hash_xy_to_ij(mangrove[k].xc, mangrove[k].yc, ij); i = ij[0]; j = ij[1]; // printf("ij = (%i, %i)\t", i, j); n = hashgrid_number[i*HASHY + j]; @@ -507,8 +827,8 @@ void update_hashgrid(int* mangrove_hashx, int* mangrove_hashy, int* hashgrid_num if (m < HASHX*HASHY*HASHMAX) hashgrid_mangroves[m] = k; else printf("Too many mangroves in hash cell, try increasing HASHMAX\n"); hashgrid_number[i*HASHY + j]++; - mangrove_hashx[k] = i; - mangrove_hashy[k] = j; + mangrove[k].hashx = i; + mangrove[k].hashy = j; if (n > max) max = n; // printf("Placed mangrove %i at (%i,%i) in hashgrid\n", k, ij[0], ij[1]); @@ -528,10 +848,9 @@ void animation() int i, j, k, n, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q; static int imin, imax; static short int first = 1; - - double *circleenergy, *circley_wrapped, *anchor_x, *anchor_y, *vx, *vy, *circlerad_initial, *mass_inverse; - short int *circle_attached; - int *mangrove_hashx, *mangrove_hashy, *hashgrid_number, *hashgrid_mangroves; + t_mangrove *mangrove; + int *hashgrid_number, *hashgrid_mangroves; + t_hashgrid *hashgrid; /* Since NX and NY are big, it seemed wiser to use some memory allocation here */ for (i=0; i= YMAX) y -= circlerad[i]; - if (y <= YMIN) y += circlerad[i]; + /* to avoid having to recode init_circle_config, would be more elegant in C++ */ + mangrove[i].xc = circles[i].xc; + mangrove[i].yc = circles[i].yc; + mangrove[i].radius = circles[i].radius; + mangrove[i].active = circles[i].active; + + mangrove[i].energy = 0.0; + y = mangrove[i].yc; + if (y >= YMAX) y -= mangrove[i].radius; + if (y <= YMIN) y += mangrove[i].radius; // if (y >= YMAX) y -= (YMAX - YMIN); // if (y <= YMIN) y += (YMAX - YMIN); - circley_wrapped[i] = y; -// circleactive[i] = 1; + mangrove[i].yc_wrapped = y; +// mangrove[i].active = 1; - if (RANDOM_RADIUS) circlerad[i] = circlerad[i]*(0.75 + 0.5*((double)rand()/RAND_MAX)); + if (RANDOM_RADIUS) mangrove[i].radius = mangrove[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); - circlerad_initial[i] = circlerad[i]; - circle_attached[i] = 1; - mass_inverse[i] = MU*MU/(MANGROVE_MASS*circlerad[i]*circlerad[i]); + mangrove[i].radius_initial = mangrove[i].radius; + mangrove[i].attached = 1; + mangrove[i].mass_inv = MU*MU/(MANGROVE_MASS*mangrove[i].radius*mangrove[i].radius); if (MOVE_MANGROVES) { - anchor_x[i] = circlex[i]; - anchor_y[i] = circley_wrapped[i]; -// anchor_y[i] = circley[i]; + mangrove[i].anchorx = mangrove[i].xc; + mangrove[i].anchory = mangrove[i].yc_wrapped; +// mangrove[i].anchory = mangrove[i].yc; } if (INERTIA) { - vx[i] = 0.0; - vy[i] = 0.0; + mangrove[i].vx = 0.0; + mangrove[i].vy = 0.0; } } /* initialise hash table for interacting mangroves */ - if (REPELL_MANGROVES) update_hashgrid(mangrove_hashx, mangrove_hashy, hashgrid_number, hashgrid_mangroves); + if (REPELL_MANGROVES) update_hashgrid(mangrove, hashgrid_number, hashgrid_mangroves); if (first) /* compute box limits where circles are reset */ { /* find leftmost and rightmost circle */ for (i=0; i xright)) xright = circlex[i] + circlerad[i]; + if ((mangrove[i].active)&&(mangrove[i].xc + mangrove[i].radius > xright)) xright = mangrove[i].xc + mangrove[i].radius; xy_to_ij(xleft, 0.0, ij); imin = ij[0] - 10; @@ -670,34 +978,34 @@ void animation() /* move mangroves */ - if (MOVE_MANGROVES) for (j=0; j ", j, circlex[j], circley[j]); +// if (j%NGRIDY == 0) printf("circle %i (%.3lg, %.3lg) -> ", j, mangrove[j].xc, mangrove[j].yc); /* compute force of wave */ dx = DT_MANGROVE*KWAVE*gradient[0]; dy = DT_MANGROVE*KWAVE*gradient[1]; /* compute force of spring */ - if (circle_attached[j]) + if (mangrove[j].attached) { - dx += DT_MANGROVE*(-KSPRING*(circlex[j] - anchor_x[j])); - dy += DT_MANGROVE*(-KSPRING*(circley_wrapped[j] - anchor_y[j])); + dx += DT_MANGROVE*(-KSPRING*(mangrove[j].xc - mangrove[j].anchorx)); + dy += DT_MANGROVE*(-KSPRING*(mangrove[j].yc_wrapped - mangrove[j].anchory)); } /* compute repelling force from other mangroves */ if (REPELL_MANGROVES) { /* determine neighboring grid points */ - i0 = mangrove_hashx[j]; + i0 = mangrove[j].hashx; iminus = i0 - 1; if (iminus < 0) iminus = 0; iplus = i0 + 1; if (iplus >= HASHX) iplus = HASHX-1; - j0 = mangrove_hashy[j]; + j0 = mangrove[j].hashy; jminus = j0 - 1; if (jminus < 0) jminus = 0; jplus = j0 + 1; if (jplus >= HASHY) jplus = HASHY-1; @@ -706,9 +1014,9 @@ void animation() for (p=iminus; p<= iplus; p++) for (q=jminus; q<= jplus; q++) for (k=0; k L_DETACH) circle_attached[j] = 0; - if (length*mass_inverse[j] > L_DETACH) circle_attached[j] = 0; +// if (length > L_DETACH) mangrove[j].attached = 0; + if (length*mangrove[j].mass_inv > L_DETACH) mangrove[j].attached = 0; } if (dx > DXMAX) dx = DXMAX; @@ -735,36 +1043,36 @@ void animation() if (INERTIA) { - vx[j] += (dx - DAMP_MANGROVE*vx[j])*mass_inverse[j]; - vy[j] += (dy - DAMP_MANGROVE*vy[j])*mass_inverse[j]; - circlex[j] += vx[j]*DT_MANGROVE; - circley[j] += vy[j]*DT_MANGROVE; - circley_wrapped[j] += vy[j]*DT_MANGROVE; + mangrove[j].vx += (dx - DAMP_MANGROVE*mangrove[j].vx)*mangrove[j].mass_inv; + mangrove[j].vy += (dy - DAMP_MANGROVE*mangrove[j].vy)*mangrove[j].mass_inv; + mangrove[j].xc += mangrove[j].vx*DT_MANGROVE; + mangrove[j].yc += mangrove[j].vy*DT_MANGROVE; + mangrove[j].yc_wrapped += mangrove[j].vy*DT_MANGROVE; // if (j%NGRIDY == 0) // printf("circle %.i: (dx,dy) = (%.3lg,%.3lg), (vx,vy) = (%.3lg,%.3lg)\n", -// j, circlex[j]-anchor_x[j], circley[j]-anchor_y[j], vx[j], vy[j]); +// j, mangrove[j].xc-mangrove[j].anchorx, mangrove[j].yc-mangrove[j].anchory, mangrove[j].vx, mangrove[j].vy); } else { - circlex[j] += dx*mass_inverse[j]*DT_MANGROVE; - circley[j] += dy*mass_inverse[j]*DT_MANGROVE; - circley_wrapped[j] += dy*mass_inverse[j]*DT_MANGROVE; + mangrove[j].xc += dx*mangrove[j].mass_inv*DT_MANGROVE; + mangrove[j].yc += dy*mangrove[j].mass_inv*DT_MANGROVE; + mangrove[j].yc_wrapped += dy*mangrove[j].mass_inv*DT_MANGROVE; } - if (circlex[j] <= XMIN) circlex[j] = XMIN; - if (circlex[j] >= XMAX) circlex[j] = XMAX; - if (circley_wrapped[j] <= YMIN) circley_wrapped[j] = YMIN; - if (circley_wrapped[j] >= YMAX) circley_wrapped[j] = YMAX; + if (mangrove[j].xc <= XMIN) mangrove[j].xc = XMIN; + if (mangrove[j].xc >= XMAX) mangrove[j].xc = XMAX; + if (mangrove[j].yc_wrapped <= YMIN) mangrove[j].yc_wrapped = YMIN; + if (mangrove[j].yc_wrapped >= YMAX) mangrove[j].yc_wrapped = YMAX; -// if (j%NGRIDY == 0) printf("(%.3lg, %.3lg)\n", circlex[j], circley[j]); - +// if (j%NGRIDY == 0) printf("(%.3lg, %.3lg)\n", mangrove[j].xc, mangrove[j].yc); + redraw = 1; } /* test for debugging */ if (1) for (j=0; j 0.1*MANGROVE_EMAX) @@ -773,10 +1081,10 @@ void animation() printf("Flooring dissipation!\n"); } - if (circleactive[j]) + if (mangrove[j].active) { - circleenergy[j] += dissip; - ej = circleenergy[j]; + mangrove[j].energy += dissip; + ej = mangrove[j].energy; // printf("ej = %.3f\n", ej); if (ej <= MANGROVE_EMAX) { @@ -788,35 +1096,35 @@ void animation() else hue = MANGROVE_HUE_MIN; hsl_to_rgb(hue, 0.9, 0.5, rgb); // if (j%NGRIDY == 0) printf("Circle %i, energy %.5lg, hue %.5lg\n", j, ej, hue); - draw_colored_circle(circlex[j], circley[j], circlerad[j], NSEG, rgb); + draw_colored_circle(mangrove[j].xc, mangrove[j].yc, mangrove[j].radius, NSEG, rgb); /* shrink mangrove */ if ((ERODE_MANGROVES)&&(ej > 0.0)) { - circlerad[j] = circlerad_initial[j]*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX)); + mangrove[j].radius = mangrove[j].radius_initial*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX)); redraw = 1; } - else circlerad[j] = circlerad_initial[j]; + else mangrove[j].radius = mangrove[j].radius_initial; } else /* remove mangrove */ { - circleactive[j] = 0; + mangrove[j].active = 0; /* reinitialize table xy_in */ redraw = 1; } } else if (RECOVER_MANGROVES) /* allow disabled mangroves to recover */ { - circleenergy[j] -= 0.15*dissip; - printf("Circle %i energy %.3lg\n", j, circleenergy[j]); - if (circleenergy[j] < 0.0) + mangrove[j].energy -= 0.15*dissip; + printf("Circle %i energy %.3lg\n", j, mangrove[j].energy); + if (mangrove[j].energy < 0.0) { printf("Reactivating circle %i?\n", j); /* THE PROBLEM occurs when circleactive[0] is set to 1 again */ - if (j>0) circleactive[j] = 1; - circlerad[j] = circlerad_initial[j]; - circleenergy[j] = -MANGROVE_EMAX; + if (j>0) mangrove[j].active = 1; + mangrove[j].radius = mangrove[j].radius_initial; + mangrove[j].energy = -MANGROVE_EMAX; /* reinitialize table xy_in */ redraw = 1; } @@ -824,11 +1132,19 @@ void animation() } } + /* for compatibility with draw_billiard, may be improvable */ + for (j=0; j 0.0) // { -// circlerad[j] -= MU*ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX); -// if (circlerad[j] < 0.0) circlerad[j] = 0.0; -// circlerad[j] = circlerad_initial[j]*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX)); +// mangrove[j].radius -= MU*ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX); +// if (mangrove[j].radius < 0.0) mangrove[j].radius = 0.0; +// mangrove[j].radius = mangrove[j].radius_initial*(1.0 - ej*ej/(MANGROVE_EMAX*MANGROVE_EMAX)); // redraw = 1; // } -// else circlerad[j] = circlerad_initial[j]; +// else mangrove[j].radius = mangrove[j].radius_initial; // } // else /* remove mangrove */ // { -// circleactive[j] = 0; +// mangrove[j].active = 0; /* reinitialize table xy_in */ // redraw = 1; // } // } // else /* allow disabled mangroves to recover */ // { -// circleenergy[j] -= 0.15*dissip; -// printf("ej = %.3f\n", circleenergy[j]); -// circlerad[j] += 0.005*MU; -// if (circlerad[j] > MU) circlerad[j] = MU; -// if ((circleenergy[j] < 0.0)&&(circlerad[j] > 0.0)) -// if (circleenergy[j] < 0.0) +// mangrove[j].energy -= 0.15*dissip; +// printf("ej = %.3f\n", mangrove[j].energy); +// mangrove[j].radius += 0.005*MU; +// if (mangrove[j].radius > MU) mangrove[j].radius = MU; +// if ((mangrove[j].energy < 0.0)&&(mangrove[j].radius > 0.0)) +// if (mangrove[j].energy < 0.0) // { -// circleactive[j] = 1; -// circlerad[j] = circlerad[j]*(0.75 + 0.5*((double)rand()/RAND_MAX)); -// circlerad[j] = circlerad_initial[j]; -// circleenergy[j] = -MANGROVE_EMAX; +// mangrove[j].active = 1; +// mangrove[j].radius = mangrove[j].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); +// mangrove[j].radius = mangrove[j].radius_initial; +// mangrove[j].energy = -MANGROVE_EMAX; /* reinitialize table xy_in */ // redraw = 1; // } // } -// printf("Circle %i, energy %.5lg\n", j, circleenergy[j]); +// printf("Circle %i, energy %.5lg\n", j, mangrove[j].energy); // } printf("Updating hashgrid\n"); - if (REPELL_MANGROVES) update_hashgrid(mangrove_hashx, mangrove_hashy, hashgrid_number, hashgrid_mangroves); + if (REPELL_MANGROVES) update_hashgrid(mangrove, hashgrid_number, hashgrid_mangroves); printf("Drawing billiard\n"); @@ -941,22 +1257,10 @@ void animation() free(psi_tmp[i]); free(xy_in[i]); } - - free(circleenergy); - free(circley_wrapped); - free(anchor_x); - free(anchor_y); - free(vx); - free(vy); - free(circlerad_initial); - free(mass_inverse); - free(circle_attached); - - free(mangrove_hashx); - free(mangrove_hashy); + free(mangrove); + free(hashgrid_number); free(hashgrid_mangroves); - } diff --git a/particle_billiard.c b/particle_billiard.c index 5a183c3..a4cb77a 100644 --- a/particle_billiard.c +++ b/particle_billiard.c @@ -32,43 +32,41 @@ #define MOVIE 0 /* set to 1 to generate movie */ -// #define WINWIDTH 1280 /* window width */ -#define WINWIDTH 720 /* window width */ +#define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ -#define XMIN -1.4 -#define XMAX 1.4 /* x interval */ -#define YMIN -1.4 -#define YMAX 1.4 /* y interval for 9/16 aspect ratio */ -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -// #define YMIN -1.125 -// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.125 +#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ #define SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */ /* Choice of the billiard table, see global_particles.c */ -#define B_DOMAIN 16 /* choice of domain shape */ +#define B_DOMAIN 30 /* choice of domain shape */ -#define CIRCLE_PATTERN 2 /* pattern of circles */ +#define CIRCLE_PATTERN 1 /* pattern of circles */ +#define POLYLINE_PATTERN 1 /* pattern of polyline */ -#define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ +#define ABSORBING_CIRCLES 1 /* 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 NMAXCIRCLES 100000 /* total number of circles (must be at least NCX*NCY for square grid) */ +#define NMAXPOLY 100000 /* total number of sides of polygonal line */ // #define NCX 10 /* number of circles in x direction */ // #define NCY 10 /* number of circles in y direction */ #define NCX 30 /* number of circles in x direction */ #define NCY 20 /* number of circles in y direction */ #define NPOISSON 500 /* number of points for Poisson C_RAND_POISSON arrangement */ #define NGOLDENSPIRAL 2000 /* max number of points for C_GOLDEN_SPIRAL arrandement */ +#define SDEPTH 1 /* Sierpinski gastket depth */ -// #define LAMBDA 1.4 /* parameter controlling shape of domain */ -// #define MU 0.2 /* second parameter controlling shape of billiard */ -#define LAMBDA 1.3 /* parameter controlling shape of domain */ -#define MU 0.3 /* second parameter controlling shape of billiard */ +#define LAMBDA 1.8 /* parameter controlling shape of domain */ +#define MU 0.01 /* second parameter controlling shape of billiard */ +// #define LAMBDA 0.3 /* parameter controlling shape of domain */ +// #define MU 0.7 /* second parameter controlling shape of billiard */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NPOLY 4 /* number of sides of polygon */ +#define NPOLY 6 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define DRAW_BILLIARD 1 /* set to 1 to draw billiard */ #define DRAW_CONSTRUCTION_LINES 0 /* set to 1 to draw additional construction lines for billiard */ @@ -79,17 +77,20 @@ /* Simulation parameters */ -#define NPART 100 /* number of particles */ +#define NPART 5000 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ #define LMAX 0.01 /* minimal segment length triggering resampling */ #define DMIN 0.02 /* minimal distance to boundary for triggering resampling */ #define CYCLE 1 /* set to 1 for closed curve (start in all directions) */ -#define SHOWTRAILS 1 /* set to 1 to keep trails of the particles */ +#define SHOWTRAILS 0 /* set to 1 to keep trails of the particles */ +#define SHOWZOOM 1 /* set to 1 to show zoom on specific area */ +#define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print number of particles */ #define TEST_ACTIVE 1 /* set to 1 to test whether particle is in billiard */ -#define NSTEPS 1300 /* number of frames of movie */ -#define TIME 2500 /* time between movie frames, for fluidity of real-time simulation */ -#define DPHI 0.00001 /* integration step */ +#define NSTEPS 1000 /* number of frames of movie */ +#define TIME 1500 /* time between movie frames, for fluidity of real-time simulation */ +// #define DPHI 0.00001 /* integration step */ +#define DPHI 0.000005 /* integration step */ #define NVID 150 /* number of iterations between images displayed on screen */ /* Decreasing TIME accelerates the animation and the movie */ @@ -100,15 +101,15 @@ /* Colors and other graphical parameters */ -#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ -#define NCOLORS 64 /* number of colors */ +#define NCOLORS 16 /* number of colors */ #define COLORSHIFT 0 /* hue of initial color */ #define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */ #define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */ #define NSEG 100 /* number of segments of boundary */ -#define LENGTH 0.03 /* length of velocity vectors */ -#define BILLIARD_WIDTH 2 /* width of billiard */ +#define LENGTH 0.025 /* length of velocity vectors */ +#define BILLIARD_WIDTH 3 /* width of billiard */ #define PARTICLE_WIDTH 2 /* width of particles */ #define FRONT_WIDTH 3 /* width of wave front */ @@ -116,7 +117,7 @@ #define COLOR_OUTSIDE 0 /* set to 1 for colored outside */ #define OUTER_COLOR 270.0 /* color outside billiard */ #define PAINT_INT 0 /* set to 1 to paint interior in other color (for polygon/Reuleaux) */ -#define PAINT_EXT 1 /* set to 1 to paint exterior */ +#define PAINT_EXT 0 /* set to 1 to paint exterior */ #define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ @@ -249,6 +250,131 @@ void init_line_config(double x0, double y0, double x1, double y1, double angle, } } +void draw_zoom(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTMAX], double x_target, double y_target, double width) +/* draw zoom around target (for laser in room of mirrors) */ +{ + int i; + double x1, y1, x2, y2, xb, yb, cosphi, sinphi, rgb[3], shiftx = 0.0, shifty = 0.65, tradius, phi, zoomwidth = 0.4; + + glEnable(GL_LINE_SMOOTH); + glColor3f(1.0, 1.0, 1.0); + + /* draw zoomed area */ + glLineWidth(BILLIARD_WIDTH/2); + + x1 = x_target - width; + y1 = y_target - width; + x2 = x_target + width; + y2 = y_target + width; + glBegin(GL_LINE_LOOP); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); + + /* draw zoom boundary */ + glLineWidth(BILLIARD_WIDTH*2); + + x1 = shiftx - zoomwidth; + y1 = shifty - zoomwidth; + x2 = shiftx + zoomwidth; + y2 = shifty + zoomwidth; + glBegin(GL_LINE_LOOP); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); + + /* draw billiard boundaries in zoom */ + glLineWidth(BILLIARD_WIDTH*2); + + if (y_target + width > 1.0) + { + yb = shifty + 0.5*(1.0 - y_target)/width; + glBegin(GL_LINE_STRIP); + glVertex2d(x1, yb); + glVertex2d(x2, yb); + glVertex2d(x2, yb + 0.02); + glVertex2d(x1, yb + 0.02); + glEnd(); + } + /* other boundaries not yet implemented */ + + /* draw target in zoom */ + glLineWidth(BILLIARD_WIDTH*2); + + glColor3f(0.0, 0.8, 0.2); + glBegin(GL_LINE_LOOP); + tradius = zoomwidth*MU/width; + for (i=0; i<=NSEG; i++) + { + phi = (double)i*DPI/(double)NSEG; + x1 = shiftx + tradius*cos(phi); + y1 = shifty + tradius*sin(phi); + glVertex2d(x1, y1); + } + glEnd (); + +// glLineWidth(PARTICLE_WIDTH*2); + + for (i=0; i 1.0)) + { + if (x1 > 0.0) xb = 1.0; + else xb = -1.0; + y2 = y1 + (xb - x1)*(y2 - y1)/(x2 - x1); + x2 = xb; + } + else + if ((vabs(x1) > 1.0)&&(vabs(x2) < 1.0)) + { + if (x2 > 0.0) xb = 1.0; + else xb = -1.0; + y1 = y2 + (xb - x2)*(y1 - y2)/(x1 - x2); + x1 = xb; + } + if ((vabs(y1) < 1.0)&&(vabs(y2) > 1.0)) + { + if (y1 > 0.0) yb = 1.0; + else yb = -1.0; + x2 = x1 + (yb - y1)*(x2 - x1)/(y2 - y1); + y2 = yb; + } + else + if ((vabs(y1) > 1.0)&&(vabs(y2) < 1.0)) + { + if (y2 > 0.0) yb = 1.0; + else yb = -1.0; + x1 = x2 + (yb - y2)*(x1 - x2)/(y1 - y2); + y1 = yb; + } + +// if ((active[i])&&(vabs(x1) < 1.0)&&(vabs(y1) < 1.0)&&(vabs(x2) < 1.0)&&(vabs(y2) < 1.0)) + if (((active[i])&&(vabs(x1) < 1.0)&&(vabs(y1) < 1.0))||((vabs(x2) < 1.0)&&(vabs(y2) < 1.0))) + { + rgb_color_scheme(color[i], rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); + + glBegin(GL_LINE_STRIP); + glVertex2d(shiftx + zoomwidth*SCALING_FACTOR*x1, shifty + zoomwidth*SCALING_FACTOR*y1); + glVertex2d(shiftx + zoomwidth*SCALING_FACTOR*x2, shifty + zoomwidth*SCALING_FACTOR*y2); + glEnd (); + } + + } +} + void draw_config_showtrails(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTMAX]) /* draw the particles */ @@ -312,13 +438,15 @@ void draw_config_showtrails(int color[NPARTMAX], double *configs[NPARTMAX], int // } } if (DRAW_BILLIARD) draw_billiard(); + + if (SHOWZOOM) draw_zoom(color, configs, active, 0.95, 0.0, 0.1); } void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTMAX]) /* draw the particles */ { - int i; + int i, c; double x1, y1, x2, y2, cosphi, sinphi, rgb[3]; glutSwapBuffers(); @@ -333,12 +461,13 @@ void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPAR { if (configs[i][2]<0.0) { - vbilliard(configs[i]); + c = vbilliard(configs[i]); if (!RAINBOW_COLOR) { color[i]++; if (color[i] >= NCOLORS) color[i] -= NCOLORS; } + if ((ABSORBING_CIRCLES)&&(c < 0)) active[i] = 0; } configs[i][2] += DPHI; @@ -414,6 +543,8 @@ void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPAR if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3]; } if (DRAW_BILLIARD) draw_billiard(); + + if (SHOWZOOM) draw_zoom(color, configs, active, 0.95, 0.0, 0.1); } @@ -430,6 +561,7 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a { // printf("reflecting particle %i\n", i); c = vbilliard(configs[i]); + if ((ABSORBING_CIRCLES)&&(c < 0)) active[i] = 0; // if (c>=0) color[i]++; if ((!RAINBOW_COLOR)&&(c>=0)) color[i]++; if (!RAINBOW_COLOR) @@ -451,31 +583,66 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a // draw_config(color, configs); } +void print_part_number(double *configs[NPARTMAX], int active[NPARTMAX], double x, double y) +{ + char message[50]; + int i, n_active_particles = 0; + double rgb[3]; + + /* count active particles, using the fact that absorbed particles have been given dummy coordinates */ + for (i=0; i= 0) ncol++; - if ((c >= 0)&&(circlecolor[c] == 0)) nobst++; - circlecolor[c]++; + if ((c >= 0)&&(circles[c].color == 0)) nobst++; + circles[c].color++; /* take care of circles doubled because of periodic boundary conditions */ - if ((circleactive[c])&&(B_DOMAIN == D_CIRCLES_IN_TORUS)&&(partner_circle[c] != c)) circlecolor[partner_circle[c]]++; + if ((circles[c].active)&&(B_DOMAIN == D_CIRCLES_IN_TORUS)&&(circles[c].partner != c)) circles[circles[c].partner].color++; - newcircle[c] = 10; + circles[c].new = 10; /* update free path statistics */ if (ncol > 1) /* disregard very first collision */ @@ -519,9 +516,9 @@ void print_particle_numbers(double *configs[NPARTMAX]) { /* find leftmost and rightmost circle */ for (i=0; i xright)) xright = circlex[i] + circlerad[i]; + if ((circles[i].active)&&(circles[i].xc + circles[i].radius > xright)) xright = circles[i].xc + circles[i].radius; first = 0; @@ -574,9 +571,9 @@ void draw_statistics() /* histogram of number of collisions per peg */ for (i=0; i 0.99) circleactive[i] = 0; - if (vabs(circlex[i]) + circlerad[i] > 0.99*LAMBDA) circleactive[i] = 0; + if (vabs(circles[i].xc) + circles[i].radius > 0.99) circles[i].active = 0; + if (vabs(circles[i].xc) + circles[i].radius > 0.99*LAMBDA) circles[i].active = 0; } else if (B_DOMAIN == D_CIRCLES_IN_TORUS) for (i=0; i 1.0) circleactive[i] = 0; - if (vabs(circlex[i]) - circlerad[i] > LAMBDA) circleactive[i] = 0; + if (vabs(circles[i].yc) - circles[i].radius > 1.0) circles[i].active = 0; + if (vabs(circles[i].xc) - circles[i].radius > LAMBDA) circles[i].active = 0; } // for (i=0; i nmaxpeg) nmaxpeg = circlecolor[j]; + if (circles[j].color > nmaxpeg) nmaxpeg = circles[j].color; // sprintf(message, "max hits per peg: %d", nmaxpeg); // glColor3f(1.0, 1.0, 1.0); diff --git a/schrodinger.c b/schrodinger.c index c4cdecc..03b3252 100644 --- a/schrodinger.c +++ b/schrodinger.c @@ -39,42 +39,39 @@ /* General geometrical parameters */ -// #define WINWIDTH 1280 /* window width */ -#define WINWIDTH 720 /* window width */ +#define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ // #define NX 1280 /* number of grid points on x axis */ -#define NX 720 /* number of grid points on x axis */ -#define NY 720 /* number of grid points on y axis */ -// #define NX 640 /* number of grid points on x axis */ -// #define NY 360 /* number of grid points on y axis */ +// #define NX 720 /* number of grid points on x axis */ +#define NX 640 /* number of grid points on x axis */ +#define NY 360 /* number of grid points on y axis */ /* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */ /* but will multiply run time by 4 */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ -#define YMIN -2.0 -#define YMAX 2.0 /* y interval for 9/16 aspect ratio */ -// #define YMIN -1.125 -// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +#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, see list in global_pdes.c */ -#define B_DOMAIN 19 /* choice of domain shape */ +#define B_DOMAIN 10 /* choice of domain shape */ #define CIRCLE_PATTERN 0 /* 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 RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.0 /* parameter controlling the dimensions of domain */ -#define MU 1.75 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.1 /* parameter controlling the dimensions of domain */ +#define MU 0.03 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ -#define MDEPTH 3 /* depth of computation of Menger gasket */ +#define MDEPTH 5 /* 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 */ @@ -98,8 +95,8 @@ /* Physical patameters of wave equation */ -// #define DT 0.00000005 #define DT 0.00000001 +// #define DT 0.00000001 // #define DT 0.000000005 // #define DT 0.000000005 #define HBAR 1.0 @@ -110,9 +107,9 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1400 /* number of frames of movie */ -#define NVID 2000 /* number of iterations between images displayed on screen */ -// #define NVID 1200 /* number of iterations between images displayed on screen */ +#define NSTEPS 2500 /* number of frames of movie */ +// #define NVID 2000 /* number of iterations between images displayed on screen */ +#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 */ @@ -120,6 +117,7 @@ #define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ +#define END_FRAMES 100 /* still frames at end of movie */ /* For debugging purposes only */ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ @@ -133,22 +131,28 @@ /* Color schemes, see list in global_pdes.c */ -#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ -#define COLOR_SCHEME 1 /* choice of color scheme */ +#define COLOR_SCHEME 3 /* choice of color scheme */ #define SCALE 1 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ +#define E_SCALE 150.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 150.0 /* mean value of hue for color scheme C_HUE */ -#define HUEAMP -150.0 /* amplitude of variation of hue for color scheme C_HUE */ +#define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */ +#define HUEAMP 180.0 /* amplitude of variation of hue for color scheme C_HUE */ + +#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #include "global_pdes.c" #include "sub_wave.c" @@ -249,7 +253,7 @@ void draw_wave(double *phi[NX], double *psi[NX], short int *xy_in[NX], double sc glEnd (); } -void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], +void evolve_wave_half_old(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 */ @@ -342,6 +346,168 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX // printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]); } +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 */ +{ + int i, j, iplus, iminus, jplus, jminus; + double delta1, delta2, x, y; + + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta1,delta2,x,y) + for (i=1; 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; + } + } + } +} + 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 real part, psi is imaginary part */ @@ -392,6 +558,13 @@ void renormalise_field(double *phi[NX], double *psi[NX], short int *xy_in[NX], d } +void draw_color_bar(int plot, double range) +{ + if (ROTATE_COLOR_SCHEME) draw_color_scheme(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range); + else draw_color_scheme(1.7, YMIN + 0.1, 1.9, YMAX - 0.1, plot, -range, range); +} + + void animation() { double time, scale, dx, var; @@ -408,6 +581,10 @@ void animation() psi_tmp[i] = (double *)malloc(NY*sizeof(double)); xy_in[i] = (short int *)malloc(NY*sizeof(short int)); } + + /* initialise polyline for von Koch and simular domains */ + npolyline = init_polyline(MDEPTH, polyline); +// for (i=0; i= 1) + if (circles[k].new >= 1) { - rgb_color_scheme_lum(circlecolor[k], 0.85, rgb1); - newcircle[k]--; + rgb_color_scheme_lum(circles[k].color, 0.85, rgb1); + circles[k].new--; } - else rgb_color_scheme(circlecolor[k], rgb1); - draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb1); + else rgb_color_scheme(circles[k].color, rgb1); + draw_colored_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG, rgb1); } } init_billiard_color(); - for (k=0; k= 1) +// { +// rgb_color_scheme_lum(circlecolor[k], 0.85, rgb1); +// newcircle[k]--; +// } +// else rgb_color_scheme(circlecolor[k], rgb1); +// draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb1); +// } +// } +// init_billiard_color(); +// for (k=0; k= 1) + if (circles[k].new >= 1) { - rgb_color_scheme_lum(circlecolor[k], 0.85, rgb1); - newcircle[k]--; + rgb_color_scheme_lum(circles[k].color, 0.85, rgb1); + circles[k].new--; } - else rgb_color_scheme(circlecolor[k], rgb1); - draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb1); + else rgb_color_scheme(circles[k].color, rgb1); + draw_colored_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG, rgb1); } } init_billiard_color(); - for (k=0; k= 1) + for (k=0; k= 1) { - if (circleactive[k] == 2) + if (circles[k].active == 2) { // hsl_to_rgb(150.0, 0.9, 0.4, rgb); // glColor3f(rgb[0], rgb[1], rgb[2]); glColor3f(0.0, 1.0, 0.0); rgb[0] = 0.0; rgb[1] = 0.9; rgb[2] = 0.0; - draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb); + draw_colored_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG, rgb); } - else draw_circle(circlex[k], circley[k], circlerad[k], NSEG); + else draw_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG); init_billiard_color(); } @@ -1157,7 +1232,7 @@ void draw_billiard() /* draws the billiard boundary */ 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_circle(x_shooter, y_shooter, circles[ncircles-1].radius, NSEG); } init_billiard_color(); @@ -1172,16 +1247,16 @@ void draw_billiard() /* draws the billiard boundary */ } case (D_CIRCLES_IN_GENUSN): { -// for (k=0; k= 1) +// for (k=0; k= 1) for (k=0; k= 1)&&(in_polygon(circlex[k], circley[k], 1.0, NPOLY, APOLY))) + if ((circles[k].active >= 1)&&(in_polygon(circles[k].xc, circles[k].yc, 1.0, NPOLY, APOLY))) { - if (circleactive[k] == 2) + if (circles[k].active == 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); + draw_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG); init_billiard_color(); } @@ -1191,7 +1266,7 @@ void draw_billiard() /* draws the billiard boundary */ 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_circle(x_shooter, y_shooter, circles[ncircles-1].radius, NSEG); } /* draw polygon */ @@ -1211,30 +1286,30 @@ void draw_billiard() /* draws the billiard boundary */ case (D_CIRCLES_IN_TORUS): { rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; - for (k=0; k= 1) + if (circles[k].new >= 1) { - rgb_color_scheme_lum(circlecolor[k], 0.85, rgb1); - newcircle[k]--; + rgb_color_scheme_lum(circles[k].color, 0.85, rgb1); + circles[k].new--; } - else rgb_color_scheme(circlecolor[k], rgb1); - draw_colored_circle(circlex[k], circley[k], circlerad[k], NSEG, rgb1); + else rgb_color_scheme(circles[k].color, rgb1); + draw_colored_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG, rgb1); } } init_billiard_color(); - for (k=0; k= 1) + for (k=0; k= 1) { - if (circleactive[k] == 2) + if (circles[k].active == 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); + draw_circle(circles[k].xc, circles[k].yc, circles[k].radius, NSEG); init_billiard_color(); } @@ -1244,7 +1319,7 @@ void draw_billiard() /* draws the billiard boundary */ 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_circle(x_shooter, y_shooter, circles[ncircles-1].radius, NSEG); } init_billiard_color(); @@ -1257,6 +1332,30 @@ void draw_billiard() /* draws the billiard boundary */ glEnd(); break; } + case (D_POLYLINE): + { + for (k=0; k margin) /* there is an intersection with circle i */ @@ -3977,7 +4076,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ 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]); + phiint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); nt++; } @@ -4066,8 +4165,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ angle = conf[0] - (double)ncirc*DPI; - pos[0] = circlex[ncirc] + circlerad[ncirc]*cos(angle); - pos[1] = circley[ncirc] + circlerad[ncirc]*sin(angle); + pos[0] = circles[ncirc].xc + circles[ncirc].radius*cos(angle); + pos[1] = circles[ncirc].yc + circles[ncirc].radius*sin(angle); *alpha = angle + PID + conf[1]; @@ -4099,10 +4198,10 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ c0 = cos(alpha); s0 = sin(alpha); - for (i=0; i margin) /* there is an intersection with circle i */ @@ -4115,7 +4214,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ 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]); + phiint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); /* test wether intersection is in rectangle */ if ((vabs(xint[nt]) < LAMBDA)&&(vabs(yint[nt]) < 1.0)) nt++; @@ -4209,8 +4308,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ angle = conf[0] - (double)ncirc*DPI; - pos[0] = circlex[ncirc] + circlerad[ncirc]*cos(angle); - pos[1] = circley[ncirc] + circlerad[ncirc]*sin(angle); + pos[0] = circles[ncirc].xc + circles[ncirc].radius*cos(angle); + pos[1] = circles[ncirc].yc + circles[ncirc].radius*sin(angle); *alpha = angle + PID + conf[1]; @@ -4246,10 +4345,10 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ length = 2.0*sin(0.5*omega); cw = cos(omega*0.5); - for (i=0; i margin) /* there is an intersection with circle i */ @@ -4262,7 +4361,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ 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]); + phiint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); /* test wether intersection is in polygon */ if (in_polygon(xint[nt], yint[nt], 1.0, NPOLY, APOLY)) nt++; @@ -4350,8 +4449,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ angle = conf[0] - (double)ncirc*DPI; - pos[0] = circlex[ncirc] + circlerad[ncirc]*cos(angle); - pos[1] = circley[ncirc] + circlerad[ncirc]*sin(angle); + pos[0] = circles[ncirc].xc + circles[ncirc].radius*cos(angle); + pos[1] = circles[ncirc].yc + circles[ncirc].radius*sin(angle); *alpha = angle + PID + conf[1]; @@ -4383,10 +4482,10 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ c0 = cos(alpha); s0 = sin(alpha); - for (i=0; i margin) /* there is an intersection with circle i */ @@ -4399,7 +4498,7 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ 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]); + phiint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); /* test wether intersection is in rectangle */ if ((vabs(xint[nt]) < LAMBDA)&&(vabs(yint[nt]) < 1.0)) nt++; @@ -4477,6 +4576,180 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ } +/****************************************************************************************/ +/* billiard in poylgonal line */ +/****************************************************************************************/ + + int pos_polyline(double conf[2], double pos[2], double *alpha) + /* determine position on boundary of domain */ + /* position varies between 0 and nsides */ + /* returns number of hit setment */ + { + double len, rlarge = 1000.0; + int nside; + + nside = (int)conf[0]; + if (nside >= nsides) nside = nsides - 1; + + if (nside >= 0) /* position is on a side of the polygonal line */ + { + len = conf[0] - (double)nside; + + pos[0] = polyline[nside].x1 + len*(polyline[nside].x2 - polyline[nside].x1); + pos[1] = polyline[nside].y1 + len*(polyline[nside].y2 - polyline[nside].y1); + + *alpha = polyline[nside].angle + conf[1]; + return(nside); + } + else /* position is on an absorbing circle */ + { + pos[0] = rlarge; + pos[1] = 0.0; + + *alpha = 0.0; + return(-1); + } + } + + + int vpolyline_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double c0, s0, a, b, c, t, dx, delta, s, xi, yi, margin = 1.0e-12, tmin, rlarge = 1000.0; + double tval[nsides + ncircles], xint[nsides + ncircles], yint[nsides + ncircles], sint[nsides + ncircles]; + int i, nt = 0, nsegment[nsides + ncircles], ntmin; + + c0 = cos(alpha); + s0 = sin(alpha); + + for (i=0; i margin) /* there is an intersection with the line containing segment i */ + { + t = -(a*pos[0] + b*pos[1] + c)/delta; + if (t > margin) + { + xi = pos[0] + t*c0; + yi = pos[1] + t*s0; + dx = polyline[i].x2 - polyline[i].x1; + + if (vabs(dx) > margin) s = (xi - polyline[i].x1)/dx; + else s = (yi - polyline[i].y1)/(polyline[i].y2 - polyline[i].y1); +// printf("s = %.2f\n", s); + + if ((s >= 0.0)&&(s <= 1.0)) + /* the intersection is on the segment */ + { + nsegment[nt] = i; + tval[nt] = t; + sint[nt] = s; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; +// printf("s = %.2f, x = %.2f, y = %.2f\n", s, xint[nt], yint[nt]); + nt++; + } + } + } + } + + if (ABSORBING_CIRCLES) for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nsegment[nt] = -1-i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + sint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); + + nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= 0) + { + config[0] = (double)nsegment[ntmin] + sint[ntmin]; + config[1] = polyline[nsegment[ntmin]].angle - alpha; + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + } + /* set dummy coordinates if circles are absorbing */ + else if ((ABSORBING_CIRCLES)&&(nsegment[ntmin] < 0)) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + 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]; + + +// print_config(config); + + return(nsegment[ntmin]); + } + else /* there is no intersection - set dummy values */ + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + config[2] = 0.0; + config[3] = rlarge; + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = rlarge*cos(alpha); + config[7] = rlarge*sin(alpha); + + return(-1); + } + } + + int vpolyline(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_polyline(config, pos, &alpha); + + c = vpolyline_xy(config, alpha, pos); + + return(c); + } + /****************************************************************************************/ /* general billiard */ /****************************************************************************************/ @@ -4580,6 +4853,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(pos_circles_in_torus(conf, pos, alpha)); break; } + case (D_POLYLINE): + { + return(pos_polyline(conf, pos, alpha)); + break; + } default: { printf("Function pos_billiard not defined for this billiard \n"); @@ -4688,6 +4966,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(vcircles_in_torus_xy(config, alpha, pos)); break; } + case (D_POLYLINE): + { + return(vpolyline_xy(config, alpha, pos)); + break; + } default: { printf("Function vbilliard_xy not defined for this billiard \n"); @@ -4837,6 +5120,13 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(vcircles_in_torus(config)); break; } + case (D_POLYLINE): + { + c = pos_polyline(config, pos, &alpha); + + return(vpolyline(config)); + break; + } default: { printf("Function vbilliard not defined for this billiard \n"); @@ -5009,7 +5299,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ { condition = 1; for (k=0; k YMAX) y -= height; - circlerad[n] = MU; - circleactive[n] = 1; + circles[n].radius = MU; + circles[n].active = 1; } break; @@ -5161,8 +5460,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ case (C_GOLDEN_SPIRAL): { ncircles = 1; - circlex[0] = 0.0; - circley[0] = 0.0; + circles[0].xc = 0.0; + circles[0].yc = 0.0; gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ phi = 0.0; @@ -5179,17 +5478,17 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU)) { - circlex[ncircles] = x; - circley[ncircles] = y; + circles[ncircles].xc = x; + circles[ncircles].yc = y; ncircles++; } } for (i=0; i YMIN - MU)) circleactive[i] = 1; + if ((circles[i].yc < YMAX + MU)&&(circles[i].yc > YMIN - MU)) circles[i].active = 1; } break; } @@ -5201,10 +5500,10 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ for (j = 0; j < NCY; j++) { n = NCY*i + j; - circlex[n] = ((double)(i-NCX/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; + circles[n].xc = ((double)(i-NCX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].yc = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].radius = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); + circles[n].active = 1; } break; } @@ -5213,10 +5512,10 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ ncircles = NPOISSON; for (n = 0; n < NPOISSON; n++) { - circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); - circley[n] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; - circlerad[n] = MU; - circleactive[n] = 1; + circles[n].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[n].yc = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; + circles[n].radius = MU; + circles[n].active = 1; } break; } @@ -5224,8 +5523,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ { printf("Generating Poisson disc sample\n"); /* generate first circle */ - circlex[0] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); - circley[0] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; + circles[0].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[0].yc = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; active_poisson[0] = 1; n_p_active = 1; ncircles = 1; @@ -5242,24 +5541,24 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ { r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0); phi = DPI*(double)rand()/RAND_MAX; - x = circlex[i] + r*cos(phi); - y = circley[i] + r*sin(phi); + x = circles[i].xc + r*cos(phi); + y = circles[i].yc + r*sin(phi); // printf("Testing new circle at (%.3f,%.3f)\t", x, y); far = 1; for (k=0; k= dpoisson*dpoisson); + far = far*((x - circles[k].xc)*(x - circles[k].xc) + (y - circles[k].yc)*(y - circles[k].yc) >= dpoisson*dpoisson); /* new circle is in domain */ far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > YMIN); } if (far) /* accept new circle */ { printf("New circle at (%.3f,%.3f) accepted\n", x, y); - circlex[ncircles] = x; - circley[ncircles] = y; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = x; + circles[ncircles].xc = y; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; active_poisson[ncircles] = 1; ncircles++; n_p_active++; @@ -5296,35 +5595,25 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ yy[2] = -yy[0]; yy[3] = -yy[1]; -// 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]; + circles[4*i + j].xc = xx[i]; + circles[4*i + j].yc = yy[j]; } - circlex[ncircles - 1] = x_target; - circley[ncircles - 1] = y_target; + circles[ncircles - 1].xc = x_target; + circles[ncircles - 1].yc = y_target; for (i=0; i NMAXPOLY) + { + printf("NMAXPOLY has to be increased to at least %i\n", nsides); + nsides = NMAXPOLY; + } + + for (i=0; i DPI) angle -= DPI; + while (angle < 0.0) angle += DPI; + } + + for (i=0; i YMAX) y -= height; - circlerad[n] = MU; - circleactive[n] = 1; - circlecolor[n] = 0; - newcircle[n] = 0; + circles[n].radius = MU; + circles[n].active = 1; + circles[n].color = 0; + circles[n].new = 0; } break; } case (C_GOLDEN_SPIRAL): { ncircles = 1; - circlex[0] = 0.0; - circley[0] = 0.0; + circles[0].xc = 0.0; + circles[0].yc = 0.0; gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ phi = 0.0; @@ -154,19 +156,19 @@ void init_circle_config_pinball(int circle_pattern) if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU)) { - circlex[ncircles] = x; - circley[ncircles] = y; + circles[ncircles].xc = x; + circles[ncircles].yc = y; ncircles++; } } for (i=0; i YMIN - MU)) circleactive[i] = 1; + if ((circles[i].yc < YMAX + MU)&&(circles[i].yc > YMIN - MU)) circles[i].active = 1; } break; } @@ -178,12 +180,12 @@ void init_circle_config_pinball(int circle_pattern) for (j = 0; j < NCY; j++) { n = NCY*i + j; - circlex[n] = ((double)(i-NCX/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; - circlecolor[n] = 0; - newcircle[n] = 0; + circles[n].xc = ((double)(i-NCX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].yc = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].radius = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); + circles[n].active = 1; + circles[n].color = 0; + circles[n].new = 0; } break; } @@ -192,52 +194,52 @@ void init_circle_config_pinball(int circle_pattern) ncircles = NPOISSON; for (n = 0; n < NPOISSON; n++) { - circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); - circley[n] = (BOXYMAX - BOXYMIN)*(double)rand()/RAND_MAX + BOXYMIN; - circlerad[n] = MU; - circleactive[n] = 1; - circlecolor[n] = 0; - newcircle[n] = 0; - double_circle[n] = 0; - partner_circle[n] = n; + circles[n].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[n].yc = (BOXYMAX - BOXYMIN)*(double)rand()/RAND_MAX + BOXYMIN; + circles[n].radius = MU; + circles[n].active = 1; + circles[n].color = 0; + circles[n].new = 0; + circles[n].double_circle = 0; + circles[n].partner = n; /* take care of periodic boundary conditions */ if (B_DOMAIN == D_CIRCLES_IN_TORUS) { /* inactivate circles in corners */ - if ((vabs(circlex[n]) > LAMBDA - MU)&&((circley[n] < BOXYMIN + MU)||(circley[n] > BOXYMAX - MU))) - circleactive[n] = 0; + if ((vabs(circles[n].xc) > LAMBDA - MU)&&((circles[n].yc < BOXYMIN + MU)||(circles[n].yc > BOXYMAX - MU))) + circles[n].active = 0; - if (circlex[n] < - LAMBDA + MU) + if (circles[n].xc < - LAMBDA + MU) { - circlex[ncircles] = circlex[n] + 2.0*LAMBDA; - circley[ncircles] = circley[n]; - partner_circle[ncircles] = n; - partner_circle[n] = ncircles; + circles[ncircles].xc = circles[n].xc + 2.0*LAMBDA; + circles[ncircles].yc = circles[n].yc; + circles[ncircles].partner = n; + circles[n].partner = ncircles; ncircles++; } - else if (circlex[n] > LAMBDA - MU) + else if (circles[n].xc > LAMBDA - MU) { - circlex[ncircles] = circlex[n] - 2.0*LAMBDA; - circley[ncircles] = circley[n]; - partner_circle[ncircles] = n; - partner_circle[n] = ncircles; + circles[ncircles].xc = circles[n].xc - 2.0*LAMBDA; + circles[ncircles].yc = circles[n].yc; + circles[ncircles].partner = n; + circles[n].partner = ncircles; ncircles++; } - if (circley[n] < BOXYMIN + MU) + if (circles[n].yc < BOXYMIN + MU) { - circlex[ncircles] = circlex[n]; - circley[ncircles] = circley[n] + BOXYMAX - BOXYMIN; - partner_circle[ncircles] = n; - partner_circle[n] = ncircles; + circles[ncircles].xc = circles[n].xc; + circles[ncircles].yc = circles[n].yc + BOXYMAX - BOXYMIN; + circles[ncircles].partner = n; + circles[n].partner = ncircles; ncircles++; } - else if (circley[n] > BOXYMAX - MU) + else if (circles[n].yc > BOXYMAX - MU) { - circlex[ncircles] = circlex[n]; - circley[ncircles] = circley[n] - BOXYMAX + BOXYMIN; - partner_circle[ncircles] = n; - partner_circle[n] = ncircles; + circles[ncircles].xc = circles[n].xc; + circles[ncircles].yc = circles[n].yc - BOXYMAX + BOXYMIN; + circles[ncircles].partner = n; + circles[n].partner = ncircles; ncircles++; } } @@ -245,12 +247,12 @@ void init_circle_config_pinball(int circle_pattern) printf("%i circles\n", ncircles); if (B_DOMAIN == D_CIRCLES_IN_TORUS) for (n = NPOISSON; n < ncircles; n++) { -// printf("circle %i at (%.3f, %.3f)\n", n, circlex[n], circley[n]); - circlerad[n] = MU; - if (circleactive[partner_circle[n]]) circleactive[n] = 1; - circlecolor[n] = 0; - newcircle[n] = 0; - double_circle[n] = 1; +// printf("circle %i at (%.3f, %.3f)\n", n, circles[n].xc, circles[n].yc); + circles[n].radius = MU; + if (circles[circles[n].partner].active) circles[n].active = 1; + circles[n].color = 0; + circles[n].new = 0; + circles[n].double_circle = 1; } break; @@ -259,9 +261,11 @@ void init_circle_config_pinball(int circle_pattern) { printf("Generating Poisson disc sample\n"); /* generate first circle */ - circlex[0] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); - circley[0] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; + circles[0].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[0].yc = (BOXYMAX - BOXYMIN)*(double)rand()/RAND_MAX + BOXYMIN; active_poisson[0] = 1; + circles[0].double_circle = 0; + circles[0].partner = 0; n_p_active = 1; ncircles = 1; @@ -277,30 +281,119 @@ void init_circle_config_pinball(int circle_pattern) { r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0); phi = DPI*(double)rand()/RAND_MAX; - x = circlex[i] + r*cos(phi); - y = circley[i] + r*sin(phi); + x = circles[i].xc + r*cos(phi); + y = circles[i].yc + r*sin(phi); // printf("Testing new circle at (%.3f,%.3f)\t", x, y); far = 1; for (k=0; k= dpoisson*dpoisson); - /* new circle is in domain */ - far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > YMIN); - } + far = far*((x - xk)*(x - xk) + (y - yk)*(y - yk) >= dpoisson*dpoisson); + far = far*((x - xk)*(x - xk) + (y - yk + hy)*(y - yk + hy) >= dpoisson*dpoisson); + far = far*((x - xk)*(x - xk) + (y - yk - hy)*(y - yk - hy) >= dpoisson*dpoisson); + far = far*((x - xk + 2.0*LAMBDA)*(x - xk + 2.0*LAMBDA) + (y - yk)*(y - yk) >= dpoisson*dpoisson); + far = far*((x - xk - 2.0*LAMBDA)*(x - xk - 2.0*LAMBDA) + (y - yk)*(y - yk) >= dpoisson*dpoisson); + } + /* new circle is in domain */ + far = far*(vabs(x) < LAMBDA + 0.0)*(y < BOXYMAX + 0.0)*(y > BOXYMIN - 0.0); + + /* exclude circles in corners */ + if ((x > LAMBDA - MU)&&((y < BOXYMIN + MU)||(y > BOXYMAX - MU))) + far = 0; + if (far) /* accept new circle */ { printf("New circle at (%.3f,%.3f) accepted\n", x, y); - circlex[ncircles] = x; - circley[ncircles] = y; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; - circlecolor[ncircles] = 0; - newcircle[ncircles] = 0; + circles[ncircles].xc = x; + circles[ncircles].yc = y; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; + circles[ncircles].color = 0; + circles[ncircles].new = 0; active_poisson[ncircles] = 1; + circles[ncircles].double_circle = 0; + circles[ncircles].partner = ncircles; ncircles++; n_p_active++; naccepted++; + + /* take care of periodic boundary conditions */ + if (B_DOMAIN == D_CIRCLES_IN_TORUS) + { + n = ncircles - 1; + if (x < - LAMBDA + MU) + { + circles[ncircles].xc = x + 2.0*LAMBDA; + circles[ncircles].yc = y; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; + circles[ncircles].color = 0; + circles[ncircles].new = 0; + active_poisson[ncircles] = 1; + circles[n].double_circle = 1; + circles[ncircles].double_circle = 0; + circles[n].partner = ncircles; + circles[ncircles].partner = n; + ncircles++; + n_p_active++; + naccepted++; + } + else if (x > LAMBDA - MU) + { + circles[ncircles].xc = x - 2.0*LAMBDA; + circles[ncircles].yc = y; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; + circles[ncircles].color = 0; + circles[ncircles].new = 0; + active_poisson[ncircles] = 1; + circles[n].double_circle = 1; + circles[ncircles].double_circle = 0; + circles[n].partner = ncircles; + circles[ncircles].partner = n; + ncircles++; + n_p_active++; + naccepted++; + } + if (y < BOXYMIN + MU) + { + circles[ncircles].xc = x; + circles[ncircles].yc = y + BOXYMAX - BOXYMIN; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; + circles[ncircles].color = 0; + circles[ncircles].new = 0; + active_poisson[ncircles] = 1; + circles[n].double_circle = 1; + circles[ncircles].double_circle = 0; + circles[n].partner = ncircles; + circles[ncircles].partner = n; + ncircles++; + n_p_active++; + naccepted++; + } + else if (y > BOXYMAX - MU) + { + circles[ncircles].xc = x; + circles[ncircles].yc = y - BOXYMAX + BOXYMIN; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; + circles[ncircles].color = 0; + circles[ncircles].new = 0; + active_poisson[ncircles] = 1; + circles[n].double_circle = 1; + circles[ncircles].double_circle = 0; + circles[n].partner = ncircles; + circles[ncircles].partner = n; + ncircles++; + n_p_active++; + naccepted++; + } + } } // else printf("Rejected\n"); } @@ -314,6 +407,9 @@ void init_circle_config_pinball(int circle_pattern) } printf("Generated %i circles\n", ncircles); + +// for (i=0; i= 1.0) return(0); + + omega = DPI/((double)polygon.nsides); + cw = cos(omega*0.5); + for (k=0; k YMIN - MU)) circleactive[n] = 1; - else circleactive[n] = 0; + if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; + else circles[n].active = 0; } break; } @@ -410,10 +465,10 @@ void init_circle_config() for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; - 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; + circles[n].xc = ((double)(i-NGRIDX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].yc = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].radius = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); + circles[n].active = 1; } break; } @@ -425,12 +480,12 @@ void init_circle_config() for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; - circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy; - circley[n] = YMIN + ((double)j + 0.5)*dy; - circlerad[n] = MU; + circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; + circles[n].yc = YMIN + ((double)j + 0.5)*dy; + circles[n].radius = MU; p = (double)rand()/RAND_MAX; - if (p < P_PERCOL) circleactive[n] = 1; - else circleactive[n] = 0; + if (p < P_PERCOL) circles[n].active = 1; + else circles[n].active = 0; } break; } @@ -439,10 +494,10 @@ void init_circle_config() ncircles = NPOISSON; for (n = 0; n < NPOISSON; n++) { - circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); - circley[n] = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; - circlerad[n] = MU; - circleactive[n] = 1; + circles[n].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[n].yc = (YMAX - YMIN)*(double)rand()/RAND_MAX + YMIN; + circles[n].radius = MU; + circles[n].active = 1; } break; } @@ -455,10 +510,10 @@ void init_circle_config() n = 5*i + j; phi = (double)i*DPI/40.0; r = LAMBDA*0.5*(1.0 + (double)j/5.0); - circlex[n] = r*cos(phi); - circley[n] = r*sin(phi); - circlerad[n] = MU; - circleactive[n] = 1; + circles[n].xc = r*cos(phi); + circles[n].yc = r*sin(phi); + circles[n].radius = MU; + circles[n].active = 1; } break; } @@ -476,10 +531,10 @@ void init_circle_config() n = 5*i + j; phi = (double)i*DPI/40.0; r = LAMBDA*sa[j]; - circlex[n] = r*cos(phi); - circley[n] = r*sin(phi); - circlerad[n] = LAMBDA*ra[j]; - circleactive[n] = 1; + circles[n].xc = r*cos(phi); + circles[n].yc = r*sin(phi); + circles[n].radius = LAMBDA*ra[j]; + circles[n].active = 1; } break; } @@ -500,22 +555,22 @@ void init_circle_config() for (i = 0; i < 4; i++) for (j = 0; j < 4; j++) { - circlex[4*i + j] = xx[i]; - circley[4*i + j] = yy[j]; + circles[4*i + j].xc = xx[i]; + circles[4*i + j].yc = yy[j]; } - circlex[ncircles - 1] = X_TARGET; - circley[ncircles - 1] = Y_TARGET; + circles[ncircles - 1].xc = X_TARGET; + circles[ncircles - 1].yc = Y_TARGET; for (i=0; i= dpoisson*dpoisson); + far = far*((x - circles[k].xc)*(x - circles[k].xc) + (y - circles[k].yc)*(y - circles[k].yc) >= dpoisson*dpoisson); /* new circle is in domain */ far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > YMIN); } if (far) /* accept new circle */ { printf("New circle at (%.3f,%.3f) accepted\n", x, y); - circlex[ncircles] = x; - circley[ncircles] = y; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = x; + circles[ncircles].yc = y; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; active_poisson[ncircles] = 1; ncircles++; n_p_active++; @@ -587,32 +642,32 @@ void init_circle_config() dx = 2.0*LAMBDA/((double)ncircles); for (n = 0; n < ncircles; n++) { - circlex[n] = -LAMBDA + n*dx; - circley[n] = y; + circles[n].xc = -LAMBDA + n*dx; + circles[n].yc = y; y += height*gamma; if (y > YMAX) y -= height; - circlerad[n] = MU; - circleactive[n] = 1; + circles[n].radius = MU; + circles[n].active = 1; } /* test for circles that overlap top or bottom boundary */ ncirc0 = ncircles; for (n=0; n < ncirc0; n++) { - if (circley[n] + circlerad[n] > YMAX) + if (circles[n].yc + circles[n].radius > YMAX) { - circlex[ncircles] = circlex[n]; - circley[ncircles] = circley[n] - height; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = circles[n].xc; + circles[ncircles].yc = circles[n].yc - height; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; ncircles ++; } - else if (circley[n] - circlerad[n] < YMIN) + else if (circles[n].yc - circles[n].radius < YMIN) { - circlex[ncircles] = circlex[n]; - circley[ncircles] = circley[n] + height; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = circles[n].xc; + circles[ncircles].yc = circles[n].yc + height; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; ncircles ++; } } @@ -621,8 +676,8 @@ void init_circle_config() case (C_GOLDEN_SPIRAL): { ncircles = 1; - circlex[0] = 0.0; - circley[0] = 0.0; + circles[0].xc = 0.0; + circles[0].yc = 0.0; gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ phi = 0.0; @@ -639,18 +694,18 @@ void init_circle_config() if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU)) { - circlex[ncircles] = x; - circley[ncircles] = y; + circles[ncircles].xc = x; + circles[ncircles].yc = y; ncircles++; } } for (i=0; i YMIN - MU)) circleactive[i] = 1; -// printf("i = %i, circlex = %.3lg, circley = %.3lg\n", i, circlex[i], circley[i]); + if ((circles[i].yc < YMAX + MU)&&(circles[i].yc > YMIN - MU)) circles[i].active = 1; +// printf("i = %i, circlex = %.3lg, circley = %.3lg\n", i, circles[i].xc, circles[i].yc); } break; } @@ -663,37 +718,37 @@ void init_circle_config() 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; + circles[n].xc = ((double)(i-NGRIDX/2) + 0.5)*dy; /* is +0.5 needed? */ + circles[n].yc = YMIN + ((double)j - 0.5)*dy; + if (((i+NGRIDX)%4 == 2)||((i+NGRIDX)%4 == 3)) circles[n].yc += 0.5*dy; + circles[n].radius = MU; /* activate only circles that intersect the domain */ - if ((circley[n] < YMAX + MU)&&(circley[n] > YMIN - MU)) circleactive[n] = 1; - else circleactive[n] = 0; + if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; + else circles[n].active = 0; } break; } case (C_ONE): { - circlex[ncircles] = 0.0; - circley[ncircles] = 0.0; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = 0.0; + circles[ncircles].yc = 0.0; + circles[ncircles].radius = MU; + circles[ncircles].active = 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; + circles[ncircles].xc = 0.0; + circles[ncircles].yc = 0.0; + circles[ncircles].radius = MU; + circles[ncircles].active = 2; ncircles += 1; - circlex[ncircles] = 0.0; - circley[ncircles] = 0.0; - circlerad[ncircles] = 2.0*MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = 0.0; + circles[ncircles].yc = 0.0; + circles[ncircles].radius = 2.0*MU; + circles[ncircles].active = 1; ncircles += 1; break; } @@ -709,7 +764,30 @@ void init_circle_config() } } +void init_polygon_config(t_polygon polygons[NMAXCIRCLES]) +/* initialise the polygon configuration, for billiard shape D_CIRCLES */ +/* uses init_circle_config, this is where C++ would be more elegant */ +{ + int i; + t_circle circle[NMAXCIRCLES]; + + init_circle_config(circle); + for (i=0; ix = z.x - 2.0*zperp.x; + zprime->y = z.y - 2.0*zperp.y; + + return(1); +} -void compute_isospectral_coordinates(int type, double xshift, double yshift, double scaling, double vertex[9][2]) +int compute_tokarsky_coordinates(double xshift, double yshift, double scaling, + t_vertex polyline[NMAXPOLY]) +/* compute positions of vertices of tokarsky room */ +{ + int i; + double pos[2]; + + polyline[0].x = 0.0; polyline[0].y = 2.0; + polyline[1].x = 1.0; polyline[1].y = 3.0; + polyline[2].x = 1.0; polyline[2].y = 4.0; + polyline[3].x = 2.0; polyline[3].y = 4.0; + polyline[4].x = 2.0; polyline[4].y = 3.0; + polyline[5].x = 3.0; polyline[5].y = 3.0; + polyline[6].x = 3.0; polyline[6].y = 2.0; + polyline[7].x = 5.0; polyline[7].y = 2.0; + polyline[8].x = 5.0; polyline[8].y = 3.0; + polyline[9].x = 6.0; polyline[9].y = 3.0; + + polyline[10].x = 6.0; polyline[10].y = 4.0; + polyline[11].x = 7.0; polyline[11].y = 3.0; + polyline[12].x = 8.0; polyline[12].y = 3.0; + polyline[13].x = 8.0; polyline[13].y = 2.0; + polyline[14].x = 7.0; polyline[14].y = 2.0; + polyline[15].x = 7.0; polyline[15].y = 1.0; + polyline[16].x = 6.0; polyline[16].y = 0.0; + polyline[17].x = 6.0; polyline[17].y = 1.0; + polyline[18].x = 5.0; polyline[18].y = 1.0; + polyline[19].x = 4.0; polyline[19].y = 0.0; + + polyline[20].x = 4.0; polyline[20].y = 1.0; + polyline[21].x = 3.0; polyline[21].y = 1.0; + polyline[22].x = 2.0; polyline[22].y = 0.0; + polyline[23].x = 2.0; polyline[23].y = 1.0; + polyline[24].x = 1.0; polyline[24].y = 1.0; + polyline[25].x = 1.0; polyline[25].y = 2.0; + + for (i=0; i<26; i++) + { + polyline[i].x = (polyline[i].x + xshift)*scaling; + polyline[i].y = (polyline[i].y + yshift)*scaling; + xy_to_pos(polyline[i].x, polyline[i].y, pos); + polyline[i].posi = pos[0]; + polyline[i].posj = pos[1]; + } + return(26); +} + +void compute_isospectral_coordinates(int type, int ishift, double xshift, double yshift, double scaling, + t_vertex polyline[NMAXPOLY]) /* compute positions of vertices of isospectral billiards */ /* central triangle has coordinates (0,0), (1,0) and (LAMBDA,MU) fed into affine transformation */ /* defined by (xshift - 0.5), (yshift - 0.25) and scaling*/ { - vertex[0][0] = (xshift - 0.5)*scaling; - vertex[0][1] = (yshift - 0.25)*scaling; + int i; + double pos[2]; - vertex[1][0] = (0.5 + xshift)*scaling; - vertex[1][1] = (yshift - 0.25)*scaling; + polyline[ishift].x = (xshift - 0.5)*scaling; + polyline[ishift].y = (yshift - 0.25)*scaling; - vertex[2][0] = (LAMBDA + xshift - 0.5)*scaling; - vertex[2][1] = (MU + yshift - 0.25)*scaling; + polyline[ishift+1].x = (0.5+xshift)*scaling; + polyline[ishift+1].y = (yshift - 0.25)*scaling; - axial_symmetry(vertex[0], vertex[2], vertex[1], vertex[3]); - axial_symmetry(vertex[0], vertex[1], vertex[2], vertex[4]); - axial_symmetry(vertex[1], vertex[2], vertex[0], vertex[5]); + polyline[ishift+2].x = (LAMBDA+xshift - 0.5)*scaling; + polyline[ishift+2].y = (MU+yshift - 0.25)*scaling; + + axial_symmetry_tvertex(polyline[ishift], polyline[ishift+2], polyline[ishift+1], &polyline[ishift+3]); + axial_symmetry_tvertex(polyline[ishift], polyline[ishift+1], polyline[ishift+2], &polyline[ishift+4]); + axial_symmetry_tvertex(polyline[ishift+1], polyline[ishift+2], polyline[ishift], &polyline[ishift+5]); if (type == 0) { - axial_symmetry(vertex[0], vertex[3], vertex[2], vertex[6]); - axial_symmetry(vertex[1], vertex[4], vertex[0], vertex[7]); - axial_symmetry(vertex[2], vertex[5], vertex[1], vertex[8]); + axial_symmetry_tvertex(polyline[ishift], polyline[ishift+3], polyline[ishift+2], &polyline[ishift+6]); + axial_symmetry_tvertex(polyline[ishift+1], polyline[ishift+4], polyline[ishift], &polyline[ishift+7]); + axial_symmetry_tvertex(polyline[ishift+2], polyline[ishift+5], polyline[ishift+1], &polyline[ishift+8]); } else { - axial_symmetry(vertex[2], vertex[3], vertex[0], vertex[6]); - axial_symmetry(vertex[0], vertex[4], vertex[1], vertex[7]); - axial_symmetry(vertex[1], vertex[5], vertex[2], vertex[8]); + axial_symmetry_tvertex(polyline[ishift+2], polyline[ishift+3], polyline[ishift], &polyline[ishift+6]); + axial_symmetry_tvertex(polyline[ishift], polyline[ishift+4], polyline[ishift+1], &polyline[ishift+7]); + axial_symmetry_tvertex(polyline[ishift+1], polyline[ishift+5], polyline[ishift+2], &polyline[ishift+8]); + } + + for (i=ishift; i NMAXPOLY) + { + printf("NMAXPOLY needs to be increased to %i\n", nsides); + nsides = NMAXPOLY; + } + + for (i=0; i DPI) angle -= DPI; + while (angle < 0.0) angle += DPI; + + xy_to_pos(x*MU, y*MU, pos); + polyline[i].posi = pos[0]; + polyline[i].posj = pos[1]; + } + + return(nsides); +} + + +int init_polyline(int depth, t_vertex polyline[NMAXPOLY]) +/* initialise variable polyline, for certain polygonal domain shapes */ +{ + switch (B_DOMAIN) { + case (D_TOKARSKY): + { + return(compute_tokarsky_coordinates(-4.0, -2.0, (XMAX - XMIN)/8.4, polyline)); + } + case (D_ISOSPECTRAL): + { + compute_isospectral_coordinates(0, 0, ISO_XSHIFT_LEFT, ISO_YSHIFT_LEFT, ISO_SCALE, polyline); + compute_isospectral_coordinates(1, 9, ISO_XSHIFT_RIGHT, ISO_YSHIFT_RIGHT, ISO_SCALE, polyline); + return(18); + } + case (D_HOMOPHONIC): + { + compute_homophonic_coordinates(0, 0, ISO_XSHIFT_LEFT, ISO_YSHIFT_LEFT, ISO_SCALE, polyline); + compute_homophonic_coordinates(1, 22, ISO_XSHIFT_RIGHT, ISO_YSHIFT_RIGHT, ISO_SCALE, polyline); + return(44); + } + case (D_VONKOCH): + { + return(compute_vonkoch_coordinates(depth, polyline)); + } + case (D_VONKOCH_HEATED): + { + return(compute_vonkoch_coordinates(depth, polyline)); + } + default: + { + return(0); + } } } @@ -785,6 +1123,15 @@ void isospectral_initial_point(double x, double y, double left[2], double right[ right[1] = (y + ISO_YSHIFT_RIGHT)*ISO_SCALE; } +void homophonic_initial_point(double xleft, double yleft, double xright, double yright, double left[2], double right[2]) +/* compute initial coordinates in isospectral billiards */ +{ + left[0] = (xleft + ISO_XSHIFT_LEFT)*ISO_SCALE; + left[1] = (yleft + ISO_YSHIFT_LEFT)*ISO_SCALE; + right[0] = (xright + ISO_XSHIFT_RIGHT)*ISO_SCALE; + right[1] = (yright + ISO_YSHIFT_RIGHT)*ISO_SCALE; +} + int xy_in_triangle(double x, double y, double z1[2], double z2[2], double z3[2]) /* returns 1 iff (x,y) is inside the triangle with vertices z1, z2, z3 */ { @@ -799,15 +1146,28 @@ int xy_in_triangle(double x, double y, double z1[2], double z2[2], double z3[2]) else return(0); } +int xy_in_triangle_tvertex(double x, double y, t_vertex z1, t_vertex z2, t_vertex z3) +/* returns 1 iff (x,y) is inside the triangle with vertices z1, z2, z3 */ +{ + double v1, v2, v3; + + /* compute wedge products */ + v1 = (z2.x - z1.x)*(y - z1.y) - (z2.y - z1.y)*(x - z1.x); + v2 = (z3.x - z2.x)*(y - z2.y) - (z3.y - z2.y)*(x - z2.x); + v3 = (z1.x - z3.x)*(y - z3.y) - (z1.y - z3.y)*(x - z3.x); + + if ((v1 >= 0.0)&&(v2 >= 0.0)&&(v3 >= 0.0)) return(1); + else return(0); +} + int xy_in_billiard(double x, double y) /* returns 1 if (x,y) represents a point in the billiard */ // double x, y; { double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width; - int i, j, k, k1, k2, condition; - static double vertex[9][2], wertex[9][2]; - static int first = 1; + int i, j, k, k1, k2, condition, m; + static int first = 1, nsides; switch (B_DOMAIN) { case (D_RECTANGLE): @@ -1041,39 +1401,90 @@ int xy_in_billiard(double x, double y) } case (D_ISOSPECTRAL): { - if (first) - { - compute_isospectral_coordinates(0, ISO_XSHIFT_LEFT, ISO_YSHIFT_LEFT, ISO_SCALE, vertex); - compute_isospectral_coordinates(1, ISO_XSHIFT_RIGHT, ISO_YSHIFT_RIGHT, ISO_SCALE, wertex); - first = 0; - } /* 1st triangle */ - condition = xy_in_triangle(x, y, vertex[0], vertex[1], vertex[2]); - condition += xy_in_triangle(x, y, vertex[0], vertex[4], vertex[1]); - condition += xy_in_triangle(x, y, vertex[1], vertex[5], vertex[2]); - condition += xy_in_triangle(x, y, vertex[0], vertex[2], vertex[3]); - condition += xy_in_triangle(x, y, vertex[1], vertex[4], vertex[7]); - condition += xy_in_triangle(x, y, vertex[2], vertex[5], vertex[8]); - condition += xy_in_triangle(x, y, vertex[0], vertex[3], vertex[6]); + condition = xy_in_triangle_tvertex(x, y, polyline[0], polyline[1], polyline[2]); + condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[4], polyline[1]); + condition += xy_in_triangle_tvertex(x, y, polyline[1], polyline[5], polyline[2]); + condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[2], polyline[3]); + condition += xy_in_triangle_tvertex(x, y, polyline[1], polyline[4], polyline[7]); + condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[5], polyline[8]); + condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[3], polyline[6]); /* 2nd triangle */ - condition += xy_in_triangle(x, y, wertex[0], wertex[1], wertex[2]); - condition += xy_in_triangle(x, y, wertex[0], wertex[4], wertex[1]); - condition += xy_in_triangle(x, y, wertex[1], wertex[5], wertex[2]); - condition += xy_in_triangle(x, y, wertex[0], wertex[2], wertex[3]); - condition += xy_in_triangle(x, y, wertex[0], wertex[7], wertex[4]); - condition += xy_in_triangle(x, y, wertex[1], wertex[8], wertex[5]); - condition += xy_in_triangle(x, y, wertex[2], wertex[6], wertex[3]); + condition += xy_in_triangle_tvertex(x, y, polyline[9], polyline[10], polyline[11]); + condition += xy_in_triangle_tvertex(x, y, polyline[9], polyline[13], polyline[10]); + condition += xy_in_triangle_tvertex(x, y, polyline[10], polyline[14], polyline[11]); + condition += xy_in_triangle_tvertex(x, y, polyline[9], polyline[11], polyline[12]); + condition += xy_in_triangle_tvertex(x, y, polyline[9], polyline[16], polyline[13]); + condition += xy_in_triangle_tvertex(x, y, polyline[10], polyline[17], polyline[14]); + condition += xy_in_triangle_tvertex(x, y, polyline[11], polyline[15], polyline[12]); + return(condition >= 1); + } + case (D_HOMOPHONIC): + { + /* conditions could be summarised in larger triangles, but this is to keep */ + /* the option of using triangles with other angles than 30-60-90 */ + + /* 1st triangle */ + condition = xy_in_triangle_tvertex(x, y, polyline[2], polyline[0], polyline[1]); + condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[1], polyline[3]); + condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[21], polyline[1]); + condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[10], polyline[21]); + condition += xy_in_triangle_tvertex(x, y, polyline[10], polyline[11], polyline[21]); + condition += xy_in_triangle_tvertex(x, y, polyline[11], polyline[13], polyline[21]); + condition += xy_in_triangle_tvertex(x, y, polyline[11], polyline[12], polyline[13]); + condition += xy_in_triangle_tvertex(x, y, polyline[13], polyline[14], polyline[21]); + condition += xy_in_triangle_tvertex(x, y, polyline[14], polyline[20], polyline[21]); + condition += xy_in_triangle_tvertex(x, y, polyline[14], polyline[15], polyline[20]); + condition += xy_in_triangle_tvertex(x, y, polyline[15], polyline[19], polyline[20]); + + condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[4], polyline[5]); + condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[5], polyline[7]); + condition += xy_in_triangle_tvertex(x, y, polyline[5], polyline[6], polyline[7]); + condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[7], polyline[8]); + condition += xy_in_triangle_tvertex(x, y, polyline[2], polyline[8], polyline[0]); + condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[8], polyline[9]); + condition += xy_in_triangle_tvertex(x, y, polyline[0], polyline[9], polyline[10]); + + condition += xy_in_triangle_tvertex(x, y, polyline[15], polyline[16], polyline[19]); + condition += xy_in_triangle_tvertex(x, y, polyline[16], polyline[17], polyline[18]); + condition += xy_in_triangle_tvertex(x, y, polyline[16], polyline[18], polyline[19]); + + /* 2nd triangle */ + condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+0], polyline[22+1]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+1], polyline[22+3]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+0], polyline[22+21], polyline[22+1]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+0], polyline[22+10], polyline[22+21]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+10], polyline[22+11], polyline[22+21]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+13], polyline[22+21]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+12], polyline[22+13]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+13], polyline[22+14], polyline[22+21]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+14], polyline[22+20], polyline[22+21]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+14], polyline[22+15], polyline[22+20]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+15], polyline[22+19], polyline[22+20]); + + condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+3], polyline[22+5]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+3], polyline[22+4], polyline[22+5]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+5], polyline[22+6]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+6], polyline[22+8]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+2], polyline[22+8], polyline[22+9]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+6], polyline[22+7], polyline[22+8]); + + condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+10], polyline[22+16]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+16], polyline[22+18]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+11], polyline[22+18], polyline[22+12]); + condition += xy_in_triangle_tvertex(x, y, polyline[22+16], polyline[22+17], polyline[22+18]); + return(condition >= 1); } case (D_CIRCLES): { for (i = 0; i < ncircles; i++) - if (circleactive[i]) + if (circles[i].active) { - x1 = circlex[i]; - y1 = circley[i]; - r2 = circlerad[i]*circlerad[i]; + x1 = circles[i].xc; + y1 = circles[i].yc; + r2 = circles[i].radius*circles[i].radius; if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(0); } return(1); @@ -1081,16 +1492,36 @@ int xy_in_billiard(double x, double y) case (D_CIRCLES_IN_RECT): /* returns 2 inside circles, 0 outside rectangle */ { for (i = 0; i < ncircles; i++) - if (circleactive[i]) + if (circles[i].active) { - x1 = circlex[i]; - y1 = circley[i]; - r2 = circlerad[i]*circlerad[i]; + x1 = circles[i].xc; + y1 = circles[i].yc; + r2 = circles[i].radius*circles[i].radius; if ((x-x1)*(x-x1) + (y-y1)*(y-y1) < r2) return(2); } if ((vabs(x) >= LAMBDA)||(vabs(y) >= 1.0)) return(0); else return(1); } + case (D_POLYGONS): + { + for (i = 0; i < ncircles; i++) + if ((polygons[i].active)&&(in_tpolygon(x, y, polygons[i]))) return(0); + return(1); + } + case (D_VONKOCH): + { + condition = xy_in_triangle_tvertex(x, y, polyline[0], polyline[npolyline/3], polyline[2*npolyline/3]); + m = 1; + k = 1; + for (i = 0; i < MDEPTH; i++) + { + m = m*4; + for (j = 0; j < npolyline/m; j++) + condition += xy_in_triangle_tvertex(x, y, polyline[j*m + k], polyline[j*m + 2*k], polyline[j*m + 3*k]); + k = k*4; + } + return(condition >= 1); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -1226,6 +1657,25 @@ int xy_in_billiard(double x, double y) // else if ((vabs(x) > XMAX - 0.01)||(vabs(y) > YMAX - 0.01)) return(2); else return(1); } + case (D_VONKOCH_HEATED): + { + if (x*x + y*y > LAMBDA*LAMBDA) return(2); + + x1 = x; + y1 = y; + condition = xy_in_triangle_tvertex(x1, y1, polyline[0], polyline[npolyline/3], polyline[2*npolyline/3]); + m = 1; + k = 1; + for (i = 0; i < MDEPTH; i++) + { + m = m*4; + for (j = 0; j < npolyline/m; j++) + condition += xy_in_triangle_tvertex(x1, y1, polyline[j*m + k], polyline[j*m + 2*k], polyline[j*m + 3*k]); + k = k*4; + } + if (condition > 0) return(0); + else return(1); + } default: { printf("Function ij_in_billiard not defined for this billiard \n"); @@ -1244,34 +1694,18 @@ int ij_in_billiard(int i, int j) return(xy_in_billiard(xy[0], xy[1])); } -void toka_lineto(double x1, double y1) -/* draws boundary segments of Tokarsky billiard */ -{ - double ratio, x, y, pos[2]; - - ratio = (XMAX - XMIN)/8.4; - x = ratio*(x1 - 4.0); - y = ratio*(y1 - 2.0); - xy_to_pos(x, y, pos); - glVertex2d(pos[0], pos[1]); -} - -void iso_lineto(double z[2]) +void tvertex_lineto(t_vertex z) /* draws boundary segments of isospectral billiard */ { - double pos[2]; - - xy_to_pos(z[0], z[1], pos); - glVertex2d(pos[0], pos[1]); + glVertex2d(z.posi, z.posj); } void draw_billiard() /* draws the billiard boundary */ { double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax; - static double vertex[9][2], wertex[9][2]; int i, j, k, k1, k2, mr2; - static int first = 1; + static int first = 1, nsides; if (BLACK) glColor3f(1.0, 1.0, 1.0); else glColor3f(0.0, 0.0, 0.0); @@ -1802,32 +2236,7 @@ void draw_billiard() /* draws the billiard boundary */ case (D_TOKARSKY): { glBegin(GL_LINE_LOOP); - toka_lineto(0.0, 2.0); - toka_lineto(1.0, 3.0); - toka_lineto(1.0, 4.0); - toka_lineto(2.0, 4.0); - toka_lineto(2.0, 3.0); - toka_lineto(3.0, 3.0); - toka_lineto(3.0, 2.0); - toka_lineto(5.0, 2.0); - toka_lineto(5.0, 3.0); - toka_lineto(6.0, 3.0); - toka_lineto(6.0, 4.0); - toka_lineto(7.0, 3.0); - toka_lineto(8.0, 3.0); - toka_lineto(8.0, 2.0); - toka_lineto(7.0, 2.0); - toka_lineto(7.0, 1.0); - toka_lineto(6.0, 0.0); - toka_lineto(6.0, 1.0); - toka_lineto(5.0, 1.0); - toka_lineto(4.0, 0.0); - toka_lineto(4.0, 1.0); - toka_lineto(3.0, 1.0); - toka_lineto(2.0, 0.0); - toka_lineto(2.0, 1.0); - toka_lineto(1.0, 1.0); - toka_lineto(1.0, 2.0); + for (i=0; i= COL_TURBO) color_scheme_asym(COLOR_SCHEME, value, 1.0, 1, rgb); + else color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); + break; + } + case (P_MEAN_ENERGY): + { + value = dy_e*(double)(j - jmin)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, value, 1.0, 1, rgb); + else color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); + break; + } + case (P_PHASE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); + break; + } + } + glColor3f(rgb[0], rgb[1], rgb[2]); + if (ROTATE_COLOR_SCHEME) + { + glVertex2i(j, imin); + glVertex2i(j, imax); + glVertex2i(j+1, imax); + glVertex2i(j+1, imin); + } + else + { + glVertex2i(imin, j); + glVertex2i(imax, j); + glVertex2i(imax, j+1); + glVertex2i(imin, j+1); + } + } + glEnd (); + + glColor3f(1.0, 1.0, 1.0); + glLineWidth(BOUNDARY_WIDTH); + draw_rectangle(x1, y1, x2, y2); +} + + diff --git a/sub_wave_comp.c b/sub_wave_comp.c index 4c4309e..e994bcc 100644 --- a/sub_wave_comp.c +++ b/sub_wave_comp.c @@ -3,7 +3,7 @@ short int circletop[NMAXCIRCLES]; /* set to 1 if circle is in top half */ -void init_circle_config_half(int pattern, int top) +void init_circle_config_half(int pattern, int top, t_circle circles[NMAXCIRCLES]) /* initialise the arrays circlex, circley, circlerad and circleactive */ /* for billiard shape D_CIRCLES */ { @@ -21,12 +21,12 @@ void init_circle_config_half(int pattern, int top) { 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; + circles[n].xc = ((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; + if (top) circles[n].yc = ymean + y; + else circles[n].yc = ymean - y; + circles[n].radius = MU; + circles[n].active = 1; circletop[n] = top; } ncircles += NGRIDX*NGRIDY/2; @@ -40,13 +40,13 @@ void init_circle_config_half(int pattern, int top) for (j = 0; j < NGRIDY/2+1; j++) { n = ncircles + (NGRIDY+1)*i/2 + j; - circlex[n] = ((double)(i-NGRIDX/2))*dy; + circles[n].xc = ((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; + if (top) circles[n].yc = ymean + 0.5*dy + y; + else circles[n].yc = ymean - 0.5*dy - y; + circles[n].radius = MU; + circles[n].active = 1; circletop[n] = top; } ncircles += NGRIDX*(NGRIDY+1)/2; @@ -59,15 +59,15 @@ void init_circle_config_half(int pattern, int top) 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; + circles[n].xc = ((double)(i-NGRIDX/2) + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + if (NGRIDX%2 == 0) circles[n].xc += 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; + if (top) circles[n].yc = ymean + y; + else circles[n].yc = ymean - y; + circles[n].radius = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); + circles[n].active = 1; circletop[n] = top; - printf("n = %i, x = %.3lg\n", n, circlex[n]); + printf("n = %i, x = %.3lg\n", n, circles[n].xc); } ncircles += NGRIDX*NGRIDY/2; break; @@ -79,14 +79,14 @@ void init_circle_config_half(int pattern, int top) for (j = 0; j < NGRIDY/2; j++) { n = ncircles + NGRIDY*i/2 + j; - circlex[n] = ((double)(i-NGRIDX/2) + 0.5)*dy; + circles[n].xc = ((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; + circles[n].yc = y; + circles[n].radius = MU; p = (double)rand()/RAND_MAX; - if (p < P_PERCOL) circleactive[n] = 1; - else circleactive[n] = 0; + if (p < P_PERCOL) circles[n].active = 1; + else circles[n].active = 0; circletop[n] = top; } ncircles += NGRIDX*NGRIDY/2; @@ -97,15 +97,15 @@ void init_circle_config_half(int pattern, int top) for (i = 0; i < NPOISSON/2; i++) { n = ncircles + i; - circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[n].xc = 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; + circles[n].yc = y; + circles[n].radius = MU; + circles[n].active = 1; circletop[n] = top; - printf("n = %i, x = %.3lg\n", n, circlex[n]); + printf("n = %i, x = %.3lg\n", n, circles[n].xc); } ncircles += NPOISSON/2; break; @@ -118,12 +118,12 @@ void init_circle_config_half(int pattern, int top) 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); + circles[n].xc = 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; + circles[n].yc = y; + circles[n].radius = MU; + circles[n].active = 1; circletop[n] = top; } ncircles += 200; @@ -142,25 +142,25 @@ void init_circle_config_half(int pattern, int top) 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); + circles[n].xc = r*cos(phi); + circles[n].yc = r*sin(phi); if (top) y = r*sin(phi); else y = -r*sin(phi); - circley[n] = y; + circles[n].yc = y; - circlerad[n] = LAMBDA*ra[j]; - circleactive[n] = 1; + circles[n].radius = LAMBDA*ra[j]; + circles[n].active = 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; + circles[ncircles].xc = 0.0; + circles[ncircles].yc = 0.0; + circles[ncircles].radius = MU; + circles[ncircles].active = 2; circletop[ncircles] = top; ncircles += 1; break; @@ -170,12 +170,12 @@ void init_circle_config_half(int pattern, int top) printf("Generating Poisson disc sample\n"); /* generate first circle */ n = ncircles; - circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + circles[n].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); if (top) y = ymean + (YMAX-ymean)*(double)rand()/RAND_MAX; else y = ymean + (YMIN-ymean)*(double)rand()/RAND_MAX; - circley[n] = y; - circlerad[n] = MU; - circleactive[n] = 1; + circles[n].yc = y; + circles[n].radius = MU; + circles[n].active = 1; circletop[n] = top; active_poisson[n] = 1; n_p_active = 1; @@ -187,21 +187,21 @@ void init_circle_config_half(int pattern, int top) i = rand()%(ncirc0); n = ncircles + i; while (!active_poisson[ncircles + i]) i = rand()%(ncirc0); -// printf("Starting from circle %i at (%.3f,%.3f)\n", i, circlex[i], circley[i]); +// printf("Starting from circle %i at (%.3f,%.3f)\n", i, circles[i].xc, circles[i].yc); /* generate new candidates */ naccepted = 0; for (j=0; j= dpoisson*dpoisson); + far = far*((x - circles[k].xc)*(x - circles[k].xc) + (y - circles[k].yc)*(y - circles[k].yc) >= dpoisson*dpoisson); /* new circle is in domain */ if (top) far = far*(vabs(x) < LAMBDA)*(y < YMAX)*(y > 0.0); else far = far*(vabs(x) < LAMBDA)*(y > YMIN)*(y < 0.0); @@ -210,12 +210,12 @@ void init_circle_config_half(int pattern, int top) { printf("New circle at (%.3f,%.3f) accepted\n", x, y); nnew = ncircles + ncirc0; - circlex[nnew] = x; - circley[nnew] = y; - circlerad[nnew] = MU; - circleactive[nnew] = 1; + circles[nnew].xc = x; + circles[nnew].yc = y; + circles[nnew].radius = MU; + circles[nnew].active = 1; active_poisson[nnew] = 1; - circleactive[nnew] = 1; +// circleactive[nnew] = 1; circletop[nnew] = top; ncirc0++; n_p_active++; @@ -248,8 +248,8 @@ void init_circle_config_half(int pattern, int top) else y = ymean - 0.5*height; for (n = 0; n < 150; n++) { - circlex[ncircles + n] = -LAMBDA + n*dx; - circley[ncircles + n] = y; + circles[ncircles + n].xc = -LAMBDA + n*dx; + circles[ncircles + n].yc = y; if (top) { @@ -263,8 +263,8 @@ void init_circle_config_half(int pattern, int top) } - circlerad[ncircles + n] = MU; - circleactive[ncircles + n] = 1; + circles[ncircles + n].radius = MU; + circles[ncircles + n].active = 1; circletop[ncircles] = top; } @@ -279,20 +279,20 @@ void init_circle_config_half(int pattern, int top) for (n=0; n < 150; n++) { - if (circley[ncirc0 + n] + circlerad[ncirc0 + n] > ytop) + if (circles[ncirc0+n].yc + circles[ncirc0 + n].radius > ytop) { - circlex[ncircles] = circlex[ncirc0 + n]; - circley[ncircles] = circley[ncirc0 + n] - height; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = circles[ncirc0 + n].xc; + circles[ncircles].yc = circles[ncirc0+n].yc - height; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; ncircles ++; } - else if (circley[ncirc0 + n] - circlerad[ncirc0 + n] < ybottom) + else if (circles[ncirc0+n].yc - circles[ncirc0 + n].radius < ybottom) { - circlex[ncircles] = circlex[ncirc0 + n]; - circley[ncircles] = circley[ncirc0 + n] + height; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = circles[ncirc0 + n].xc; + circles[ncircles].yc = circles[ncirc0+n].yc + height; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; ncircles ++; } } @@ -301,10 +301,10 @@ void init_circle_config_half(int pattern, int top) } case (C_GOLDEN_SPIRAL): { - circlex[ncircles] = 0.0; - circley[ncircles] = 0.0; - circlerad[ncircles] = MU; - circleactive[ncircles] = 1; + circles[ncircles].xc = 0.0; + circles[ncircles].yc = 0.0; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; circletop[ncircles] = top; gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ @@ -322,13 +322,13 @@ void init_circle_config_half(int pattern, int top) if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU)) { - circlex[ncircles] = x; - circley[ncircles] = y; - circlerad[ncircles] = MU; - if (((top)&&(circley[ncircles] < YMAX + MU)&&(circley[ncircles] > ymean - MU)) - ||((!top)&&(circley[ncircles] < ymean + MU)&&(circley[ncircles] > YMIN - MU))) + circles[ncircles].xc = x; + circles[ncircles].yc = y; + circles[ncircles].radius = MU; + if (((top)&&(circles[ncircles].yc < YMAX + MU)&&(circles[ncircles].yc > ymean - MU)) + ||((!top)&&(circles[ncircles].yc < ymean + MU)&&(circles[ncircles].yc > YMIN - MU))) { - circleactive[ncircles] = 1; + circles[ncircles].active = 1; circletop[ncircles] = top; ncircles++; } @@ -338,30 +338,30 @@ void init_circle_config_half(int pattern, int top) } case (C_ONE): { - circlex[ncircles] = 0.0; - circley[ncircles] = 0.0; - if (top) circlerad[ncircles] = MU; - else circlerad[ncircles] = MUB; - circleactive[ncircles] = 1; + circles[ncircles].xc = 0.0; + circles[ncircles].yc = 0.0; + if (top) circles[ncircles].radius = MU; + else circles[ncircles].radius = MUB; + circles[ncircles].active = 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; + circles[ncircles].xc = 0.0; + circles[ncircles].yc = 0.0; + if (top) circles[ncircles].radius = MU; + else circles[ncircles].radius = MUB; + circles[ncircles].active = 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; + circles[ncircles].xc = 0.0; + circles[ncircles].yc = 0.0; + if (top) circles[ncircles].radius = 2.0*MU; + else circles[ncircles].radius = 2.0*MUB; + circles[ncircles].active = 1; circletop[ncircles] = top; ncircles += 1; break; @@ -379,21 +379,21 @@ void init_circle_config_half(int pattern, int top) } -void init_circle_config_comp() +void init_circle_config_comp(t_circle circles[NMAXCIRCLES]) /* 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); + init_circle_config_half(CIRCLE_PATTERN, 1, circles); + init_circle_config_half(CIRCLE_PATTERN_B, 0, circles); } -void init_circle_config_energy() +void init_circle_config_energy(t_circle circles[NMAXCIRCLES]) /* initialise the arrays circlex, circley, circlerad and circleactive */ /* for billiard shape D_CIRCLES */ { ncircles = 0; - init_circle_config_half(CIRCLE_PATTERN, 0); + init_circle_config_half(CIRCLE_PATTERN, 0, circles); } int xy_in_billiard_half(double x, double y, int domain, int pattern, int top) @@ -419,15 +419,15 @@ int xy_in_billiard_half(double x, double y, int domain, int pattern, int top) case D_CIRCLES: { for (i = 0; i < ncircles; i++) - if (circleactive[i] != 0) + if (circles[i].active != 0) { - /* choose specific type according to value of circleactive[i] */ - if (circleactive[i] == 1) type = 0; - else type = circleactive[i]; + /* choose specific type according to value of circles[i].active */ + if (circles[i].active == 1) type = 0; + else type = circles[i].active; - x1 = circlex[i]; - y1 = circley[i]; - r2 = circlerad[i]*circlerad[i]; + x1 = circles[i].xc; + y1 = circles[i].yc; + r2 = circles[i].radius*circles[i].radius; 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); } @@ -535,14 +535,14 @@ void draw_billiard_half(int domain, int pattern, int top) /* draws the bill { glLineWidth(2); for (i = 0; i < ncircles; i++) - if ((circleactive[i])&&(circletop[i] == top)) + if ((circles[i].active)&&(circletop[i] == top)) { glBegin(GL_LINE_STRIP); 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); + x = circles[i].xc + circles[i].radius*cos(phi); + y = circles[i].yc + circles[i].radius*sin(phi); if ((top)&&(circletop[i])&&(y < 0.0)) y = 0.0; if ((!top)&&(!circletop[i])&&(y > 0.0)) y = 0.0; xy_to_pos(x, y, pos); @@ -681,9 +681,9 @@ void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX], { /* find leftmost and rightmost circle */ for (i=0; i xright)) xright = circlex[i] + circlerad[i]; + if ((circles[i].active)&&(circles[i].xc + circles[i].radius > xright)) xright = circles[i].xc + circles[i].radius; xy_to_ij(xleft, 0.0, ij); ileft = ij[0]; diff --git a/wave_billiard.c b/wave_billiard.c index d6cbe6f..ad685c8 100644 --- a/wave_billiard.c +++ b/wave_billiard.c @@ -62,35 +62,38 @@ /* Choice of the billiard table */ -#define B_DOMAIN 32 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 38 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 7 /* pattern of circles, see list in global_pdes.c */ +#define CIRCLE_PATTERN 2 /* pattern of circles or polygons, 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 RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.0 /* parameter controlling the dimensions of domain */ -#define MU 1.0 /* parameter controlling the dimensions of domain */ -#define NPOLY 7 /* 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 LAMBDA 0.9 /* parameter controlling the dimensions of domain */ +#define MU 0.03 /* parameter controlling the dimensions of domain */ +#define NPOLY 6 /* number of sides of polygon */ +#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define MDEPTH 5 /* 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 16 /* number of grid point for grid of disks */ -#define NGRIDY 20 /* number of grid point for grid of disks */ +#define NGRIDX 10 /* number of grid point for grid of disks */ +#define NGRIDY 12 /* number of grid point for grid of disks */ +// #define NGRIDX 16 /* number of grid point for grid of disks */ +// #define NGRIDY 20 /* number of grid point for grid of disks */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 #define X_TARGET 0.4 #define Y_TARGET 0.7 /* shooter and target positions in laser fight */ -#define ISO_XSHIFT_LEFT -1.65 -#define ISO_XSHIFT_RIGHT 0.4 -#define ISO_YSHIFT_LEFT -0.05 -#define ISO_YSHIFT_RIGHT -0.05 -#define ISO_SCALE 0.85 /* coordinates for isospectral billiards */ +#define ISO_XSHIFT_LEFT -2.9 +#define ISO_XSHIFT_RIGHT 1.4 +#define ISO_YSHIFT_LEFT -0.15 +#define ISO_YSHIFT_RIGHT -0.15 +#define ISO_SCALE 0.5 /* coordinates for isospectral billiards */ /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard below */ @@ -104,9 +107,11 @@ #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.02 /* Courant number in medium B */ +#define COURANTB 0.01 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMAB 5.0e-3 /* damping factor in wave equation */ +// #define GAMMAB 5.0e-3 /* damping factor in wave equation */ +// #define GAMMAB 1.0e-2 /* 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-7 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ @@ -119,12 +124,13 @@ /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 2 +#define B_COND 3 /* Parameters for length and speed of simulation */ -#define NSTEPS 5000 /* number of frames of movie */ -#define NVID 50 /* number of iterations between images displayed on screen */ +// #define NSTEPS 500 /* number of frames of movie */ +#define NSTEPS 1500 /* number of frames of movie */ +#define NVID 40 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ #define INITIAL_TIME 0 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ @@ -134,31 +140,35 @@ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ #define MID_FRAMES 20 /* number of still frames between parts of two-part movie */ -#define END_FRAMES 100 /* number of still frames at end of movie */ +#define END_FRAMES 50 /* number of still frames at end of movie */ /* Parameters of initial condition */ #define INITIAL_AMP 0.75 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */ +// #define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ +// #define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ +#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ #define PLOT 1 -#define PLOT_B 3 /* plot type for second movie */ +#define PLOT_B 0 /* plot type for second movie */ /* Color schemes */ -#define COLOR_PALETTE 13 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ #define COLOR_SCHEME 3 /* 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 0.15 /* sensitivity of color on wave amplitude */ +// #define SLOPE 0.25 /* sensitivity of color on wave amplitude */ +#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ +// #define E_SCALE 150.0 /* scaling factor for energy representation */ #define E_SCALE 100.0 /* scaling factor for energy representation */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -166,9 +176,14 @@ #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 180.0 /* mean value of hue for color scheme C_HUE */ +// #define HUEMEAN 210.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ +// #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 1 /* set to 1 to draw color scheme horizontally */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ @@ -188,8 +203,7 @@ 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], +void evolve_wave_half_old(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 */ @@ -306,7 +320,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX /* add oscillating boundary condition on the left */ if ((i == 0)&&(OSCILLATE_LEFT)) phi_out[i][j] = AMPLITUDE*cos((double)time*OMEGA); - psi_out[i][j] = x; +// psi_out[i][j] = x; if (FLOOR) { @@ -322,12 +336,271 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[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]) +/* time step of field evolution */ +/* phi is value of field at time t, psi at time t-1 */ +/* this version of the function has been rewritten in order to minimize the number of if-branches */ +{ + int i, j, iplus, iminus, jplus, jminus; + double delta, x, y, c, cc, gamma; + static long time = 0; + static double tc[NX][NY], tcc[NX][NY], tgamma[NX][NY]; + static short int first = 1; + + time++; + + /* initialize tables with wave speeds and dissipation */ + if (first) + { + for (i=0; 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; + } + } + } +} + + 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); +// evolve_wave_half_old(phi, psi, phi_tmp, psi_tmp, xy_in); +// evolve_wave_half_old(phi_tmp, psi_tmp, phi, psi, xy_in); +} + + +void draw_color_bar(int plot, double range) +{ + if (ROTATE_COLOR_SCHEME) draw_color_scheme(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range); + else draw_color_scheme(1.7, YMIN + 0.1, 1.9, YMAX - 0.1, plot, -range, range); } @@ -358,7 +631,13 @@ void animation() } /* initialise positions and radii of circles */ - if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) init_circle_config(); + if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) init_circle_config(circles); + else if (B_DOMAIN == D_POLYGONS) init_polygon_config(polygons); + printf("Polygons initialized\n"); + + /* initialise polyline for von Koch and simular domains */ + npolyline = init_polyline(MDEPTH, polyline); + for (i=0; i= COL_TURBO) color_scheme_asym(COLOR_SCHEME, value, 1.0, 1, rgb); - else color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); - break; - } - case (P_MEAN_ENERGY): - { - value = dy_e*(double)(j - jmin); - if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, value, 1.0, 1, rgb); - else color_scheme(COLOR_SCHEME, value, 1.0, 1, rgb); - break; - } - } - glColor3f(rgb[0], rgb[1], rgb[2]); - glVertex2i(imin, j); - glVertex2i(imax, j); - glVertex2i(imax, j+1); - glVertex2i(imin, j+1); - } - glEnd (); - - glColor3f(1.0, 1.0, 1.0); - draw_rectangle(x1, y1, x2, y2); -} - diff --git a/wave_comparison.c b/wave_comparison.c index 6539b91..26e8eb9 100644 --- a/wave_comparison.c +++ b/wave_comparison.c @@ -67,6 +67,7 @@ #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 RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ #define MU 0.03 /* parameter controlling the dimensions of domain */ @@ -174,6 +175,12 @@ #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 DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 4.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ + + /* For debugging purposes only */ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 5.0 /* max value of wave amplitude */ @@ -191,7 +198,7 @@ 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], +void evolve_wave_half_old(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 */ @@ -338,6 +345,368 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX // printf("phi(0,0) = %.3lg, psi(0,0) = %.3lg\n", phi[NX/2][NY/2], psi[NX/2][NY/2]); } +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; + static double tc[NX][NY], tcc[NX][NY], tgamma[NX][NY]; + static short int first = 1; + + time++; + + /* initialize tables with wave speeds and dissipation */ + if (first) + { + for (i=0; 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 */ @@ -368,7 +737,7 @@ void animation() /* initialise positions and radii of circles */ printf("initializing circle configuration\n"); - if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES)) init_circle_config_comp(); + if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES)) init_circle_config_comp(circles); courant2 = COURANT*COURANT; courantb2 = COURANTB*COURANTB; diff --git a/wave_energy.c b/wave_energy.c index af88bac..86d15eb 100644 --- a/wave_energy.c +++ b/wave_energy.c @@ -66,11 +66,12 @@ #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 2 /* pattern of circles, see list in global_pdes.c */ #define CIRCLE_PATTERN_B 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 RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ #define LAMBDA 0.75 /* parameter controlling the dimensions of domain */ #define MU 0.03 /* parameter controlling the dimensions of domain */ @@ -123,6 +124,7 @@ /* Boundary conditions, see list in global_pdes.c */ +// #define B_COND 2 #define B_COND 3 /* Parameters for length and speed of simulation */ @@ -152,14 +154,14 @@ /* Color schemes */ -#define COLOR_PALETTE 0 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ -#define COLOR_SCHEME 1 /* choice of color scheme, see list in global_pdes.c */ +#define COLOR_SCHEME 3 /* 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 50.0 /* sensitivity of color on wave amplitude */ +#define SLOPE 10.0 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define E_SCALE 500.0 /* scaling factor for energy representation */ @@ -170,6 +172,11 @@ #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 DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 4.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ + /* For debugging purposes only */ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 5.0 /* max value of wave amplitude */ @@ -227,16 +234,29 @@ void draw_wave_energy(double *phi[NX], double *psi[NX], short int *xy_in[NX], do 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); + if (((TWOSPEEDS)&&(xy_in[i][j] != 2))||(xy_in[i][j] == 1)) { + switch (PLOT) { + case (P_AMPLITUDE): + { + /* make wave luminosity larger inside obstacles */ + if (!(xy_in[i][j])) color_scheme_lum(COLOR_SCHEME, phi[i][j], scale, time, 0.7, rgb); + else color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); + break; + } + case (P_ENERGY): + { + energy = compute_energy(phi, psi, xy_in, i, j); + /* adjust energy to color palette */ + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym(COLOR_SCHEME, energy, scale, time, rgb); + else color_scheme(COLOR_SCHEME, energy, scale, time, rgb); + break; + } + case (P_MIXED): + { + if (j > NY/2) color_scheme(COLOR_SCHEME, phi[i][j], scale, time, rgb); + else color_scheme(COLOR_SCHEME, compute_energy(phi, psi, xy_in, i, j), scale, time, rgb); + break; + } } glColor3f(rgb[0], rgb[1], rgb[2]); @@ -360,7 +380,7 @@ void draw_wave_energy(double *phi[NX], double *psi[NX], short int *xy_in[NX], do /* animation part */ /*********************/ -void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], +void evolve_wave_half_old(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 */ @@ -393,28 +413,22 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX { iplus = (i+1); if (iplus == NX) iplus = NX-1; iminus = (i-1); if (iminus == -1) iminus = 0; - jplus = (j+1); - if (jplus == jmid) jplus = jmid-1; - jminus = (j-1); - if (jminus == -1) jminus = 0; + jplus = (j+1); if (jplus == jmid) jplus = jmid-1; + jminus = (j-1); if (jminus == -1) jminus = 0; } else if (B_COND == BC_PERIODIC) { iplus = (i+1) % NX; - iminus = (i-1) % NX; - if (iminus < 0) iminus += NX; + iminus = (i-1) % NX; if (iminus < 0) iminus += NX; jplus = (j+1) % jmid; - jminus = (j-1) % jmid; - if (jminus < 0) jminus += jmid; + jminus = (j-1) % jmid; if (jminus < 0) 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; - jplus = (j+1); - if (jplus >= jmid) jplus -= jmid; - jminus = (j-1); - if (jminus < 0) jminus += jmid; + jplus = (j+1); if (jplus >= jmid) jplus -= jmid; + jminus = (j-1); if (jminus < 0) jminus += jmid; } /* imposing linear wave on top and bottom by making Laplacian 1d */ @@ -486,10 +500,261 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[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]) +/* 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; + static double tc[NX][NY/2], tcc[NX][NY/2], tgamma[NX][NY/2]; + static short int first = 1; + + time++; + + /* initialize tables with wave speeds and dissipation */ + if (first) + { + for (i=0; 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_old(phi, psi, phi_tmp, psi_tmp, xy_in); +// evolve_wave_half_old(phi_tmp, psi_tmp, phi, psi, xy_in); evolve_wave_half(phi, psi, phi_tmp, psi_tmp, xy_in); evolve_wave_half(phi_tmp, psi_tmp, phi, psi, xy_in); } @@ -515,8 +780,9 @@ void animation() /* initialise positions and radii of circles */ printf("initializing circle configuration\n"); - if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES)) init_circle_config_energy(); -// if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES)) init_circle_config_comp(); + if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN_B == D_CIRCLES)) init_circle_config_energy(circles); + else if (B_DOMAIN == D_POLYGONS) init_polygon_config(polygons); + courant2 = COURANT*COURANT; courantb2 = COURANTB*COURANTB;