diff --git a/drop_billiard.c b/drop_billiard.c index cf0a7cd..12a2364 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -42,7 +42,7 @@ /* Choice of the billiard table, see global_particles.c */ -#define B_DOMAIN 15 /* choice of domain shape */ +#define B_DOMAIN 1 /* choice of domain shape */ #define CIRCLE_PATTERN 0 /* pattern of circles */ #define POLYLINE_PATTERN 1 /* pattern of polyline */ @@ -59,7 +59,7 @@ #define NGOLDENSPIRAL 2000 /* max number of points for C_GOLDEN_SPIRAL arrandement */ #define SDEPTH 1 /* Sierpinski gastket depth */ -#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */ +#define LAMBDA 1.5 /* 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 */ @@ -72,13 +72,14 @@ #define DRAW_BILLIARD 1 /* set to 1 to draw billiard */ #define DRAW_CONSTRUCTION_LINES 1 /* set to 1 to draw additional construction lines for billiard */ #define PERIODIC_BC 0 /* set to 1 to enforce periodic boundary conditions when drawing particles */ +#define PENROSE_RATIO 2.5 /* parameter controlling the shape of small ellipses in Penrose room */ -#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */ +#define RESAMPLE 1 /* set to 1 if particles should be added when dispersion too large */ #define NPART 20000 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ -#define NSTEPS 5000 /* number of frames of movie */ +#define NSTEPS 5250 /* number of frames of movie */ #define TIME 150 /* time between movie frames, for fluidity of real-time simulation */ #define DPHI 0.0001 /* integration step */ #define NVID 75 /* number of iterations between images displayed on screen */ @@ -90,25 +91,25 @@ /* simulation parameters */ -#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 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 LMAX 0.1 /* minimal segment length triggering resampling */ +#define LPERIODIC 0.5 /* lines longer than this are not drawn (useful for Sinai billiard) */ +#define LCUT 10000000.0 /* controls the max size of segments not considered as being cut */ +#define DMIN 0.1 /* 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 */ +#define ORDER_COLORS 0 /* set to 1 if colors should be drawn in order */ /* color and other graphical parameters */ -#define COLOR_PALETTE 1 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ -#define NCOLORS 14 /* number of colors */ +#define NCOLORS 16 /* 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 0 /* set to 1 for black background */ +#define BLACK 1 /* set to 1 for black background */ #define COLOR_OUTSIDE 0 /* set to 1 for colored outside */ #define OUTER_COLOR 300.0 /* color outside billiard */ #define PAINT_INT 0 /* set to 1 to paint interior in other color (for polygon) */ @@ -119,7 +120,7 @@ #define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 100 /* final sleeping time */ -#define END_FRAMES 0 /* number of frames at end of movie */ +#define END_FRAMES 25 /* number of frames at end of movie */ #define PI 3.141592654 #define DPI 6.283185307 @@ -469,11 +470,11 @@ void animation() for (i=0; i INITIAL_TIME + NSTEPS - FINAL_CONSTANT_PHASE) beta = BETA*BETA_FACTOR; - else + else if (i < INITIAL_TIME + t1) { beta = BETA*exp(bexponent*(double)(i - INITIAL_TIME)); if (!NO_OSCILLATION) beta = beta*2.0/(1.0 + cos(omega*(double)(i - INITIAL_TIME))); } + else if (i < INITIAL_TIME + t2) beta = BETA*BETA_FACTOR; + else if (i < INITIAL_TIME + t3) + { + beta = BETA*exp(bexp2*(double)(i - INITIAL_TIME - t3)); + } + else beta = BETA; printf("beta = %.3lg\n", beta); return(beta); } @@ -720,10 +740,12 @@ void evolve_wall(double fboundary) } -void evolve_segments(t_segment segment[NMAXSEGMENTS]) +void evolve_segments(t_segment segment[NMAXSEGMENTS], int time) { int i, nactive = 0, group; - double fx[2] = {0.0, 0.0}, fy[2] = {0.0, 0.0}, x, y, padding = 3.0*MU, mass2 = SEGMENTS_MASS*TWO_CIRCLES_RADIUS_RATIO; + double fx[2] = {0.0, 0.0}, fy[2] = {0.0, 0.0}, x, y, padding = 3.0*MU, mass2 = SEGMENTS_MASS; + + if (SEGMENT_PATTERN == S_TWO_CIRCLES_EXT) mass2 = SEGMENTS_MASS*TWO_CIRCLES_RADIUS_RATIO; for (group=0; group<2; group++) { @@ -745,6 +767,11 @@ void evolve_segments(t_segment segment[NMAXSEGMENTS]) if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y); else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding); } + else if (BOUNDARY_COND == BC_REFLECT_ABS) /* add force from simulation boundary */ + { + y = 0.5*(segment[i].y1 + segment[i].y2); + if (y < YMIN) fy[group] += KSPRING_BOUNDARY*(YMIN - y); + } if (group == 0) fy[group] -= GRAVITY*SEGMENTS_MASS; else fy[group] -= GRAVITY*mass2; } @@ -778,12 +805,31 @@ void evolve_segments(t_segment segment[NMAXSEGMENTS]) xsegments[1] += vxsegments[1]*DT_PARTICLE; ysegments[1] += vysegments[1]*DT_PARTICLE; } + + /* add some damping if y coordinate is small (for lunar landing) */ + if (DAMP_SEGS_AT_NEGATIVE_Y) + for (group=0; group<2; group++) + if (ysegments[group] < 0.1) + { + vysegments[group] *= exp(-DAMPING*DT_PARTICLE); + vxsegments[group] *= exp(-DAMPING*DT_PARTICLE); + } /* to avoid numerical instabilities */ - for (group=0; group<2; group++) if (xsegments[group] + 1.0 > BCXMAX) + for (group=0; group<2; group++) { - xsegments[group] = BCXMAX - 1.0; - vxsegments[group] = 0.0; + if (xsegments[group] + 1.0 > BCXMAX) + { + xsegments[group] = BCXMAX - 1.0; + vxsegments[group] = 0.0; + } + if ((RELEASE_ROCKET_AT_DEACTIVATION)&&((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_ABSORBING))) + { +// ysegments[group] = SEGMENTS_Y0; + if (time < SEGMENT_DEACTIVATION_TIME) vysegments[group] = 0.0; + else if ((ysegments[group] < SEGMENTS_Y0)&&(vysegments[group] < 0.0)) + vysegments[group] = -0.5*vysegments[group]; + } } translate_segments(segment, xsegments, ysegments); @@ -876,8 +922,9 @@ void animation() thermostat_on = thermostat_schedule(i); printf("Termostat: %i\n", thermostat_on); } - if ((ADD_FIXED_SEGMENTS)&&(DEACTIVATE_SEGMENT)&&(i > INITIAL_TIME + SEGMENT_DEACTIVATION_TIME)) - segment[nsegments-1].active = 0; + /* deactivate some segments */ + if ((ADD_FIXED_SEGMENTS)&&(DEACTIVATE_SEGMENT)&&(i == INITIAL_TIME + SEGMENT_DEACTIVATION_TIME + 1)) + for (j=0; j FMAX) particle[j].fx = FMAX; @@ -947,7 +998,7 @@ void animation() if (i < INITIAL_TIME + WALL_TIME) evolve_wall(fboundary); else xwall = 0.0; } - if ((MOVE_BOUNDARY)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segments(segment); + if ((MOVE_BOUNDARY)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segments(segment, i); } /* end of for (n=0; n= INITIAL_TIME) { - save_frame_lj(); - if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) + if (TIME_LAPSE_FIRST) { - save_frame_lj_counter(NSTEPS + END_FRAMES + (i - INITIAL_TIME)/TIME_LAPSE_FACTOR); + if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) + { + save_frame_lj(); + } + save_frame_lj_counter(NSTEPS/TIME_LAPSE_FACTOR + MID_FRAMES + i - INITIAL_TIME); + } + else + { + save_frame_lj(); + if ((TIME_LAPSE)&&((i - INITIAL_TIME)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) + { + save_frame_lj_counter(NSTEPS + END_FRAMES + (i - INITIAL_TIME)/TIME_LAPSE_FACTOR); + } } } else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); @@ -1120,10 +1182,13 @@ void animation() else if (PRINT_SEGMENTS_SPEEDS) print_segments_speeds(vxsegments, vysegments); glutSwapBuffers(); } - for (i=0; i #include -#define MOVIE 0 /* set to 1 to generate movie */ +#define MOVIE 1 /* set to 1 to generate movie */ /* General geometrical parameters */ @@ -68,25 +68,22 @@ /* Boundary conditions, see list in global_pdes.c */ -#define LATTICE 10 +#define LATTICE 2 -#define FLOOD_LEFT_BOUNDARY 1 /* set to 1 to flood cells on left boundary */ -#define FIND_ALL_CLUSTERS 0 /* set to 1 to find all open clusters */ +#define FLOOD_LEFT_BOUNDARY 0 /* set to 1 to flood cells on left boundary */ +#define FIND_ALL_CLUSTERS 1 /* set to 1 to find all open clusters */ -#define PLOT_CLUSTER_SIZE 1 /* set to 1 to add a plot for the size of the percolation cluster */ +#define PLOT_CLUSTER_SIZE 0 /* set to 1 to add a plot for the size of the percolation cluster */ #define PLOT_CLUSTER_NUMBER 0 /* set to 1 to add a graph of the number of clusters */ -#define PLOT_CLUSTER_HISTOGRAM 0 /* set to 1 to add a histogram of the number of clusters */ -#define PRINT_LARGEST_CLUSTER_SIZE 0 /* set to 1 to print size of largest cluster */ +#define PLOT_CLUSTER_HISTOGRAM 1 /* set to 1 to add a histogram of the number of clusters */ +#define PRINT_LARGEST_CLUSTER_SIZE 1 /* set to 1 to print size of largest cluster */ #define MAX_CLUSTER_NUMBER 6 /* vertical scale of the cluster number plot */ -#define HISTO_BINS 30 /* number of bins in histrogram */ +#define HISTO_BINS 30 /* number of bins in histogram */ -#define NSTEPS 1270 /* number of frames of movie */ - -#define DEBUG 0 /* set to 1 for some debugging features */ -#define DEBUG_SLEEP_TIME 1 /* sleep time between frames when debugging */ -#define TEST_GRAPH 0 /* set to 1 to test graph connectivity matrix */ -#define TEST_START 2210 /* start position of connectivity test */ +#define NSTEPS 100 /* number of frames of movie */ +// #define NSTEPS 700 /* number of frames of movie */ +// #define NSTEPS 830 /* number of frames of movie */ #define PAUSE 200 /* number of frames after which to pause */ #define PSLEEP 2 /* sleep time during pause */ @@ -102,19 +99,19 @@ #define BLACK 1 /* background */ -#define COLOR_CLUSTERS_BY_SIZE 0 /* set to 1 to link cluster color to their size */ +#define COLOR_CLUSTERS_BY_SIZE 1 /* set to 1 to link cluster color to their size */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define HUE_CLOSED 330.0 /* color hue of closed cells */ +#define HUE_CLOSED 350.0 /* color hue of closed cells */ #define HUE_OPEN 200.0 /* color hue of open (dry) cells */ #define HUE_FLOODED 45.0 /* color hue of open flooded cells */ -#define HUE_GRAPH 150.0 /* color hue in graph of cluster size */ +#define HUE_GRAPH 250.0 /* color hue in graph of cluster size */ -#define CLUSTER_HUEMIN 0.0 /* minimal color hue of clusters */ -#define CLUSTER_HUEMAX 250.0 /* maximal color hue of clusters */ +#define CLUSTER_HUEMIN 60.0 /* minimal color hue of clusters */ +#define CLUSTER_HUEMAX 300.0 /* maximal color hue of clusters */ #define N_CLUSTER_COLORS 20 /* number of different colors of clusters */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -124,6 +121,13 @@ #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 */ +/* debugging options */ +#define VERBOSE 0 /* set to 1 to print more messages in shell */ +#define DEBUG 0 /* set to 1 for some debugging features */ +#define DEBUG_SLEEP_TIME 1 /* sleep time between frames when debugging */ +#define TEST_GRAPH 0 /* set to 1 to test graph connectivity matrix */ +#define TEST_START 2210 /* start position of connectivity test */ + #define ADD_PLOT ((PLOT_CLUSTER_SIZE)||(PLOT_CLUSTER_NUMBER)||(PLOT_CLUSTER_HISTOGRAM)) #define FIND_CLUSTER_SIZES ((COLOR_CLUSTERS_BY_SIZE)||(PLOT_CLUSTER_HISTOGRAM)) @@ -134,8 +138,9 @@ void animation(int size) { - int i, j, k, s, nx, ny, ncells, nmaxcells, maxsize, nopen, nflooded, nstack, nclusters, maxclustersize = 0, maxclusterlabel; + int i, j, k, s, nx, ny, nmaxcells, maxsize, nopen, nflooded, nstack, nclusters, maxclustersize = 0, maxclusterlabel; int *plot_cluster_number, *cluster_sizes; + int ncells; double p, *plot_cluster_size; t_perco *cell; t_perco **pstack; @@ -147,10 +152,12 @@ void animation(int size) cell = (t_perco *)malloc(nmaxcells*sizeof(t_perco)); if (PLOT_CLUSTER_SIZE) plot_cluster_size = (double *)malloc(NSTEPS*sizeof(double)); if (PLOT_CLUSTER_NUMBER) plot_cluster_number = (int *)malloc(NSTEPS*sizeof(double)); - if (FIND_CLUSTER_SIZES) cluster_sizes = (int *)malloc(2*nmaxcells*sizeof(int)); +// if (FIND_CLUSTER_SIZES) + cluster_sizes = (int *)malloc(2*nmaxcells*sizeof(int)); ncells = init_cell_lattice(cell, nx, ny); - printf("nx = %i, ny = %i, ncells = %i\n", nx, ny, ncells); + printf("nx = %i, ny = %i, ncells = %i, maxcells = %i\n", nx, ny, ncells, nmaxcells); + pstack = (t_perco* *)malloc(ncells*sizeof(struct t_perco *)); if (TEST_GRAPH) test_neighbours(TEST_START, cell, nx, ny, size, ncells); @@ -164,11 +171,17 @@ void animation(int size) init_cell_state(cell, p, ncells, (i == 0)); - if (FLOOD_LEFT_BOUNDARY) nstack = init_flooded_cells(cell, nx, ny, pstack); + if (FLOOD_LEFT_BOUNDARY) nstack = init_flooded_cells(cell, ncells, nx, ny, pstack); nopen = count_open_cells(cell, ncells); - if (FLOOD_LEFT_BOUNDARY) nflooded = find_percolation_cluster(cell, ncells, pstack, nstack); - + printf("Flooded cells, %i open cells, nstack = %i\n", nopen, nstack); + + if (FLOOD_LEFT_BOUNDARY) + { + nflooded = find_percolation_cluster(cell, ncells, pstack, nstack); + printf("Found percolation cluster with %i flooded cells\n", nflooded); + } + if (FIND_ALL_CLUSTERS) { nclusters = find_all_clusters(cell, ncells, (i == 0), &maxclusterlabel); @@ -183,7 +196,7 @@ void animation(int size) // print_cluster_sizes(cell, ncells, cluster_sizes); - draw_configuration(cell, cluster_sizes, nx, ny, size, ncells); + draw_configuration(cell, cluster_sizes, ncells, nx, ny, size, ncells); print_p(p); if (PRINT_LARGEST_CLUSTER_SIZE) print_largest_cluster_size(maxclustersize); @@ -226,7 +239,8 @@ void animation(int size) if (PLOT_CLUSTER_SIZE) free(plot_cluster_size); if (PLOT_CLUSTER_NUMBER) free(plot_cluster_number); - if (FIND_CLUSTER_SIZES) free(cluster_sizes); +// if (FIND_CLUSTER_SIZES) + free(cluster_sizes); } @@ -249,12 +263,12 @@ void display(void) // animation(128); // animation(64); - animation(32); - animation(16); - animation(8); +// animation(32); +// animation(16); +// animation(8); animation(4); - animation(2); - animation(1); +// animation(2); +// animation(1); sleep(SLEEP2); diff --git a/rde.c b/rde.c index cfeeba7..e1e8914 100644 --- a/rde.c +++ b/rde.c @@ -39,54 +39,66 @@ #include #include -#define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define MOVIE 1 /* set to 1 to generate movie */ +#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ /* General geometrical parameters */ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1000 /* window height */ -#define NX 480 /* number of grid points on x axis */ -#define NY 240 /* number of grid points on y axis */ - -#define XMIN -2.0 -#define XMAX 2.0 /* x interval */ -#define YMIN -1.041666667 -#define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ - -// #define WINWIDTH 1280 /* window width */ -// #define WINHEIGHT 720 /* window height */ -// -// // #define NX 160 /* number of grid points on x axis */ -// // #define NY 90 /* number of grid points on y axis */ -// #define NX 320 /* number of grid points on x axis */ -// #define NY 180 /* 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 1280 /* number of grid points on x axis */ -// // #define NY 720 /* number of grid points on y axis */ +// #define NX 600 /* number of grid points on x axis */ +// #define NY 300 /* number of grid points on y axis */ +// // #define NX 480 /* number of grid points on x axis */ +// // #define NY 240 /* number of grid points on y axis */ +// // #define NX 1920 /* number of grid points on x axis */ +// // #define NY 1000 /* number of grid points on y axis */ // // #define XMIN -2.0 // #define XMAX 2.0 /* x interval */ -// #define YMIN -1.125 -// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +// #define YMIN -1.041666667 +// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ + +// #define WINWIDTH 1280 /* window width */ +// #define WINHEIGHT 720 /* window height */ + +// #define NX 200 /* number of grid points on x axis */ +// #define NY 200 /* number of grid points on y axis */ +#define NX 500 /* number of grid points on x axis */ +#define NY 500 /* 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 1280 /* number of grid points on x axis */ +// #define NY 720 /* number of grid points on y axis */ + +#define XMIN -1.8 +#define XMAX 1.8 /* x interval */ +#define YMIN -1.8 +#define YMAX 1.8 /* y interval for 9/16 aspect ratio */ /* Choice of simulated equation */ #define RDE_EQUATION 5 /* choice of reaction term, see list in global_3d.c */ #define NFIELDS 2 /* number of fields in reaction-diffusion equation */ #define NLAPLACIANS 2 /* number of fields for which to compute Laplacian */ +// #define RDE_EQUATION 4 /* choice of reaction term, see list in global_3d.c */ +// #define NFIELDS 3 /* number of fields in reaction-diffusion equation */ +// #define NLAPLACIANS 3 /* number of fields for which to compute Laplacian */ -#define ADD_POTENTIAL 1 /* set to 1 to add a potential (for Schrodiner equation) */ -#define POTENTIAL 2 /* type of potential, see list in global_3d.c */ +#define ADD_POTENTIAL 1 /* set to 1 to add a potential (for Schrodinger equation) */ +#define POTENTIAL 1 /* type of potential, see list in global_3d.c */ +#define ADD_MAGNETIC_FIELD 1 /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */ + +#define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */ #define JULIA_SCALE 0.5 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 99 /* pattern of circles, see list in global_pdes.c */ @@ -96,8 +108,8 @@ #define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ #define MU 1.0 /* parameter controlling the dimensions of domain */ -#define NPOLY 6 /* number of sides of polygon */ -#define APOLY 0.333333333 /* angle by which to turn polygon, in units of Pi/2 */ +#define NPOLY 5 /* number of sides of polygon */ +#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 7 /* depth of computation of Menger gasket */ #define MRATIO 5 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ @@ -111,7 +123,7 @@ #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_LEFT -1.65 #define ISO_XSHIFT_RIGHT 0.4 #define ISO_YSHIFT_LEFT -0.05 #define ISO_YSHIFT_RIGHT -0.05 @@ -123,6 +135,9 @@ /* Physical patameters of wave equation */ #define DT 0.00000002 +// #define DT 0.00000003 +// #define DT 0.000000011 +// #define DT 0.00000001 #define VISCOSITY 2.0 @@ -133,7 +148,7 @@ #define DELTA 0.1 /* time scale separation */ #define FHNA 1.0 /* parameter in FHN equation */ #define FHNC -0.01 /* parameter in FHN equation */ -#define K_HARMONIC 0.2 /* spring constant of harmonic potential */ +#define K_HARMONIC 0.5 /* spring constant of harmonic potential */ #define K_COULOMB 0.5 /* constant in Coulomb potential */ #define BZQ 0.0008 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */ @@ -167,13 +182,15 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1150 /* number of frames of movie */ -#define NVID 850 /* number of iterations between images displayed on screen */ +// #define NSTEPS 500 /* number of frames of movie */ +#define NSTEPS 1100 /* number of frames of movie */ +#define NVID 500 /* number of iterations between images displayed on screen */ +// #define NVID 1100 /* number of iterations between images displayed on screen */ #define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */ #define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */ #define MAX_DT 0.024 /* maximal value of integration step */ #define NSEG 100 /* number of segments of boundary */ -#define BOUNDARY_WIDTH 4 /* width of billiard boundary */ +#define BOUNDARY_WIDTH 5 /* width of billiard boundary */ #define PAUSE 100 /* number of frames after which to pause */ #define PSLEEP 2 /* sleep time during pause */ @@ -188,19 +205,19 @@ #define PLOT_3D 1 /* controls whether plot is 2D or 3D */ -#define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ +#define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ /* Plot type - color scheme */ -#define CPLOT 30 +#define CPLOT 32 +// #define CPLOT 32 #define CPLOT_B 31 /* Plot type - height of 3D plot */ -// #define ZPLOT 30 /* z coordinate in 3D plot */ -// #define ZPLOT_B 32 /* z coordinate in second 3D plot */ -#define ZPLOT 30 /* z coordinate in 3D plot */ +#define ZPLOT 32 /* z coordinate in 3D plot */ +// #define ZPLOT 32 /* z coordinate in 3D plot */ #define ZPLOT_B 30 /* z coordinate in second 3D plot */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */ @@ -209,8 +226,8 @@ #define WRAP_ANGLE 1 /* experimental: wrap angle to [0, 2Pi) for interpolation in angle schemes */ #define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ #define DRAW_OUTSIDE_GRAY 0 /* experimental - draw outside of billiard in gray */ -#define ADD_POTENTIAL_TO_Z 1 /* set to 1 to add the external potential to z-coordinate of plot */ -#define ADD_POT_CONSTANT 0.5 /* constant in front of added potential */ +#define ADD_POTENTIAL_TO_Z 0 /* set to 1 to add the external potential to z-coordinate of plot */ +#define ADD_POT_CONSTANT 1.0 /* constant in front of added potential */ #define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */ @@ -237,8 +254,8 @@ /* Color schemes, see list in global_pdes.c */ -#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 0 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ @@ -248,12 +265,12 @@ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ -#define VSCALE_AMPLITUDE 7.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_AMPLITUDE 30.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define CURL_SCALE 0.000015 /* scaling factor for curl representation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ -#define SLOPE_SCHROD_LUM 15.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ -#define MIN_SCHROD_LUM 0.075 /* minimal luminosity in color scheme Z_ARGUMENT*/ +#define SLOPE_SCHROD_LUM 50.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ +#define MIN_SCHROD_LUM 0.2 /* minimal luminosity in color scheme Z_ARGUMENT*/ #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 */ @@ -263,11 +280,11 @@ #define HUEAMP -359.0 /* amplitude of variation of hue for color scheme C_HUE */ #define E_SCALE 100.0 /* scaling factor for energy representation */ #define LOG_SCALE 1.0 /* scaling factor for energy log representation */ -#define LOG_SHIFT 0.0 +#define LOG_SHIFT 0.0 #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 3.0 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ /* only for compatibility with wave_common.c */ @@ -281,6 +298,13 @@ #define VSCALE_ENERGY 200.0 /* additional scaling factor for color scheme P_3D_ENERGY */ #define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ +#define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */ +#define AMPLITUDE 0.8 /* amplitude of periodic excitation */ +#define ACHIRP 0.2 /* acceleration coefficient in chirp */ +#define DAMPING 0.0 /* damping of periodic excitation */ +#define COMPARISON 0 /* set to 1 to compare two different patterns (beta) */ +#define B_DOMAIN_B 20 /* second domain shape, for comparisons */ +#define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */ /* end of constants added only for compatibility with wave_common.c */ @@ -288,14 +312,14 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {8.0, 8.0, 6.0}; /* location of observer for REP_PROJ_3D representation */ +double observer[3] = {8.0, 8.0, 12.0}; /* location of observer for REP_PROJ_3D representation */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ -#define Z_SCALING_FACTOR 0.75 /* overall scaling factor of z axis for REP_PROJ_3D representation */ -#define XY_SCALING_FACTOR 2.2 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define Z_SCALING_FACTOR 1.25 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define XY_SCALING_FACTOR 1.8 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ #define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D 0.2 /* overall y shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */ /* For debugging purposes only */ @@ -319,7 +343,7 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO double potential(int i, int j) /* compute potential (e.g. for Schrödinger equation) */ { - double x, y, xy[2], r, small = 2.0e-1, kx, ky; + double x, y, xy[2], r, small = 2.0e-1, kx, ky, lx = XMAX - XMIN, r1, r2, r3; ij_to_xy(i, j, xy); x = xy[0]; @@ -343,12 +367,42 @@ double potential(int i, int j) ky = 2.0*DPI/(YMAX - YMIN); return(-K_HARMONIC*cos(kx*x)*cos(ky*y)); } + case (POT_FERMIONS): + { + r = sqrt((x-y)*(x-y) + small*small); + return (-K_COULOMB/r); + } + case (POT_FERMIONS_PERIODIC): + { + r1 = sqrt((x-y)*(x-y) + small*small); + r2 = sqrt((x-lx-y)*(x-lx-y) + small*small); + r3 = sqrt((x+lx-y)*(x+lx-y) + small*small); +// r = r/3.0; + return (-0.5*K_COULOMB*(1.0/r1 + 1.0/r2 + 1.0/r3)); + } default: { return(0.0); } } -} +} + + +void compute_vector_potential(int i, int j, double *ax, double *ay) +/* initialize the vector potential, for Schrodinger equation in a magnetic field */ +{ + double x, y, xy[2], b; + + ij_to_xy(i, j, xy); + x = xy[0]; + y = xy[1]; + + b = sqrt(K_HARMONIC); + /* magnetic field strength b is chosen such that b^2/4 = K_HARMONIC */ + + *ax = b*y; + *ay = -b*x; +} void initialize_potential(double potential_field[NX*NY]) @@ -364,18 +418,40 @@ void initialize_potential(double potential_field[NX*NY]) } } -void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY]) +void initialize_vector_potential(double vpotential_field[2*NX*NY]) +/* initialize the potential field, e.g. for the Schrödinger equation */ +{ + int i, j; + + #pragma omp parallel for private(i,j) + for (i=0; i INITYMIN - MU)) particles[n].active = 1; + if ((particles[n].yc < INITYMAX + MU)&&(particles[n].yc > INITYMIN - MU)&&(particles[n].xc < INITXMAX + MU)&&(particles[n].xc > INITXMIN - MU)) particles[n].active = 1; else particles[n].active = 0; } break; @@ -1030,20 +1031,92 @@ void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, d else printf("Warning: NMAXSEGMENTS too small\n"); } -double nozzle_width(double x, double width) +double nozzle_width(double x, double width, int nozzle_shape) /* width of bell-shaped nozzle */ { double lam = 0.5*LAMBDA; if (x >= 0.0) return(width); - else return(sqrt(width*width - 0.5*x)); + else switch (nozzle_shape) { + case (NZ_STRAIGHT): return(width); + case (NZ_BELL): return(sqrt(width*width - 0.5*x)); + case (NZ_GLAS): return(sqrt(width*width - 1.2*x) + 1.0*x); + case (NZ_CONE): return(width - (sqrt(width*width + 0.5) - width)*x); + case (NZ_TRUMPET): return(width + (sqrt(width*width + LAMBDA)-width)*x*x); + default: return(0.0); + } +} + + +void add_rocket_to_segments(t_segment segment[NMAXSEGMENTS], double x0, double y0, int nozzle_shape, int nsides, int group) +/* add one or several rocket_shaped set of segments */ +{ + int i, j, cycle = 0, nsegments0; + double angle, dx, x1, y1, x2, y2, nozx, nozy; + + nsegments0 = nsegments; + + /* ellipse */ + for (i=1; i 7) segment[i].inactivate = 1; + else segment[i].inactivate = 0; + } + break; + } + case (S_DAM_WITH_HOLE_AND_RAMP): + { + add_rectangle_to_segments(DAM_WIDTH, BCYMIN - 0.5, -DAM_WIDTH, BCYMIN + 0.2, segment); + add_rectangle_to_segments(DAM_WIDTH, BCYMIN + 0.3, -DAM_WIDTH, LAMBDA, segment); + add_rectangle_to_segments(DAM_WIDTH, BCYMIN + 0.2, -DAM_WIDTH, BCYMIN + 0.3, segment); + + r = 1.0; + for (i=0; i<10; i++) + { + angle = 0.1*PID*(double)i; + dangle = 0.1*PID; + x1 = XMAX - r + (r + MU)*cos(angle); + y1 = YMIN + r - (r + MU)*sin(angle); + x2 = XMAX - r + (r + MU)*cos(angle + dangle); + y2 = YMIN + r - (r + MU)*sin(angle + dangle); + add_rotated_angle_to_segments(x1, y1, x2, y2, MU, segment); + } + + cycle = 0; + for (i=0; i 7)&&(i < 12)) segment[i].inactivate = 1; + else segment[i].inactivate = 0; + } + break; + } default: { printf("Function init_segment_config not defined for this pattern \n"); @@ -1458,6 +1628,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS]) segment[nsegments].y1 = LAMBDA*sin(((double)i + 0.5)*angle); segment[nsegments].x2 = -cos(((double)i + 0.5)*angle); segment[nsegments].y2 = -LAMBDA*sin(((double)i + 0.5)*angle); + segment[nsegments].inactivate = 1; nsegments++; } @@ -1530,7 +1701,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS]) int in_segment_region(double x, double y) /* returns 1 if (x,y) is inside region delimited by obstacle segments */ { - double angle, dx, height, width, theta, lx, ly; + double angle, dx, height, width, theta, lx, ly, x1, y1, x2, y2; if (x >= BCXMAX) return(0); if (x <= BCXMIN) return(0); @@ -1622,12 +1793,52 @@ int in_segment_region(double x, double y) else { lx = 0.7*LAMBDA; - ly = 0.5*LAMBDA; + ly = 0.7*LAMBDA; x -= lx; if (x*x/(lx*lx) + y*y/(ly*ly) < 0.95) return(1); } return(0); } + case (S_ROCKET_NOZZLE_ROTATED): + { + y1 = y - ysegments[0]; + if (y1 < YMIN + LAMBDA) return(0); + else if (y1 > YMIN + 2.4*LAMBDA) return(0); + else + { + ly = 0.7*LAMBDA; + lx = 0.7*LAMBDA; + x1 = x - xsegments[0]; + y1 -= YMIN + 1.7*LAMBDA; + if (x1*x1/(lx*lx) + y1*y1/(ly*ly) + MU*MU < 0.925) return(1); + } + return(0); + } + case (S_TWO_ROCKETS): + { + y1 = y - ysegments[0]; + y2 = y - ysegments[1]; + if ((y1 < YMIN + LAMBDA)&&(y2 < YMIN + LAMBDA)) return(0); + else if ((y1 > YMIN + 2.4*LAMBDA)&&(y2 > YMIN + 2.4*LAMBDA)) return(0); + else + { + ly = 0.7*LAMBDA; + lx = 0.7*LAMBDA; + x1 = x - xsegments[0]; + y1 -= YMIN + 1.7*LAMBDA; + if (x1*x1/(lx*lx) + y1*y1/(ly*ly) + MU*MU < 0.925) return(1); + x2 = x - xsegments[1]; + y2 -= YMIN + 1.7*LAMBDA; + if (x2*x2/(lx*lx) + y2*y2/(ly*ly) + MU*MU < 0.925) return(1); + } + return(0); + } + case (S_DAM): + { + if (vabs(x) > DAM_WIDTH) return(1); + else if (y > LAMBDA) return(1); + else return(0); + } default: return(1); } } @@ -2281,12 +2492,12 @@ int compute_repelling_force(int i, int j, double force[2], double *torque, t_par distance = module2(particle[i].deltax[j], particle[i].deltay[j]); /* for monitoring purposes */ - if (distance > distmin) - { - printf("i = %i, hashcell %i, j = %i, hashcell %i\n", i, particle[i].hashcell, k, particle[k].hashcell); - printf("X = (%.3lg, %.3lg)\n", particle[i].xc, particle[i].yc); - printf("Y = (%.3lg, %.3lg) d = %.3lg\n", particle[k].xc, particle[k].yc, distance); - } +// if (distance > distmin) +// { +// printf("i = %i, hashcell %i, j = %i, hashcell %i\n", i, particle[i].hashcell, k, particle[k].hashcell); +// printf("X = (%.3lg, %.3lg)\n", particle[i].xc, particle[i].yc); +// printf("Y = (%.3lg, %.3lg) d = %.3lg\n", particle[k].xc, particle[k].yc, distance); +// } if ((distance == 0.0)||(i == k)) { @@ -2519,11 +2730,202 @@ void compute_entropy(t_particle particle[NMAXCIRCLES], double entropy[2]) } -void draw_triangles(t_particle particle[NMAXCIRCLES]) +void compute_particle_colors(t_particle particle, int plot, double rgb[3], double rgbx[3], double rgby[3], double *radius, int *width) +{ + double ej, angle, hue, huex, huey, lum; + int i; + + switch (plot) { + case (P_KINETIC): + { + ej = particle.energy; + if (ej > 0.0) + { + hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej/PARTICLE_EMAX; + if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; + if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; + } + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_NEIGHBOURS): + { + hue = neighbour_color(particle.neighb); + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_BONDS): + { + hue = neighbour_color(particle.neighb); + *radius = particle.radius; + *width = 1; + break; + } + case (P_ANGLE): + { + angle = particle.angle; + hue = angle*particle.spin_freq/DPI; + hue -= (double)((int)hue); + huex = (DPI - angle)*particle.spin_freq/DPI; + huex -= (double)((int)huex); + angle = PI - angle; + if (angle < 0.0) angle += DPI; + huey = angle*particle.spin_freq/DPI; + huey -= (double)((int)huey); + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue); + huex = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(huex); + huey = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(huey); + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_TYPE): + { + if (particle.type <= 1) hue = HUE_TYPE0; + else if (particle.type == 2) hue = HUE_TYPE1; + else if (particle.type == 3) hue = HUE_TYPE2; + else hue = HUE_TYPE3; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_DIRECTION): + { + hue = argument(particle.vx, particle.vy); + if (hue > DPI) hue -= DPI; + if (hue < 0.0) hue += DPI; + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_DIRECT_ENERGY): + { + hue = argument(particle.vx, particle.vy); + if (hue > DPI) hue -= DPI; + if (hue < 0.0) hue += DPI; + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; + if (particle.energy < 0.1*PARTICLE_EMAX) lum = 10.0*particle.energy/PARTICLE_EMAX; + else lum = 1.0; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_ANGULAR_SPEED): + { + hue = 160.0*(1.0 + tanh(SLOPE*particle.omega)); +// printf("omega = %.3lg, hue = %.3lg\n", particle[j].omega, hue); + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_DIFF_NEIGHB): + { + hue = (double)(particle.diff_neighb+1)/(double)(particle.neighb+1); +// hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; + hue = 180.0*(1.0 + hue); + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_THERMOSTAT): + { + if (particle.thermostat) hue = 30.0; + else hue = 270.0; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_INITIAL_POS): + { + hue = particle.color_hue; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + } + + switch (plot) { + case (P_KINETIC): + { + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); + break; + } + case (P_BONDS): + { + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); + break; + } + case (P_DIRECTION): + { + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); + break; + } + case (P_DIRECT_ENERGY): + { + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); + for (i=0; i<3; i++) + { + rgb[i] *= lum; + rgbx[i] *= lum; + rgby[i] *= lum; + } + break; + } + case (P_DIFF_NEIGHB): + { + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); + hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); + break; + } + default: + { + hsl_to_rgb(hue, 0.9, 0.5, rgb); + hsl_to_rgb(hue, 0.9, 0.5, rgbx); + hsl_to_rgb(hue, 0.9, 0.5, rgby); + } + } +} + + +void draw_one_triangle(t_particle particle, int same_table[9*HASHMAX], int p, int q, int nsame) +{ + double x, y, dx, dy; + int k; + + x = particle.xc + (double)p*(BCXMAX - BCXMIN); + y = particle.yc + (double)q*(BCYMAX - BCYMIN); + + glBegin(GL_TRIANGLE_FAN); + glVertex2d(x, y); + + for (k=0; k 1) { - if (t0 <= 1) hue = HUE_TYPE0; - else if (t0 == 2) hue = HUE_TYPE1; - else if (t0 == 3) hue = HUE_TYPE2; - else hue = HUE_TYPE3; - - hsl_to_rgb(hue, 0.9, 0.3, rgb); - glColor3f(rgb[0], rgb[1], rgb[2]); - for (p=-1; p<2; p++) for (q=-1; q<2; q++) + if (plot == P_TYPE) { - x = particle[i].xc + (double)p*(BCXMAX - BCXMIN); - y = particle[i].yc + (double)q*(BCYMAX - BCYMIN); - - glBegin(GL_TRIANGLE_FAN); - glVertex2d(x, y); -// printf("%i | ", particle[i].hash_nneighb); - for (k=0; k 0.0) - { - hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej/PARTICLE_EMAX; - if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; - if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; - } - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_NEIGHBOURS): - { - hue = neighbour_color(particle[j].neighb); - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_BONDS): - { -// if (particle[j].type == 1) hue = 70.0; /* to make second particle type more visible */ -// if (particle[j].type == 1) hue = neighbour_color(7 - particle[j].neighb); -// else - hue = neighbour_color(particle[j].neighb); - radius = particle[j].radius; - width = 1; - break; - } - case (P_ANGLE): - { - angle = particle[j].angle; - hue = angle*particle[j].spin_freq/DPI; - hue -= (double)((int)hue); - huex = (DPI - angle)*particle[j].spin_freq/DPI; - huex -= (double)((int)huex); - angle = PI - angle; - if (angle < 0.0) angle += DPI; - huey = angle*particle[j].spin_freq/DPI; - huey -= (double)((int)huey); - hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; - huex = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*huex; - huey = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*huey; - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_TYPE): - { -// if (particle[j].type == 0) hue = 310.0; -// else hue = 70.0; - if (particle[j].type <= 1) hue = HUE_TYPE0; - else if (particle[j].type == 2) hue = HUE_TYPE1; - else if (particle[j].type == 3) hue = HUE_TYPE2; - else hue = HUE_TYPE3; - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_DIRECTION): - { - hue = argument(particle[j].vx, particle[j].vy); - if (hue > DPI) hue -= DPI; - if (hue < 0.0) hue += DPI; - hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue/DPI; - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_DIRECT_ENERGY): - { - hue = argument(particle[j].vx, particle[j].vy); - if (hue > DPI) hue -= DPI; - if (hue < 0.0) hue += DPI; - hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue/DPI; - if (particle[j].energy < 0.1*PARTICLE_EMAX) lum = 10.0*particle[j].energy/PARTICLE_EMAX; - else lum = 1.0; - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_ANGULAR_SPEED): - { - hue = 160.0*(1.0 + tanh(SLOPE*particle[j].omega)); -// printf("omega = %.3lg, hue = %.3lg\n", particle[j].omega, hue); - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - case (P_DIFF_NEIGHB): - { - hue = (double)(particle[j].diff_neighb+1)/(double)(particle[j].neighb+1); -// hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; - hue = 180.0*(1.0 + hue); - radius = particle[j].radius; - width = BOUNDARY_WIDTH; - break; - } - } - + compute_particle_colors(particle[j], plot, rgb, rgbx, rgby, &radius, &width); + switch (particle[j].interaction) { case (I_LJ_DIRECTIONAL): { @@ -2891,55 +3186,6 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta) default: nsides = NSEG; } - switch (plot) { - case (P_KINETIC): - { - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); - break; - } - case (P_BONDS): - { - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgbx); - hsl_to_rgb_turbo(hue, 0.9, 0.5, rgby); - break; - } - case (P_DIRECTION): - { - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); - break; - } - case (P_DIRECT_ENERGY): - { - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); - for (i=0; i<3; i++) - { - rgb[i] *= lum; - rgbx[i] *= lum; - rgby[i] *= lum; - } - break; - } - case (P_DIFF_NEIGHB): - { - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb); - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx); - hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby); - break; - } - default: - { - hsl_to_rgb(hue, 0.9, 0.5, rgb); - hsl_to_rgb(hue, 0.9, 0.5, rgbx); - hsl_to_rgb(hue, 0.9, 0.5, rgby); - } - } angle = particle[j].angle + APOLY*DPI; draw_one_particle(particle[j], particle[j].xc, particle[j].yc, radius, angle, nsides, width, rgb); @@ -3279,27 +3525,29 @@ void print_parameters(double beta, double temperature, double krepel, double len int i, j, k; double density, hue, rgb[3], logratio, x, y, meanpress[N_PRESSURES], phi, sphi, dphi, pprint, mean_temp; static double xbox, xtext, xmid, xmidtext, xxbox, xxtext, pressures[N_P_AVERAGE], meanpressure = 0.0, maxpressure = 0.0; - static double press[N_PRESSURES][N_P_AVERAGE], temp[N_T_AVERAGE]; + static double press[N_PRESSURES][N_P_AVERAGE], temp[N_T_AVERAGE], scale; static int first = 1, i_pressure, i_temp; if (first) { +// scale = (XMAX - XMIN)/4.0; + scale = (YMAX - YMIN)/2.5; if (left) { - xbox = XMIN + 0.4; - xtext = XMIN + 0.08; - xxbox = XMAX - 0.39; - xxtext = XMAX - 0.73; + xbox = XMIN + 0.45*scale; + xtext = XMIN + 0.12*scale; + xxbox = XMAX - 0.39*scale; + xxtext = XMAX - 0.73*scale; } else { - xbox = XMAX - 0.41; - xtext = XMAX - 0.73; - xxbox = XMIN + 0.4; - xxtext = XMIN + 0.08; + xbox = XMAX - 0.41*scale; + xtext = XMAX - 0.73*scale; + xxbox = XMIN + 0.4*scale; + xxtext = XMIN + 0.08*scale; } - xmid = 0.5*(XMIN + XMAX) - 0.1; - xmidtext = xmid - 0.24; + xmid = 0.5*(XMIN + XMAX) - 0.1*scale; + xmidtext = xmid - 0.24*scale; for (i=0; i 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio; - erase_area_hsl_turbo(xbox, y + 0.025, 0.37, 0.05, hue, 0.9, 0.5); -// erase_area_hsl_turbo(xmid + 0.1, y + 0.025, 0.4, 0.05, hue, 0.9, 0.5); + if (PRINT_LEFT) erase_area_hsl_turbo(xbox, y + 0.025*scale, 0.5*scale, 0.05*scale, hue, 0.9, 0.5); + else erase_area_hsl_turbo(xmid + 0.1, y + 0.025*scale, 0.45*scale, 0.05*scale, hue, 0.9, 0.5); if ((hue < 90)||(hue > 270)) glColor3f(1.0, 1.0, 1.0); else glColor3f(0.0, 0.0, 0.0); sprintf(message, "Temperature %.2f", 1.0/beta); - write_text(xtext, y, message); -// write_text(xmidtext, y, message); + if (PRINT_LEFT) write_text(xtext, y, message); + else write_text(xmidtext, y, message); // y -= 0.1; // erase_area_hsl(xxbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); @@ -3369,17 +3617,17 @@ void print_parameters(double beta, double temperature, double krepel, double len if (DECREASE_CONTAINER_SIZE) /* print density */ { density = (double)ncircles/((lengthcontainer)*(INITYMAX - INITYMIN)); - erase_area_hsl(xbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); + erase_area_hsl(xbox, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Density %.3f", density); write_text(xtext, y, message); - erase_area_hsl(xmid, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); + erase_area_hsl(xmid, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Temperature %.2f", temperature); write_text(xmidtext, y, message); - erase_area_hsl(xxbox, y + 0.025, 0.37, 0.05, 0.0, 0.9, 0.0); + erase_area_hsl(xxbox, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Pressure %.3f", meanpressure); write_text(xxtext, y, message); @@ -3387,7 +3635,7 @@ void print_parameters(double beta, double temperature, double krepel, double len } else if (INCREASE_KREPEL) /* print force constant */ { - erase_area_hsl(xbox, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + erase_area_hsl(xbox, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Force %.0f", krepel); write_text(xtext + 0.28, y, message); @@ -3436,7 +3684,7 @@ void print_parameters(double beta, double temperature, double krepel, double len hue = PARTICLE_HUE_MIN + 0.5*(PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*mean_temp/PARTICLE_EMAX; if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX; - erase_area_hsl_turbo(xbox, y + 0.025, 0.37, 0.05, hue, 0.9, 0.5); + erase_area_hsl_turbo(xbox, y + 0.025*scale, 0.37*scale, 0.05*scale, hue, 0.9, 0.5); if ((hue < 90)||(hue > 270)) glColor3f(1.0, 1.0, 1.0); else glColor3f(0.0, 0.0, 0.0); sprintf(message, "Temperature %.2f", mean_temp); @@ -3445,7 +3693,7 @@ void print_parameters(double beta, double temperature, double krepel, double len if (INCREASE_GRAVITY) { - erase_area_hsl(xmid, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + erase_area_hsl(xmid, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Gravity %.2f", gravity/GRAVITY); write_text(xmidtext + 0.1, y, message); @@ -3649,40 +3897,42 @@ void print_omega(double angular_speed) void print_segments_speeds(double vx[2], double vy[2]) { char message[100]; - double rgb[3], y = YMAX - 0.1; - static double xleftbox, xlefttext, xrightbox, xrighttext, ymin = YMIN + 0.05; + double rgb[3], y; + static double xleftbox, xlefttext, xrightbox, xrighttext, ymin = YMIN + 0.05, scale; static int first = 1; if (first) { - xleftbox = XMIN + 0.4; - xlefttext = xleftbox - 0.3; - xrightbox = XMAX - 0.39; - xrighttext = xrightbox - 0.3; + scale = (XMAX - XMIN)/4.0; + xleftbox = XMIN + 0.4*scale; + xlefttext = xleftbox - 0.3*scale; + xrightbox = XMAX - 0.39*scale; + xrighttext = xrightbox - 0.3*scale; first = 0; } - erase_area_hsl(xrightbox, y + 0.025, 0.25, 0.05, 0.0, 0.9, 0.0); + y = YMAX - 0.1*scale; + erase_area_hsl(xrightbox, y + 0.025*scale, 0.25*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Vx = %.2f", vx[0]); write_text(xrighttext + 0.1, y, message); - y -= 0.1; - erase_area_hsl(xrightbox, y + 0.025, 0.25, 0.05, 0.0, 0.9, 0.0); + y -= 0.1*scale; + erase_area_hsl(xrightbox, y + 0.025*scale, 0.25*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Vy = %.2f", vy[0]); write_text(xrighttext + 0.1, y, message); if (TWO_OBSTACLES) { - y = YMAX - 0.1; - erase_area_hsl(xleftbox, y + 0.025, 0.25, 0.05, 0.0, 0.9, 0.0); + y = YMAX - 0.1*scale; + erase_area_hsl(xleftbox, y + 0.025*scale, 0.25*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Vx = %.2f", vx[1]); write_text(xlefttext + 0.1, y, message); - y -= 0.1; - erase_area_hsl(xleftbox, y + 0.025, 0.25, 0.05, 0.0, 0.9, 0.0); + y -= 0.1*scale; + erase_area_hsl(xleftbox, y + 0.025*scale, 0.25*scale, 0.05*scale, 0.0, 0.9, 0.0); glColor3f(1.0, 1.0, 1.0); sprintf(message, "Vy = %.2f", vy[1]); write_text(xlefttext + 0.1, y, message); @@ -4182,24 +4432,33 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl x = particle[j].xc; y = particle[j].yc; - if ((x > XMAX + padding)||(x < XMIN - padding)||(y > YMAX + padding)||(y < YMIN - padding)) + if ((x > BCXMAX + padding)||(x < BCXMIN - padding)||(y > BCYMAX + padding)||(y < BCYMIN - padding)) { particle[j].active = 0; particle[j].vx = 0.0; particle[j].vy = 0.0; - particle[j].xc = XMAX + 2.0*padding; - particle[j].yc = YMAX + 2.0*padding; + particle[j].xc = BCXMAX + 2.0*padding; + particle[j].yc = BCYMAX + 2.0*padding; } -// if ((x > XMAX)&&(x < XMAX + 0.5*padding)) particle[j].fx += KSPRING_BOUNDARY*(XMAX + 0.5*padding - x); -// else if (x < XMIN) particle[j].fx -= KSPRING_BOUNDARY; -// if (y > YMAX) particle[j].fy += KSPRING_BOUNDARY; -// else if (y < YMIN) particle[j].fy -= KSPRING_BOUNDARY; + return(fperp); + } + case (BC_REFLECT_ABS): + { + /* add harmonic force outside screen */ + padding = 0.1; + x = particle[j].xc; + y = particle[j].yc; -// if (x > XMAX + padding) particle[j].fx -= KSPRING_BOUNDARY*(x - XMAX - padding); -// else if (x < XMIN - padding) particle[j].fx += KSPRING_BOUNDARY*(XMIN - padding - x); -// if (y > YMAX + padding) particle[j].fy -= KSPRING_BOUNDARY*(y - YMAX - padding); -// else if (y < YMIN - padding) particle[j].fy += KSPRING_BOUNDARY*(YMIN - padding - y); + if ((x > BCXMAX + padding)||(y > BCYMAX + padding)||(y < BCYMIN - padding)) + { + particle[j].active = 0; + particle[j].vx = 0.0; + particle[j].vy = 0.0; + particle[j].xc = BCXMAX + 2.0*padding; + particle[j].yc = BCYMAX + 2.0*padding; + } + else if (particle[j].yc < BCYMIN) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - particle[j].yc); return(fperp); } @@ -4236,6 +4495,55 @@ void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLE } +int reorder_particles(t_particle particle[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES]) +/* keep only active particles */ +{ + int i, k, new = 0, nactive = 0; + + 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 xmin)&&(y > ymin)); -// return((x + 0.5*(double)size > xmin)&&(y + 0.5*(double)size > ymin)); +} + +int in_plot_box_screencoord(double x, double y) +{ + static double xmin, ymin; + static int first = 1; + + if (first) + { + xmin = XMAX - 1.0; + ymin = YMAX - 1.0; + first = 0; + } + + return((x > xmin)&&(y > ymin)); } double size_ratio_color(int clustersize, int ncells) @@ -916,6 +983,11 @@ void draw_cell(int c, int nx, int ny, int size) glEnd(); break; } + case (PLOT_POISSON_DISC): + { + /* beta version, TO DO */ +// draw_circle(); + } } } @@ -939,11 +1011,11 @@ void test_neighbours(int start, t_perco *cell, int nx, int ny, int size, int nce -void draw_configuration(t_perco *cell, int *cluster_sizes, int nx, int ny, int size, int max_cluster_size) +void draw_configuration(t_perco *cell, int *cluster_sizes, int ncells, int nx, int ny, int size, int max_cluster_size) /* draw the configuration */ { - int i, j, ishift, k; - double rgb[3], x, y, alpha, r, fade = 0, dsize; + int i, j, ishift, k, n; + double rgb[3], x, y, alpha, r, fade = 0, dsize, radius, x2, y2; static double h, h1; static int first = 1; @@ -1142,7 +1214,7 @@ void draw_configuration(t_perco *cell, int *cluster_sizes, int nx, int ny, int s x = 0.5*((double)i + 1.25)*dsize; y = h*dsize*((double)j); - if ((ADD_PLOT)&&(in_plot_box(x, y))) fade = 1; + if ((ADD_PLOT)&&(in_plot_box(x, y+0.5*h*dsize))) fade = 1; else fade = 0; set_cell_color(cell[2*j*nx+i], cluster_sizes, fade, max_cluster_size); @@ -1164,6 +1236,41 @@ void draw_configuration(t_perco *cell, int *cluster_sizes, int nx, int ny, int s } break; } + case (PLOT_POISSON_DISC): + { + radius = sqrt((XMAX - XMIN)*(YMAX - YMIN)/(PI*(double)nx));; + + blank(); + + if (ADD_PLOT) + { + rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0; + erase_area_rgb(XMAX - 0.5, YMAX - 0.5, 0.5, 0.5, rgb); + } + + for (i=0; i 0)&&(ncircles < nmaxcells)) + { + /* randomly select an active circle */ + i = rand()%(ncircles); + while (!active_poisson[i]) i = rand()%(ncircles); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= dpoisson*dpoisson); + /* new circle is in domain */ + far = far*(x < XMAX)*(x > XMIN)*(y < YMAX)*(y > YMIN); + } + if (far) /* accept new circle */ + { + printf("New circle at (%.3f,%.3f) accepted\n", x, y); + cell[ncircles].x = x; + cell[ncircles].y = y; + cell[ncircles].active = 1; + active_poisson[ncircles] = 1; + ncircles++; + n_p_active++; + naccepted++; + } +// else printf("Rejected\n"); + } + if (naccepted == 0) /* inactivate circle i */ + { + active_poisson[i] = 0; + n_p_active--; + } + printf("%i active circles\n", n_p_active); + } + + printf("Generated %i circles\n", ncircles); + + free(active_poisson); + + return(ncircles); +} + + int cell_number(int nx, int ny) /* total number of cells in graph */ { @@ -1270,6 +1445,7 @@ int cell_number(int nx, int ny) case (BC_HEX_BOND_DIRICHLET): return(3*(int)((double)((nx+2)*(ny+2))*2.0/sqrt(3.0))); /* hex lattice requires more vertical space! */ case (BC_TRIANGLE_SITE_DIRICHLET): return((int)((double)(2*nx*ny)*2.0/sqrt(3.0))); /* triangle lattice requires more vertical space! */ + case (BC_POISSON_DISC): return(nx*ny); /* TO IMPROVE */ } } @@ -1296,6 +1472,13 @@ void compute_nxny(int size, int *nx, int *ny) *ny = (int)((double)NY*2.0/((double)size*sqrt(3.0))) + 1; break; } + case (BC_POISSON_DISC): + { + /* for Poisson disc configuration, use a 1d labelling */ + *nx = NX*NY/(size*size); + *ny = 1; + break; + } default: { *nx = NX/size; @@ -1308,7 +1491,9 @@ void compute_nxny(int size, int *nx, int *ny) int init_cell_lattice(t_perco *cell, int nx, int ny) /* initialize the graph of connected cells - returns the number of cells */ { - int i, j, iplus, iminus, ishift, n; + int i, j, k, iplus, iminus, ishift, n; + int ncells; + double dpoisson, radius; printf("Initializing cell lattice ..."); @@ -2073,6 +2258,33 @@ int init_cell_lattice(t_perco *cell, int nx, int ny) return(2*nx*ny); break; } + case (BC_POISSON_DISC): + { + dpoisson = 3.0*sqrt((XMAX - XMIN)*(YMAX - YMIN)/(PI*(double)nx)); + +// printf("nx = %i, dpoisson = %.5lg\n", nx, dpoisson); + ncells = generate_poisson_discs(cell, dpoisson, nx); + radius = 1.8*dpoisson; + + for (i=0; i YMIN + LAMBDA); else return (y > YMIN + LAMBDA + s); } + case (D_FABRY_PEROT): + { + return(vabs(x - y*LAMBDA/YMAX) > 0.5*MU); + } + case (D_LSHAPE): + { + if (vabs(x) > LAMBDA) return(0); + else if (vabs(y) > 1.0) return(0); + else if ((x > 0.0)&&(y > 0.0)) return(0); + else return(1); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -2335,6 +2360,17 @@ int xy_in_billiard(double x, double y) } } +int xy_in_billiard(double x, double y) +/* returns 1 if (x,y) represents a point in the billiard */ +{ + if (COMPARISON) + { + if (y > 0.0) return (xy_in_billiard_single_domain(x, y, B_DOMAIN, ncircles, circles)); + else return (xy_in_billiard_single_domain(x, y, B_DOMAIN_B, ncircles_b, circles_b)); + } + else return (xy_in_billiard_single_domain(x, y, B_DOMAIN, ncircles, circles)); +} + int ij_in_billiard(int i, int j) /* returns 1 if (i,j) represents a point in the billiard */ { @@ -3166,6 +3202,13 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound glEnd(); break; } + case (D_FABRY_PEROT): + { + glLineWidth(BOUNDARY_WIDTH); + draw_line(-LAMBDA - 0.5*MU, YMIN, LAMBDA - 0.5*MU, YMAX); + draw_line(-LAMBDA + 0.5*MU, YMIN, LAMBDA + 0.5*MU, YMAX); + break; + } case (D_CIRCLE_SEGMENT): { glLineWidth(BOUNDARY_WIDTH); @@ -3565,7 +3608,7 @@ void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int p xy_to_ij(x2, y2, ij_topright); rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; - erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb); +// erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb); if (ROTATE_COLOR_SCHEME) { @@ -3671,7 +3714,7 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, xy_to_ij(x2, y2, ij_topright); rgb[0] = 0.0; rgb[1] = 0.0; rgb[2] = 0.0; - erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb); +// erase_area_rgb(0.5*(x1 + x2), x2 - x1, 0.5*(y1 + y2), y2 - y1, rgb); if (ROTATE_COLOR_SCHEME) { @@ -3770,3 +3813,474 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, draw_rectangle(x1, y1, x2, y2); } + +void print_speed(double speed, int fade, double fade_value) +{ + char message[100]; + double y = YMAX - 0.1, pos[2]; + static double xleftbox, xlefttext; + static int first = 1; + + if (first) + { + xleftbox = XMIN + 0.3; + xlefttext = xleftbox - 0.45; + first = 0; + } + + erase_area_hsl(xleftbox, y + 0.025, 0.22, 0.05, 0.0, 0.9, 0.0); + if (fade) glColor3f(fade_value, fade_value, fade_value); + else glColor3f(1.0, 1.0, 1.0); + xy_to_pos(xlefttext + 0.28, y, pos); + sprintf(message, "Mach %.3lg", speed); + write_text(pos[0], pos[1], message); +} + + +void init_laplacian_coords(t_laplacian laplace[NX*NY], double phi[NX*NY]) +/* compute coordinates of neighbours to compute Laplacian */ +{ + int i, j, iplus, iminus, i1, i2, i3, j1, j2, j3, ij[2];; + + printf("Initialising Laplacian table\n"); + + /* Laplacian in the bulk */ + #pragma omp parallel for private(i,j) + for (i=1; i margin*length*olength)&&(x2*observer[0] + y2*observer[1] > margin*length1*olength)) + draw_segment_hsl(x1, y1, x2, y2, h, s, l); +} + + +void draw_segment_rgb_visible(double x1, double y1, double x2, double y2, double r, double g, double b, double margin) +/* hack to draw the billiard boundary in front of the wave */ +/* only parts of the boundary having a small enough angle with the observer vector are drawn */ +{ + double length, length1, olength; + + olength = module2(observer[0], observer[1]); + + length = module2(x1,y1); + length1 = module2(x2,y2); + if ((x1*observer[0] + y1*observer[1] > margin*length*olength)&&(x2*observer[0] + y2*observer[1] > margin*length1*olength)) + draw_segment_rgb(x1, y1, x2, y2, r, g, b); +} + void draw_billiard_3d_front(int fade, double fade_value) /* hack to draw the billiard boundary in front of the wave */ /* only parts of the boundary having a small enough angle with the observer vector are drawn */ { - 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, length, length1; + 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, length, length1, padd, margin; int i, j, k, k1, k2, mr2; static int first = 1, nsides; static double olength; @@ -805,6 +875,20 @@ void draw_billiard_3d_front(int fade, double fade_value) glEnd(); break; } + case (D_ELLIPSE): + { + glBegin(GL_LINE_LOOP); + dphi = DPI/(double)NSEG; + margin = 0.0; + for (i=0; i<=NSEG; i++) + { + phi = (double)i*dphi; + draw_segment_rgb_visible(LAMBDA*cos(phi), sin(phi), LAMBDA*cos(phi+dphi), sin(phi+dphi), fade_value, fade_value, fade_value, margin); + draw_vertex_x_y_z(LAMBDA*cos(phi), sin(phi), 0.0); + } + glEnd (); + break; + } case (D_YOUNG): { glBegin(GL_LINE_STRIP); @@ -852,6 +936,21 @@ void draw_billiard_3d_front(int fade, double fade_value) draw_polyline_visible(i%npolyline, (i+1)%npolyline, -0.5); break; } + case (D_LSHAPE): + { + padd = 0.005; + margin = 0.0; + glLineWidth(BOUNDARY_WIDTH); + draw_segment_hsl_visible(-LAMBDA - padd, -1.0 - padd, 0.0, -1.0 - padd, 0.0, 1.0, 0.5*fade_value, margin); + draw_segment_hsl_visible(-LAMBDA - padd, 1.0 + padd, 0.0, 1.0 + padd, 0.0, 1.0, 0.5*fade_value, margin); + draw_segment_hsl_visible(0.0, -1.0 - padd, LAMBDA + padd, -1.0 - padd, 220.0, 1.0, 0.5*fade_value, margin); + draw_segment_hsl_visible(0.0, padd, LAMBDA + padd, padd, 220.0, 1.0, 0.5*fade_value, 0.2); + draw_segment_hsl_visible(-LAMBDA - padd, -1.0 - padd, -LAMBDA - padd, 0.0, 60.0, 1.0, 0.5*fade_value, margin); + draw_segment_hsl_visible( LAMBDA + padd, -1.0 - padd, LAMBDA + padd, 0.0, 60.0, 1.0, 0.5*fade_value, margin); + draw_segment_hsl_visible(-LAMBDA - padd, 0.0, -LAMBDA - padd, 1.0 + padd, 180.0, 1.0, 0.5*fade_value, margin); + draw_segment_hsl_visible( padd, padd, padd, LAMBDA + padd, 180.0, 1.0, 0.5*fade_value, 0.6); + break; + } default: { break; @@ -1118,8 +1217,8 @@ void compute_wave_fields(double phi[NX*NY], double psi[NX*NY], short int xy_in[N if (COMPUTE_ENERGY) compute_energy_field(phi, psi, xy_in, wave); -// if ((zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_ENERGY)) -// compute_log_energy_field(phi, psi, xy_in, wave); + if ((zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_ENERGY)) + compute_log_energy_field(phi, psi, xy_in, wave); if ((zplot == P_3D_PHASE)||(cplot == P_3D_PHASE)) compute_phase_field(phi, psi, xy_in, wave); @@ -1132,7 +1231,7 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc { int i, j, k; double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB; - double u, v, u1, x, y, xy[2], norm2, speed; + double u, v, u1, x, y, xy[2], norm2, speed, r2, c; if (VARIABLE_IOR) { @@ -1211,6 +1310,23 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc } break; } + case (IOR_EARTH): + { + for (i=0; i 1.0) c = 0.0; + else if (r2 < 0.25*0.25) c = 0.8*COURANT; + else if (r2 < 0.58*0.58) c = COURANT*(0.68 - 0.55*r2); + else c = COURANT*(1.3 - 0.9*r2); + tc[i*NY+j] = c; + tcc[i*NY+j] = c; + tgamma[i*NY+j] = GAMMA; + } + } + break; + } default: { for (i=0; i 0.0)&&(observer_angle < PID)) { - for (j=0; j0; i--) - for (j=0; j0; j--) - for (i=0; i0)) + { + iplus = i+1; + iminus = i-1; if (iminus < 0) iminus = 0; + /* TO ADAPT */ + phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta[i*NY+NY-1] - KAPPA*x - tgamma[i*NY+NY-1]*(x-y); + } + + else switch (B_COND) { + case (BC_DIRICHLET): + { + phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta[i*NY+NY-1] - KAPPA*x - tgamma[i*NY+NY-1]*(x-y); + break; + } + case (BC_PERIODIC): + { + phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta[i*NY+NY-1] - KAPPA*x - tgamma[i*NY+NY-1]*(x-y); + break; + } + case (BC_ABSORBING): + { + phi_out[i*NY+NY-1] = x - tc[i*NY+NY-1]*(x - phi_in[i*NY+NY-2]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y); + break; + } + case (BC_VPER_HABS): + { + if (i==0) phi_out[NY-1] = x - tc[NY-1]*(x - phi_in[1*NY+NY-1]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y); + else phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta[i*NY+NY-1] - KAPPA*x - tgamma[i*NY+NY-1]*(x-y); + break; + } + } + } + } + + /* bottom boundary */ + for (i=0; i0)) + { + /* TO ADAPT */ + phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta[i*NY] - KAPPA*x - tgamma[i*NY]*(x-y); + } + + else switch (B_COND) { + case (BC_DIRICHLET): + { + phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta[i*NY] - KAPPA*x - tgamma[i*NY]*(x-y); + break; + } + case (BC_PERIODIC): + { + phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta[i*NY] - KAPPA*x - tgamma[i*NY]*(x-y); + break; + } + case (BC_ABSORBING): + { + phi_out[i*NY] = x - tc[i*NY]*(x - phi_in[i*NY+1]) - KAPPA_TOPBOT*x - GAMMA_TOPBOT*(x-y); + break; + } + case (BC_VPER_HABS): + { + if (i==0) phi_out[0] = x - tc[0]*(x - phi_in[NY]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y); + else phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta[i*NY] - KAPPA*x - tgamma[i*NY]*(x-y); + break; + } + } + } + } + + /* add oscillating boundary condition on the left corners */ + if (OSCILLATE_LEFT) + { + phi_out[0] = AMPLITUDE*cos((double)time*OMEGA); + phi_out[NY-1] = AMPLITUDE*cos((double)time*OMEGA); + } + + /* for debugging purposes/if there is a risk of blow-up */ + if (FLOOR) + { + #pragma omp parallel for private(i,j) + for (i=0; i VMAX) phi_out[i*NY+j] = VMAX; + if (phi_out[i*NY+j] < -VMAX) phi_out[i*NY+j] = -VMAX; + } + } + } + } + + free(delta); +} + + +void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out[NX*NY], + short int xy_in[NX*NY], double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY]) +/* 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, jtop, jbot; + double delta, x, y, c, cc, gamma; + static long time = 0; + static short int first = 1; + static double tb_shift; + + time++; + + if ((OSCILLATE_TOPBOT)&&(first)) + { + tb_shift = (int)((XMAX - XMIN)*(double)NX/(XMAX - XMIN)); + if (tb_shift >= NX-1) tb_shift = NY - 2; + first = 0; + } + #pragma omp parallel for private(i,j,iplus,iminus,jplus,jminus,delta,x,y) /* evolution in the bulk */ for (i=1; i0)) + { + iplus = i+1; + iminus = i-1; if (iminus < 0) iminus = 0; + delta = phi_in[iplus*NY+NY-1] + phi_in[iminus*NY+NY-1] + - 2.0*x; + phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta - KAPPA*x - tgamma[i*NY+NY-1]*(x-y); + } + + else switch (B_COND) { case (BC_DIRICHLET): { iplus = i+1; if (iplus == NX) iplus = NX-1; @@ -392,7 +639,7 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out iminus = (i-1) % NX; if (iminus < 0) iminus += NX; - delta = phi_in[iplus*NY+NY-1] + phi_in[iminus*NY+NY-1] + phi_in[i*NY+NY-2] + phi_in[i*NY] - 4.0*x; + delta = phi_in[iplus*NY+NY-1] + phi_in[iminus*NY+NY-1] + phi_in[i*NY+NY-2] + phi_in[i*NY+jbot] - 4.0*x; phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta - KAPPA*x - tgamma[i*NY+NY-1]*(x-y); break; } @@ -410,7 +657,7 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out iplus = (i+1); if (iplus == NX) iplus = NX-1; iminus = (i-1); if (iminus == -1) iminus = 0; - delta = phi_in[iplus*NY+NY-1] + phi_in[iminus*NY+NY-1] + phi_in[i*NY+NY-2] + phi_in[i*NY] - 4.0*x; + delta = phi_in[iplus*NY+NY-1] + phi_in[iminus*NY+NY-1] + phi_in[i*NY+NY-2] + phi_in[i*NY+jbot] - 4.0*x; if (i==0) phi_out[NY-1] = x - tc[NY-1]*(x - phi_in[1*NY+NY-1]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y); else phi_out[i*NY+NY-1] = -y + 2*x + tcc[i*NY+NY-1]*delta - KAPPA*x - tgamma[i*NY+NY-1]*(x-y); break; @@ -420,12 +667,22 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out } /* bottom boundary */ + if (COMPARISON) jtop = NY/2-1; + else jtop = NY-1; for (i=0; i0)) + { + iplus = i+1; + iminus = i-1; if (iminus < 0) iminus = 0; + delta = phi_in[iplus*NY] + phi_in[iminus*NY] + - 2.0*x; + phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta - KAPPA*x - tgamma[i*NY]*(x-y); + } + + else switch (B_COND) { case (BC_DIRICHLET): { iplus = i+1; if (iplus == NX) iplus = NX-1; @@ -441,7 +698,7 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out iminus = (i-1) % NX; if (iminus < 0) iminus += NX; - delta = phi_in[iplus*NY] + phi_in[iminus*NY] + phi_in[i*NY+1] + phi_in[i*NY+NY-1] - 4.0*x; + delta = phi_in[iplus*NY] + phi_in[iminus*NY] + phi_in[i*NY+1] + phi_in[i*NY+jtop] - 4.0*x; phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta - KAPPA*x - tgamma[i*NY]*(x-y); break; } @@ -459,7 +716,7 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out iplus = (i+1); if (iplus == NX) iplus = NX-1; iminus = (i-1); if (iminus == -1) iminus = 0; - delta = phi_in[iplus*NY] + phi_in[iminus*NY] + phi_in[i*NY+1] + phi_in[i*NY+NY-1] - 4.0*x; + delta = phi_in[iplus*NY] + phi_in[iminus*NY] + phi_in[i*NY+1] + phi_in[i*NY+jtop] - 4.0*x; if (i==0) phi_out[0] = x - tc[0]*(x - phi_in[NY]) - KAPPA_SIDES*x - GAMMA_SIDES*(x-y); else phi_out[i*NY] = -y + 2*x + tcc[i*NY]*delta - KAPPA*x - tgamma[i*NY]*(x-y); break; @@ -468,13 +725,113 @@ void evolve_wave_half(double phi_in[NX*NY], double psi_in[NX*NY], double phi_out } } + /* case of comparisons - top of bottom half */ + if (COMPARISON) for (i=0; i0)) + { + iplus = i+1; + iminus = i-1; if (iminus < 0) iminus = 0; + delta = phi_in[iplus][NY-1] + phi_in[iminus][NY-1] + - 2.0*x; + phi_out[i][NY-1] = -y + 2*x + tcc[i][NY-1]*delta - KAPPA*x - tgamma[i][NY-1]*(x-y); + } + + else switch (B_COND) { case (BC_DIRICHLET): { iplus = i+1; if (iplus == NX) iplus = NX-1; @@ -398,6 +415,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX break; } } +// psi_out[i][NY-1] = x; } } @@ -407,7 +425,15 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX x = phi_in[i][0]; y = psi_in[i][0]; - switch (B_COND) { + if ((OSCILLATE_TOPBOT)&&(i < tb_shift)&&(i0)) + { + iplus = i+1; + iminus = i-1; if (iminus < 0) iminus = 0; + delta = phi_in[iplus][0] + phi_in[iminus][0] + - 2.0*x; + phi_out[i][0] = -y + 2*x + tcc[i][0]*delta - KAPPA*x - tgamma[i][0]*(x-y); + } + + else switch (B_COND) { case (BC_DIRICHLET): { iplus = i+1; if (iplus == NX) iplus = NX-1; @@ -447,14 +473,15 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX break; } } +// psi_out[i][0] = x; } } /* add oscillating boundary condition on the left corners */ if (OSCILLATE_LEFT) { - phi_out[0][0] = AMPLITUDE*cos((double)time*OMEGA)*exp(-(double)time*DAMPING); - phi_out[0][NY-1] = AMPLITUDE*cos((double)time*OMEGA)*exp(-(double)time*DAMPING); + phi_out[0][0] = oscillating_bc(time); + phi_out[0][NY-1] = oscillating_bc(time); } /* for debugging purposes/if there is a risk of blow-up */ @@ -464,6 +491,8 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX { if (phi_out[i][j] > VMAX) phi_out[i][j] = VMAX; if (phi_out[i][j] < -VMAX) phi_out[i][j] = -VMAX; +// if (psi_out[i][j] > VMAX) psi_out[i][j] = VMAX; +// if (psi_out[i][j] < -VMAX) psi_out[i][j] = -VMAX; } } } @@ -474,7 +503,7 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *x /* time step of field evolution */ /* phi is value of field at time t, psi at time t-1 */ { - // For the purpose of these comments w[t], w[t-1], w[t+1] are used to refer + // For the purpose of these comments w[t], w[t-1], w[t+1] are used to refer // to phi, psi and the result respectively to avoid confusion with the // passed parameter names. // At the beginning w[t] is saved in phi, w[t-1] in psi and tmp is space @@ -518,7 +547,8 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl void animation() { - double time, scale, ratio, startleft[2], startright[2], sign, r2, xy[2], fade_value; + double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c; +// double *phi[NX], *psi[NX], *phi_tmp[NX], *psi_tmp[NX], *total_energy[NX], *color_scale[NX]; double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX]; short int *xy_in[NX]; int i, j, s, sample_left[2], sample_right[2], period = 0, fade; @@ -536,6 +566,8 @@ void animation() { phi[i] = (double *)malloc(NY*sizeof(double)); psi[i] = (double *)malloc(NY*sizeof(double)); +// phi_tmp[i] = (double *)malloc(NY*sizeof(double)); +// psi_tmp[i] = (double *)malloc(NY*sizeof(double)); tmp[i] = (double *)malloc(NY*sizeof(double)); total_energy[i] = (double *)malloc(NY*sizeof(double)); xy_in[i] = (short int *)malloc(NY*sizeof(short int)); @@ -543,8 +575,8 @@ void animation() } /* initialise positions and radii of circles */ - 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); + if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) ncircles = init_circle_config(circles); + else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config(polygons); printf("Polygons initialized\n"); /* initialise polyline for von Koch and similar domains */ @@ -553,6 +585,7 @@ void animation() courant2 = COURANT*COURANT; courantb2 = COURANTB*COURANTB; + c = COURANT*(XMAX - XMIN)/(double)NX; /* initialize color scale, for option RESCALE_COLOR_IN_CENTER */ if (RESCALE_COLOR_IN_CENTER) @@ -627,13 +660,23 @@ void animation() blank(); glColor3f(0.0, 0.0, 0.0); // draw_wave(phi, psi, xy_in, 1.0, 0, PLOT); - if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, 1.0, 0, PLOT, COLOR_PALETTE); + if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0); else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0); draw_billiard(0, 1.0); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value); + if (PRINT_SPEED) + { + a = 0.0075; + b = 0.00015; +// speed = a/((double)(NVID)*c); +// speed = 0.55*a/((double)(NVID*OSCILLATING_SOURCE_PERIOD)*c); + speed = a/((double)(3*NVID*OSCILLATING_SOURCE_PERIOD)*c); + /* the factor 3 is due to evolve_wave calling evolve_wave_half 3 times */ + print_speed(speed, 0, 1.0); + } glutSwapBuffers(); @@ -654,10 +697,11 @@ void animation() else scale = 1.0; // draw_wave(phi, psi, xy_in, scale, i, PLOT); - if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, i, PLOT, COLOR_PALETTE); + if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0); else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0); for (j=0; j= INITIAL_TIME)&&(DOUBLE_MOVIE)) { // draw_wave(phi, psi, xy_in, scale, i, PLOT_B); - if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B); + if (HIGHRES) + draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0); else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0); draw_billiard(0, 1.0); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); + if (PRINT_SPEED) print_speed(speed, 0, 1.0); glutSwapBuffers(); save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); // save_frame_counter(NSTEPS + 21 + counter); @@ -720,42 +777,50 @@ void animation() if (DOUBLE_MOVIE) { // draw_wave(phi, psi, xy_in, scale, i, PLOT); - if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE); + if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0); else draw_wave_epalette(phi, psi, total_energy, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0); draw_billiard(0, 1.0); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); + if (PRINT_SPEED) print_speed(speed, 0, 1.0); glutSwapBuffers(); } if (!FADE) for (i=0; ii3-1)) iplus = i1; +// else if ((j>=j2)&&(iplus>i2-1)) iplus = i1; +// +// iminus = i-1; +// if ((j=j2)&&(iminusj3-1)) jplus = j1; +// else if ((i>=i2)&&(jplus>j2-1)) jplus = j1; +// +// jminus = j-1; +// if ((i=i2)&&(jminus=i2-1)&&(j>=j2-1)) return(0.0); + if ((i>=i3-1)||(j>=j3-1)) return(0.0); + } + + gradientx2 = (phi[iplus*NY+j]-phi[i*NY+j])*(phi[iplus*NY+j]-phi[i*NY+j]) + (phi[i*NY+j] - phi[iminus*NY+j])*(phi[i*NY+j] - phi[iminus*NY+j]); gradienty2 = (phi[i*NY+jplus]-phi[i*NY+j])*(phi[i*NY+jplus]-phi[i*NY+j]) diff --git a/wave_comparison.c b/wave_comparison.c index d3f627a..37556a1 100644 --- a/wave_comparison.c +++ b/wave_comparison.c @@ -219,6 +219,12 @@ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 5.0 /* max value of wave amplitude */ +/* the following constants are only used by wave_billiard and wave_3d so far */ +#define COMPARISON 0 /* set to 1 to compare two different patterns */ +#define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */ +#define ACHIRP 0.2 /* acceleration coefficient in chirp */ +#define DAMPING 0.0 /* damping of periodic excitation */ +/* end of constants only used by wave_billiard and wave_3d */ #include "global_pdes.c" /* constants and global variables */ #include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ @@ -863,5 +869,4 @@ int main(int argc, char** argv) glutMainLoop(); return 0; -} - +} \ No newline at end of file diff --git a/wave_energy.c b/wave_energy.c index c23868a..4a82454 100644 --- a/wave_energy.c +++ b/wave_energy.c @@ -196,6 +196,12 @@ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 5.0 /* max value of wave amplitude */ +/* the following constants are only used by wave_billiard and wave_3d so far */ +#define COMPARISON 0 /* set to 1 to compare two different patterns */ +#define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */ +#define ACHIRP 0.2 /* acceleration coefficient in chirp */ +#define DAMPING 0.0 /* damping of periodic excitation */ +/* end of constants only used by wave_billiard and wave_3d */ #include "global_pdes.c" /* constants and global variables */ #include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ @@ -917,5 +923,4 @@ int main(int argc, char** argv) glutMainLoop(); return 0; -} - +} \ No newline at end of file