diff --git a/global_pdes.c b/global_pdes.c index 916a2dc..59bb027 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -29,6 +29,8 @@ #define D_MENGER 15 /* Menger-Sierpinski carpet */ #define D_JULIA_INT 16 /* interior of Julia set */ #define D_MENGER_ROTATED 17 /* rotated Menger-Sierpinski carpet */ +#define D_PARABOLA 18 /* parabolic domain */ +#define D_TWO_PARABOLAS 19 /* two facing parabolic antennas */ #define D_CIRCLES 20 /* several circles */ @@ -42,6 +44,8 @@ #define C_CLOAK 5 /* invisibility cloak */ #define C_CLOAK_A 6 /* first optimized invisibility cloak */ +#define C_POISSON_DISC 8 /* Poisson disc sampling */ + #define C_GOLDEN_MEAN 10 /* pattern based on vertical shifts by golden mean */ #define C_GOLDEN_SPIRAL 11 /* spiral pattern based on golden mean */ #define C_SQUARE_HEX 12 /* alternating between square and hexagonal/triangular */ diff --git a/sub_wave.c b/sub_wave.c index 9bda253..d49e366 100644 --- a/sub_wave.c +++ b/sub_wave.c @@ -122,7 +122,7 @@ void color_scheme(int scheme, double value, double scale, int time, double rgb[3 /* saturation = r, luminosity = y */ switch (scheme) { - case C_LUM: + case (C_LUM): { hue = COLORHUE + (double)time*COLORDRIFT/(double)NSTEPS; if (hue < 0.0) hue += 360.0; @@ -135,7 +135,7 @@ void color_scheme(int scheme, double value, double scale, int time, double rgb[3 hsl_to_rgb(hue, r, y, rgb); break; } - case C_HUE: + case (C_HUE): { r = 0.9; amplitude = color_amplitude(value, scale, time); @@ -315,6 +315,33 @@ void erase_area(double x, double y, double dx, double dy) glEnd(); } + +void erase_area_rgb(double x, double y, double dx, double dy, double rgb[3]) +{ + double pos[2]; + + glColor3f(rgb[0], rgb[1], rgb[2]); + glBegin(GL_QUADS); + xy_to_pos(x - dx, y - dy, pos); + glVertex2d(pos[0], pos[1]); + xy_to_pos(x + dx, y - dy, pos); + glVertex2d(pos[0], pos[1]); + xy_to_pos(x + dx, y + dy, pos); + glVertex2d(pos[0], pos[1]); + xy_to_pos(x - dx, y + dy, pos); + glVertex2d(pos[0], pos[1]); + glEnd(); +} + + +void erase_area_hsl(double x, double y, double dx, double dy, double h, double s, double l) +{ + double pos[2], rgb[3]; + + hsl_to_rgb(h, s, l, rgb); + erase_area_rgb(x, y, dx, dy, rgb); +} + void draw_rectangle(double x1, double y1, double x2, double y2) { double pos[2]; @@ -398,8 +425,9 @@ void init_circle_config() /* initialise the arrays circlex, circley, circlerad and circleactive */ /* for billiard shape D_CIRCLES */ { - int i, j, n, ncirc0; - double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma; + int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted; + double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = 3.25*MU; + short int active_poisson[NMAXCIRCLES], far; switch (CIRCLE_PATTERN) { case (C_SQUARE): @@ -517,6 +545,65 @@ void init_circle_config() } break; } + case (C_POISSON_DISC): + { + 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; + active_poisson[0] = 1; + n_p_active = 1; + ncircles = 1; + + while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES)) + { + /* randomly select an active circle */ + i = rand()%(ncircles); + while (!active_poisson[i]) i = rand()%(ncircles); +// printf("Starting from circle %i at (%.3f,%.3f)\n", i, circlex[i], circley[i]); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= 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; + active_poisson[ncircles] = 1; + ncircles++; + n_p_active++; + naccepted++; + } +// else printf("Rejected\n"); + } + if (naccepted == 0) /* inactivate circle i */ + { +// printf("No candidates work, inactivate circle %i\n", i); + active_poisson[i] = 0; + n_p_active--; + } + printf("%i active circles\n", n_p_active); + } + + printf("Generated %i circles\n", ncircles); + break; + } case (C_GOLDEN_MEAN): { ncircles = 300; @@ -590,7 +677,7 @@ void init_circle_config() if ((circley[i] < YMAX + MU)&&(circley[i] > YMIN - MU)) circleactive[i] = 1; // printf("i = %i, circlex = %.3lg, circley = %.3lg\n", i, circlex[i], circley[i]); } - break; + break; } case (C_SQUARE_HEX): { @@ -648,11 +735,11 @@ void init_circle_config() } -int xy_in_billiard(x, y) +int xy_in_billiard(double x, double y) /* returns 1 if (x,y) represents a point in the billiard */ -double x, y; +// double x, y; { - double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy; + double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width; int i, j, k, k1, k2, condition; switch (B_DOMAIN) { @@ -774,6 +861,21 @@ double x, y; } return(1); } + case (D_PARABOLA): + { + return(x > 0.25*y*y/LAMBDA - LAMBDA); + } + case (D_TWO_PARABOLAS): + { + x1 = 0.25*y*y/MU - MU - LAMBDA; + x2 = -x1; + width = 0.25*MU; + if (width > 0.2) width = 0.2; + if (vabs(y) > 1.5*MU) return(1); + else if ((x < x1 - width)||(x > x2 + width)) return(1); + else if ((x > x1)&&(x < x2)) return(1); + else return(0); + } case (D_CIRCLES): { for (i = 0; i < ncircles; i++) @@ -944,7 +1046,7 @@ int ij_in_billiard(int i, int j) 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; + double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width; int i, j, k, k1, k2, mr2; if (BLACK) glColor3f(1.0, 1.0, 1.0); @@ -1205,6 +1307,75 @@ void draw_billiard() /* draws the billiard boundary */ } break; } + case (D_PARABOLA): + { + dy = (YMAX - YMIN)/(double)NSEG; + glBegin(GL_LINE_STRIP); + + for (i = 0; i < NSEG+1; i++) + { + y = YMIN + dy*(double)i; + x = 0.25*y*y/LAMBDA - LAMBDA; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + glEnd (); + + if (FOCI) + { + glColor3f(0.3, 0.3, 0.3); + draw_circle(0.0, 0.0, r, NSEG); + } + break; + } + case (D_TWO_PARABOLAS): + { + dy = 3.0*MU/(double)NSEG; + width = 0.25*MU; + if (width > 0.2) width = 0.2; + glBegin(GL_LINE_LOOP); + for (i = 0; i < NSEG+1; i++) + { + y = -1.5*MU + dy*(double)i; + x = 0.25*y*y/MU - MU - LAMBDA; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + for (i = 0; i < NSEG+1; i++) + { + y = 1.5*MU - dy*(double)i; + x = 0.25*y*y/MU - (MU + width) - LAMBDA; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + glEnd (); + + glBegin(GL_LINE_LOOP); + for (i = 0; i < NSEG+1; i++) + { + y = -1.5*MU + dy*(double)i; + x = LAMBDA + MU - 0.25*y*y/MU; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + for (i = 0; i < NSEG+1; i++) + { + y = 1.5*MU - dy*(double)i; + x = LAMBDA + (MU + width) - 0.25*y*y/MU; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + glEnd (); + + if (FOCI) + { + glColor3f(0.3, 0.3, 0.3); + draw_circle(-LAMBDA, 0.0, r, NSEG); + draw_circle(LAMBDA, 0.0, r, NSEG); + } + + break; + } case (D_CIRCLES): { glLineWidth(BOUNDARY_WIDTH); diff --git a/sub_wave_comp.c b/sub_wave_comp.c index 4b0cb74..4c4309e 100644 --- a/sub_wave_comp.c +++ b/sub_wave_comp.c @@ -7,9 +7,10 @@ void init_circle_config_half(int pattern, int top) /* initialise the arrays circlex, circley, circlerad and circleactive */ /* for billiard shape D_CIRCLES */ { - int i, j, n, ncirc0; - double dx, dy, p, phi, r, ra[5], sa[5], height, y = 0.0, gamma, ymean, ytop, ybottom; - + int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted, nnew; + double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, ymean, ytop, ybottom, dpoisson = 3.05*MU; + short int active_poisson[NMAXCIRCLES], far; + ymean = 0.5*(YMIN + YMAX); switch (pattern) { case (C_SQUARE): @@ -164,6 +165,80 @@ void init_circle_config_half(int pattern, int top) ncircles += 1; break; } + case (C_POISSON_DISC): + { + printf("Generating Poisson disc sample\n"); + /* generate first circle */ + n = ncircles; + circlex[n] = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); + if (top) y = ymean + (YMAX-ymean)*(double)rand()/RAND_MAX; + else y = ymean + (YMIN-ymean)*(double)rand()/RAND_MAX; + circley[n] = y; + circlerad[n] = MU; + circleactive[n] = 1; + circletop[n] = top; + active_poisson[n] = 1; + n_p_active = 1; + ncirc0 = 1; + + while ((n_p_active > 0)&&(ncircles + ncirc0 < NMAXCIRCLES)) + { + /* randomly select an active circle */ + 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]); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= 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); + } + if (far) /* accept new circle */ + { + 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; + active_poisson[nnew] = 1; + circleactive[nnew] = 1; + circletop[nnew] = top; + ncirc0++; + n_p_active++; + naccepted++; + } +// else printf("Rejected\n"); + } + if (naccepted == 0) /* inactivate circle i */ + { +// printf("No candidates work, inactivate circle %i\n", ncircles + i); + active_poisson[ncircles + i] = 0; + n_p_active--; + } + printf("%i active circles\n", n_p_active); +// sleep(1); + } + + printf("Already existing: %i circles\n", ncircles); + ncircles += ncirc0; + printf("Generated %i circles\n", ncirc0); + printf("Total: %i circles\n", ncircles); + break; + } case (C_GOLDEN_MEAN): { gamma = (sqrt(5.0) - 1.0)*0.5; /* golden mean */ @@ -224,6 +299,43 @@ void init_circle_config_half(int pattern, int top) break; } + case (C_GOLDEN_SPIRAL): + { + circlex[ncircles] = 0.0; + circley[ncircles] = 0.0; + circlerad[ncircles] = MU; + circleactive[ncircles] = 1; + circletop[ncircles] = top; + + gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ + phi = 0.0; + r0 = 2.0*MU; + r = r0 + MU; + + for (i=0; i<1000; i++) + { + x = r*cos(phi); + y = r*sin(phi); + + phi += gamma; + r += MU*r0/r; + + if ((vabs(x) < LAMBDA)&&(vabs(y) < YMAX + MU)) + { + circlex[ncircles] = x; + circley[ncircles] = y; + circlerad[ncircles] = MU; + if (((top)&&(circley[ncircles] < YMAX + MU)&&(circley[ncircles] > ymean - MU)) + ||((!top)&&(circley[ncircles] < ymean + MU)&&(circley[ncircles] > YMIN - MU))) + { + circleactive[ncircles] = 1; + circletop[ncircles] = top; + ncircles++; + } + } + } + break; + } case (C_ONE): { circlex[ncircles] = 0.0; @@ -276,6 +388,13 @@ void init_circle_config_comp() init_circle_config_half(CIRCLE_PATTERN_B, 0); } +void init_circle_config_energy() +/* initialise the arrays circlex, circley, circlerad and circleactive */ +/* for billiard shape D_CIRCLES */ +{ + ncircles = 0; + init_circle_config_half(CIRCLE_PATTERN, 0); +} int xy_in_billiard_half(double x, double y, int domain, int pattern, int top) /* returns 1 if (x,y) represents a point in the billiard */ @@ -551,9 +670,12 @@ void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX], /* compute total energy in top/bottom left/right boxes */ { int i, j, ij[2]; - double energy = 0.0, pos, xleft = XMAX, xright = XMIN; + double energy = 0.0, rescale, pos, xleft = XMAX, xright = XMIN, emax = 1.2; + double energy_ij[NX][NY]; static short int first = 1; static int ileft, iright, jmid = NY/2; + static double sqremax; + if (first) /* compute box limits */ { @@ -572,48 +694,74 @@ void compute_energy_tblr(double *phi[NX], double *psi[NX], short int *xy_in[NX], printf("xleft = %.3lg, xright = %.3lg", xleft, xright); } + for (i=0; i emax)) + { + rescale = sqrt(emax/energy_ij[i][j]); + if (j%100 == 0) printf("Rescaling at (%i,%i) by %.5lg\n", i, j, rescale); + phi[i][j] = phi[i][j]*rescale; + psi[i][j] = psi[i][j]*rescale; + } + else if (energy_ij[i][j] > 0.1*emax) + { + rescale = sqrt(0.1*emax/energy_ij[i][j]); + if (j%10 == 0) printf("Rescaling at (%i,%i) by %.5lg\n", i, j, rescale); + phi[i][j] = phi[i][j]*rescale; + psi[i][j] = psi[i][j]*rescale; + } + } + /* top left box */ for (i=0; i /* Sam Leffler's libtiff library. */ #include -#define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define MOVIE 1 /* set to 1 to generate movie */ +#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ /* General geometrical parameters */ @@ -52,11 +52,6 @@ #define NX 1280 /* number of grid points on x axis */ #define NY 720 /* number of grid points on y axis */ -// #define NX 640 /* number of grid points on x axis */ -// #define NY 360 /* number of grid points on y axis */ - -/* setting NX to WINWIDTH and NY to WINHEIGHT increases resolution */ -/* but will multiply run time by 4 */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ @@ -67,15 +62,15 @@ /* Choice of the billiard table */ -#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 19 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 11 /* pattern of circles, see list in global_pdes.c */ +#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_pdes.c */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define LAMBDA 0.85 /* parameter controlling the dimensions of domain */ -#define MU 0.03 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.0 /* parameter controlling the dimensions of domain */ +#define MU 1.25 /* parameter controlling the dimensions of domain */ #define NPOLY 3 /* number of sides of polygon */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ @@ -83,7 +78,7 @@ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 15 /* number of grid point for grid of disks */ +#define NGRIDX 16 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ /* You can add more billiard tables by adapting the functions */ @@ -102,7 +97,7 @@ #define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMAB 1.0e-6 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ -#define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */ +#define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ #define KAPPA_SIDES 5.0e-4 /* "elasticity" term on absorbing boundary */ #define KAPPA_TOPBOT 0.0 /* "elasticity" term on absorbing boundary */ @@ -113,26 +108,28 @@ /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 3 +#define B_COND 2 /* Parameters for length and speed of simulation */ -#define NSTEPS 3000 /* number of frames of movie */ -#define NVID 20 /* number of iterations between images displayed on screen */ +#define NSTEPS 5000 /* number of frames of movie */ +#define NVID 25 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ -#define INITIAL_TIME 250 /* time after which to start saving frames */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ -#define PAUSE 1000 /* number of frames after which to pause */ +#define PAUSE 1000 /* number of frames after which to pause */ #define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ -#define SLEEP2 1 /* final sleeping time */ +#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 */ /* Plot type, see list in global_pdes.c */ -#define PLOT 0 +#define PLOT 1 -#define PLOT_B 1 /* plot type for second movie */ +#define PLOT_B 0 /* plot type for second movie */ /* Color schemes */ @@ -141,9 +138,9 @@ #define COLOR_SCHEME 1 /* choice of color scheme, see list in global_pdes.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ -#define SLOPE 1.5 /* sensitivity of color on wave amplitude */ +#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 2000.0 /* scaling factor for energy representation */ +#define E_SCALE 500.0 /* scaling factor for energy representation */ // #define E_SCALE 2500.0 /* scaling factor for energy representation */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -184,7 +181,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX // c = COURANT; // cc = courant2; - #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y) + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y,c,cc,gamma) for (i=0; i /* Sam Leffler's libtiff library. */ #include -#define MOVIE 0 /* set to 1 to generate movie */ +#define MOVIE 1 /* set to 1 to generate movie */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ @@ -59,11 +59,11 @@ /* Choice of the billiard table */ -#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ -#define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN_B 20 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_pdes.c */ -#define CIRCLE_PATTERN_B 10 /* pattern of circles, see list in global_pdes.c */ +#define CIRCLE_PATTERN 11 /* pattern of circles, see list in global_pdes.c */ +#define CIRCLE_PATTERN_B 8 /* 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 */ @@ -78,7 +78,7 @@ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 15 /* number of grid point for grid of disks */ +#define NGRIDX 20 /* number of grid point for grid of disks */ #define NGRIDY 20 /* number of grid point for grid of disks */ /* You can add more billiard tables by adapting the functions */ @@ -86,16 +86,22 @@ /* Physical parameters of wave equation */ -#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ -#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ -#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ +#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */ +#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ -#define OMEGA 0.0035 /* frequency of periodic excitation */ -#define AMPLITUDE 1.0 /* amplitude of periodic excitation */ +#define OMEGA 0.0 /* frequency of periodic excitation */ +#define AMPLITUDE 0.025 /* amplitude of periodic excitation */ #define COURANT 0.02 /* Courant number */ -#define COURANTB 0.0075 /* Courant number in medium B */ +#define COURANTB 0.004 /* Courant number in medium B */ +// #define COURANTB 0.005 /* Courant number in medium B */ +// #define COURANTB 0.008 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMAB 1.0e-5 /* damping factor in wave equation */ +// #define GAMMA 1.0e-8 /* damping factor in wave equation */ +#define GAMMAB 1.0e-8 /* damping factor in wave equation */ +// #define GAMMAB 1.0e-6 /* damping factor in wave equation */ +// #define GAMMAB 2.0e-4 /* damping factor in wave equation */ +// #define GAMMAB 2.5e-4 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-6 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ @@ -112,7 +118,7 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 6000 /* number of frames of movie */ +#define NSTEPS 3200 /* number of frames of movie */ #define NVID 25 /* number of iterations between images displayed on screen */ #define NSEG 100 /* number of segments of boundary */ #define INITIAL_TIME 200 /* time after which to start saving frames */ @@ -123,10 +129,11 @@ #define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ +#define END_FRAMES 100 /* number of still frames at end of movie */ /* Plot type, see list in global_pdes.c */ -#define PLOT 1 +#define PLOT 0 /* Color schemes */ @@ -135,7 +142,7 @@ #define COLOR_SCHEME 1 /* choice of color scheme, see list in global_pdes.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ -#define SLOPE 1.5 /* sensitivity of color on wave amplitude */ +#define SLOPE 50.0 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define E_SCALE 2000.0 /* scaling factor for energy representation */ @@ -148,7 +155,7 @@ /* For debugging purposes only */ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ -#define VMAX 10.0 /* max value of wave amplitude */ +#define VMAX 5.0 /* max value of wave amplitude */ #include "global_pdes.c" /* constants and global variables */ @@ -174,7 +181,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX time++; - #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y) + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y,c,cc,gamma) for (i=0; i= jmid) jplus -= jmid; + jminus = (j-1); if (jminus < 0) jminus += jmid; } else /* upper half */ @@ -293,6 +301,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; if (FLOOR) @@ -344,7 +353,8 @@ void animation() courantb2 = COURANTB*COURANTB; /* initialize wave with a drop at one point, zero elsewhere */ - int_planar_wave_comp(XMIN + 0.015, 0.0, phi, psi, xy_in); + init_wave_flat_comp(phi, psi, xy_in); +// int_planar_wave_comp(XMIN + 0.015, 0.0, phi, psi, xy_in); // int_planar_wave_comp(XMIN + 0.5, 0.0, phi, psi, xy_in); printf("initializing wave\n"); // int_planar_wave_comp(XMIN + 0.1, 0.0, phi, psi, xy_in); @@ -434,7 +444,7 @@ void animation() if (MOVIE) { - for (i=0; i<20; i++) save_frame(); + for (i=0; i