From cb2bcdb13ddf17255c5c005663d3113662b32883 Mon Sep 17 00:00:00 2001 From: Nils Berglund <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Sun, 20 Nov 2022 23:17:39 +0100 Subject: [PATCH] Add files via upload --- global_ljones.c | 2 + global_particles.c | 1 + global_pdes.c | 8 + heat.c | 1 + lennardjones.c | 115 +++++++-------- mangrove.c | 2 + particle_billiard.c | 59 +++++--- rde.c | 158 ++++++++++++++++---- schrodinger.c | 1 + sub_lj.c | 220 ++++++++++++++++++++++----- sub_maze.c | 121 ++++++++++++++- sub_part_billiard.c | 164 +++++++++++++++++++-- sub_perco_3d.c | 185 +++++++++++++++++++++++ sub_rde.c | 206 +++++++++++++++++++++++--- sub_wave.c | 351 +++++++++++++++++++++++++++++++++++++++++++- sub_wave_3d_rde.c | 8 + wave_3d.c | 101 ++++++------- wave_billiard.c | 84 ++++++----- wave_common.c | 122 ++++++++++++++- wave_comparison.c | 2 + wave_energy.c | 2 + 21 files changed, 1625 insertions(+), 288 deletions(-) create mode 100644 sub_perco_3d.c diff --git a/global_ljones.c b/global_ljones.c index 2b0b909..09d71e0 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -64,6 +64,8 @@ #define S_MAZE 15 /* segments forming a maze */ #define S_EXT_RECTANGLE 16 /* particles outside a rectangle */ #define S_DAM_BRICKS 17 /* dam made of several bricks */ +#define S_HLINE_HOLE 18 /* horizontal line with a hole in the bottom */ +#define S_EXT_CIRCLE_RECT 19 /* particles outside a circle and a rectangle */ /* particle interaction */ diff --git a/global_particles.c b/global_particles.c index 16cf1fe..642a416 100644 --- a/global_particles.c +++ b/global_particles.c @@ -79,6 +79,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define P_TREE 7 /* pine tree */ #define P_TOKA_NONSELF 8 /* Tokarsky non-self-unilluminable room */ #define P_MAZE 10 /* maze */ +#define P_MAZE_DIAG 11 /* maze with 45 degrees angles */ /* Color palettes */ diff --git a/global_pdes.c b/global_pdes.c index 0f5d99d..a453062 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -14,6 +14,9 @@ #define D_NOTHING 999 /* no boundaries */ #define D_RECTANGLE 0 /* rectangular domain */ #define D_ELLIPSE 1 /* elliptical domain */ +#define D_EXT_ELLIPSE 199 /* exterior of elliptical domain */ +#define D_EXT_ELLIPSE_CURVED 198 /* exterior of curved elliptical domain */ +#define D_EXT_ELLIPSE_CURVED_BDRY 197 /* exterior of curved elliptical domain, with horizontal boundaries */ #define D_STADIUM 2 /* stadium-shaped domain */ #define D_SINAI 3 /* Sinai billiard */ #define D_DIAMOND 4 /* diamond-shaped billiard */ @@ -65,6 +68,9 @@ #define D_WAVEGUIDE_W 52 /* W-shaped wave guide */ #define D_MAZE 53 /* maze */ #define D_MAZE_CLOSED 54 /* closed maze */ +#define D_CHESSBOARD 55 /* chess board configuration */ +#define D_TRIANGLE_TILES 56 /* triangular tiling */ +#define D_HEX_TILES 57 /* honeycomb tiling */ #define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */ @@ -140,6 +146,8 @@ #define P_MEAN_ENERGY 3 /* energy averaged over time */ #define P_LOG_ENERGY 4 /* log of energy averaged over time */ #define P_LOG_MEAN_ENERGY 5 /* log of energy averaged over time */ +#define P_ENERGY_FLUX 6 /* energy flux */ +#define P_TOTAL_ENERGY_FLUX 7 /* energy flux averaged over time */ /* For Schrodinger equation */ #define P_MODULE 10 /* plot module of wave function squared */ diff --git a/heat.c b/heat.c index e13107e..299cd6c 100644 --- a/heat.c +++ b/heat.c @@ -36,6 +36,7 @@ #include #define MOVIE 0 /* set to 1 to generate movie */ +#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ /* General geometrical parameters */ diff --git a/lennardjones.c b/lennardjones.c index e2d1a36..57b0e49 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -37,42 +37,32 @@ #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 DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ -#define TIME_LAPSE 1 /* set to 1 to add a time-lapse movie at the end */ +#define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */ /* so far incompatible with double movie */ #define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */ -#define TIME_LAPSE_FIRST 1 /* set to 1 to show time-lapse version first */ +#define TIME_LAPSE_FIRST 0 /* set to 1 to show time-lapse version first */ /* General geometrical parameters */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ -// #define XMIN -2.3 -// #define XMAX 3.7 /* x interval */ -// #define YMIN -1.6875 -// #define YMAX 1.6875 /* y interval for 9/16 aspect ratio */ +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.125 +#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define XMIN -3.3 -#define XMAX 4.7 /* x interval */ -#define YMIN -2.25 -#define YMAX 2.25 /* y interval for 9/16 aspect ratio */ +#define INITXMIN -1.95 +#define INITXMAX 1.95 /* x interval for initial condition */ +#define INITYMIN -1.1 +#define INITYMAX 0.15 /* y interval for initial condition */ -#define INITXMIN -2.5 -#define INITXMAX 2.5 /* x interval for initial condition */ -#define INITYMIN -1.7 -#define INITYMAX 0.7 /* y interval for initial condition */ - -// #define BCXMIN -3.1 -// #define BCXMAX 3.1 /* x interval for boundary condition */ -// #define BCYMIN -4.5 -// #define BCYMAX 4.5 /* y interval for boundary condition */ - -#define BCXMIN -5.1 -#define BCXMAX 6.1 /* x interval for boundary condition */ -#define BCYMIN -6.5 -#define BCYMAX 44.0 /* y interval for boundary condition */ +#define BCXMIN -2.0 +#define BCXMAX 2.0 /* x interval for boundary condition */ +#define BCYMIN -1.125 +#define BCYMAX 1.125 /* y interval for boundary condition */ #define OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ @@ -83,7 +73,7 @@ #define OBSTACLE_PATTERN 3 /* pattern of obstacles, see list in global_ljones.c */ #define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */ -#define SEGMENT_PATTERN 102 /* pattern of repelling segments, see list in global_ljones.c */ +#define SEGMENT_PATTERN 19 /* pattern of repelling segments, see list in global_ljones.c */ #define ROCKET_SHAPE 2 /* shape of rocket combustion chamber, see list in global_ljones.c */ #define ROCKET_SHAPE_B 2 /* shape of second rocket */ #define NOZZLE_SHAPE 1 /* shape of nozzle, see list in global_ljones.c */ @@ -103,18 +93,16 @@ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define PDISC_DISTANCE 2.7 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 3.0 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ -// #define MU 0.02 /* parameter controlling radius of particles */ -// #define MU 0.015 /* parameter controlling radius of particles */ +#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */ #define MU 0.009 /* parameter controlling radius of particles */ -// #define MU 0.012 /* parameter controlling radius of particles */ #define MU_B 0.018 /* parameter controlling radius of particles of second type */ -#define NPOLY 25 /* number of sides of polygon */ -#define APOLY 0.666666666 /* angle by which to turn polygon, in units of Pi/2 */ +// #define NPOLY 4 /* number of sides of polygon */ +#define NPOLY 3 /* number of sides of polygon */ +#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ @@ -134,13 +122,11 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 3300 /* number of frames of movie */ -// #define NSTEPS 2000 /* number of frames of movie */ +#define NSTEPS 1500 /* number of frames of movie */ #define NVID 100 /* number of iterations between images displayed on screen */ #define NSEG 250 /* number of segments of boundary */ -#define INITIAL_TIME 10 /* time after which to start saving frames */ -// #define OBSTACLE_INITIAL_TIME 10 /* time after which to start moving obstacle */ -#define OBSTACLE_INITIAL_TIME 200 /* time after which to start moving obstacle */ +#define INITIAL_TIME 250 /* time after which to start saving frames */ +#define OBSTACLE_INITIAL_TIME 350 /* time after which to start moving obstacle */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ #define CONTAINER_WIDTH 4 /* width of container boundary */ @@ -154,7 +140,7 @@ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 20 +#define BOUNDARY_COND 0 /* Plot type, see list in global_ljones.c */ @@ -164,6 +150,7 @@ #define DRAW_BONDS 0 /* set to 1 to draw bonds between neighbours */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ #define FILL_TRIANGLES 1 /* set to 1 to fill triangles between neighbours */ +#define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */ #define COLOR_SEG_GROUPS 1 /* set to 1 to collor segment groups differently */ /* Color schemes */ @@ -191,7 +178,7 @@ #define ENERGY_HUE_MAX 50.0 /* color of saturated particle */ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ -#define PARTICLE_EMAX 2.0e2 /* energy of particle with hottest color */ +#define PARTICLE_EMAX 3.0e2 /* energy of particle with hottest color */ #define HUE_TYPE0 70.0 /* hue of particles of type 0 */ #define HUE_TYPE1 280.0 /* hue of particles of type 1 */ #define HUE_TYPE2 .0 /* hue of particles of type 2 */ @@ -200,10 +187,11 @@ #define RANDOM_RADIUS 0 /* set to 1 for random circle radius */ #define DT_PARTICLE 3.0e-6 /* time step for particle displacement */ #define KREPEL 12.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 4.5 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST 4.0 /* Lennard-Jones equilibrium distance */ #define EQUILIBRIUM_DIST_B 3.5 /* Lennard-Jones equilibrium distance for second type of particle */ #define REPEL_RADIUS 20.0 /* radius in which repelling force acts (in units of particle radius) */ -#define DAMPING 5.0 /* damping coefficient of particles */ +#define DAMPING 10.0 /* damping coefficient of particles */ +#define INITIAL_DAMPING 35.0 /* damping coefficient of particles during initial phase */ #define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ #define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU */ #define PARTICLE_INERTIA_MOMENT 0.2 /* moment of inertia of particle */ @@ -211,7 +199,7 @@ #define V_INITIAL 0.0 /* initial velocity range */ #define OMEGA_INITIAL 10.0 /* initial angular velocity range */ -#define THERMOSTAT 1 /* set to 1 to switch on thermostat */ +#define THERMOSTAT 0 /* set to 1 to switch on thermostat */ #define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ #define SIGMA 5.0 /* noise intensity in thermostat */ #define BETA 0.02 /* initial inverse temperature */ @@ -219,7 +207,7 @@ #define KSPRING_BOUNDARY 1.0e7 /* confining harmonic potential outside simulation region */ #define KSPRING_OBSTACLE 1.0e11 /* harmonic potential of obstacles */ #define NBH_DIST_FACTOR 7.5 /* radius in which to count neighbours */ -#define GRAVITY 15.0 /* gravity acting on all particles */ +#define GRAVITY 1500.0 /* gravity acting on all particles */ #define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */ #define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */ #define GRAVITY_SCHEDULE 2 /* type of gravity schedule, see list in global_ljones.c */ @@ -240,7 +228,7 @@ #define SPIN_RANGE_B 5.0 /* range of spin-spin interaction for second type of particle */ #define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */ -#define INCREASE_BETA 1 /* set to 1 to increase BETA during simulation */ +#define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */ #define BETA_FACTOR 0.025 /* factor by which to change BETA during simulation */ #define N_TOSCILLATIONS 1.5 /* number of temperature oscillations in BETA schedule */ #define NO_OSCILLATION 1 /* set to 1 to have exponential BETA change only */ @@ -301,29 +289,29 @@ #define OMEGAMAX 100.0 /* maximal rotation speed */ #define PRINT_OMEGA 0 /* set to 1 to print angular speed */ #define PRINT_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */ -#define PRINT_SEGMENTS_SPEEDS 1 /* set to 1 to print velocity of moving segments */ +#define PRINT_SEGMENTS_SPEEDS 0 /* set to 1 to print velocity of moving segments */ #define MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */ #define SEGMENTS_MASS 40.0 /* mass of collection of segments */ #define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */ #define SEGMENT_DEACTIVATION_TIME 200 /* time at which to deactivate last segment */ #define RELEASE_ROCKET_AT_DEACTIVATION 1 /* set to 1 to limit segments velocity before segment release */ -#define SEGMENTS_X0 1.5 /* initial position of segments */ -#define SEGMENTS_Y0 0.0 /* initial position of segments */ +#define SEGMENTS_X0 1.0 /* initial position of segments */ +#define SEGMENTS_Y0 0.8 /* initial position of segments */ #define SEGMENTS_VX0 0.0 /* initial velocity of segments */ #define SEGMENTS_VY0 0.0 /* initial velocity of segments */ #define DAMP_SEGS_AT_NEGATIVE_Y 0 /* set to 1 to dampen segments when y coordinate is negative */ #define MOVE_SEGMENT_GROUPS 1 /* set to 1 to group segments into moving units */ -#define SEGMENT_GROUP_MASS 1000.0 /* mass of segment group */ -#define SEGMENT_GROUP_I 1000.0 /* moment of inertia of segment group */ +#define SEGMENT_GROUP_MASS 750.0 /* mass of segment group */ +#define SEGMENT_GROUP_I 500.0 /* moment of inertia of segment group */ #define SEGMENT_GROUP_DAMPING 0.0 /* damping of segment groups */ #define GROUP_REPULSION 1 /* set to 1 for groups of segments to repel each other */ #define KSPRING_GROUPS 1.0e11 /* harmonic potential between segment groups */ #define GROUP_WIDTH 0.05 /* interaction width of groups */ -#define GROUP_G_REPEL 1 /* set to 1 to add repulsion between centers of mass of groups */ +#define GROUP_G_REPEL 0 /* set to 1 to add repulsion between centers of mass of groups */ #define GROUP_G_REPEL_RADIUS 1.2 /* radius within which centers of mass of groups repel each other */ -#define TRACK_SEGMENT_GROUPS 1 /* set to 1 for view to track group of segments */ +#define TRACK_SEGMENT_GROUPS 0 /* set to 1 for view to track group of segments */ #define TRACK_X_PADDING 2.0 /* distance from x boundary where tracking starts */ #define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */ @@ -336,8 +324,8 @@ #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ #define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ -#define PLOT_SPEEDS 1 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */ -#define PLOT_TRAJECTORIES 1 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */ +#define PLOT_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */ +#define PLOT_TRAJECTORIES 0 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */ #define VMAX_PLOT_SPEEDS 0.6 /* vertical scale of plot of obstacle speeds */ #define EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */ @@ -362,7 +350,7 @@ #define PMAX 1000.0 /* maximal force */ #define HASHX 120 /* size of hashgrid in x direction */ -#define HASHY 450 /* size of hashgrid in y direction */ +#define HASHY 60 /* size of hashgrid in y direction */ #define HASHMAX 100 /* maximal number of particles per hashgrid cell */ #define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ @@ -636,12 +624,15 @@ int thermostat_schedule(int i) double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], double qx[NMAXCIRCLES], double qy[NMAXCIRCLES], double qangle[NMAXCIRCLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], - double beta, int *nactive, int *nsuccess, int *nmove) + double beta, int *nactive, int *nsuccess, int *nmove, int initial_phase) { - double a, totalenergy = 0.0; + double a, totalenergy = 0.0, damping; static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; int j, move; + if (initial_phase) damping = INITIAL_DAMPING; + else damping = DAMPING; + #pragma omp parallel for private(j,xi,totalenergy,a,move) for (j=0; j #define MOVIE 0 /* set to 1 to generate movie */ +#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ /* General geometrical parameters */ @@ -179,6 +180,7 @@ #define E_SCALE 2500.0 /* scaling factor for energy representation */ #define LOG_SCALE 1.0 /* scaling factor for energy log representation */ #define LOG_SHIFT 0.0 /* shift of colors on log scale */ +#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ diff --git a/particle_billiard.c b/particle_billiard.c index 884236f..08c28ac 100644 --- a/particle_billiard.c +++ b/particle_billiard.c @@ -31,7 +31,7 @@ #include /* 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 */ @@ -57,9 +57,9 @@ #define B_DOMAIN 30 /* choice of domain shape */ #define CIRCLE_PATTERN 1 /* pattern of circles */ -#define POLYLINE_PATTERN 10 /* pattern of polyline */ +#define POLYLINE_PATTERN 11 /* pattern of polyline */ -#define ABSORBING_CIRCLES 1 /* set to 1 for circular scatterers to be absorbing */ +#define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ #define NMAXCIRCLES 100000 /* total number of circles (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 100000 /* total number of sides of polygonal line */ @@ -87,7 +87,7 @@ /* Simulation parameters */ -#define NPART 1000 /* number of particles */ +#define NPART 10000 /* number of particles */ // #define NPART 2000 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ #define LMAX 0.01 /* minimal segment length triggering resampling */ @@ -100,9 +100,9 @@ #define PRINT_COLLISION_NUMBER 0 /* set to 1 to print number of collisions */ #define TEST_ACTIVE 1 /* set to 1 to test whether particle is in billiard */ -#define TEST_INITIAL_COND 1 /* set to 1 to allow only initial conditions that pass a test */ +#define TEST_INITIAL_COND 0 /* set to 1 to allow only initial conditions that pass a test */ -#define NSTEPS 6000 /* number of frames of movie */ +#define NSTEPS 4500 /* number of frames of movie */ #define TIME 1500 /* time between movie frames, for fluidity of real-time simulation */ #define DPHI 0.00002 /* integration step */ #define NVID 25 /* number of iterations between images displayed on screen */ @@ -117,17 +117,17 @@ /* Colors and other graphical parameters */ -#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ -#define NCOLORS 16 /* number of colors */ +#define NCOLORS 128 /* number of colors */ #define COLORSHIFT 0 /* hue of initial color */ -#define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */ +#define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ #define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */ #define NSEG 100 /* number of segments of boundary */ // #define LENGTH 0.01 /* length of velocity vectors */ #define LENGTH 0.04 /* length of velocity vectors */ -#define BILLIARD_WIDTH 3 /* width of billiard */ -#define PARTICLE_WIDTH 3 /* width of particles */ +#define BILLIARD_WIDTH 2 /* width of billiard */ +#define PARTICLE_WIDTH 2 /* width of particles */ #define FRONT_WIDTH 3 /* width of wave front */ #define BLACK 1 /* set to 1 for black background */ @@ -141,10 +141,12 @@ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ -#define NXMAZE 8 /* width of maze */ -#define NYMAZE 8 /* height of maze */ +#define NXMAZE 16 /* width of maze */ +#define NYMAZE 16 /* height of maze */ +// #define NXMAZE 10 /* width of maze */ +// #define NYMAZE 10 /* height of maze */ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ -#define RAND_SHIFT 58 /* seed of random number generator */ +#define RAND_SHIFT 1 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ @@ -660,16 +662,25 @@ void print_left_right_part_number(double *configs[NPARTMAX], int active[NPARTMAX { char message[50]; int i, nleft = 0, nright = 0; - double rgb[3], x1, cosphi; + double rgb[3], x1, y1, cosphi, sinphi; + static int maxnleft = 0, maxnright = 0; /* count active particles, using the fact that absorbed particles have been given dummy coordinates */ - for (i=0; i= DUMMY_ABSORBING) { cosphi = (configs[i][6] - configs[i][4])/configs[i][3]; + sinphi = (configs[i][7] - configs[i][5])/configs[i][3]; x1 = configs[i][4] + configs[i][2]*cosphi; + y1 = configs[i][5] + configs[i][2]*cosphi; if (x1 < xmin) nleft++; - else if (x1 > xmax) nright++; + else if (x1 > xmax) + { + nright++; + printf("Detected leaving particle %i at (%.2f, %2f)\n", i, x1, y1); + } } + if (nleft > maxnleft) maxnleft = nleft; + if (nright > maxnright) maxnright = nright; hsl_to_rgb(0.0, 0.0, 0.0, rgb); @@ -677,11 +688,11 @@ void print_left_right_part_number(double *configs[NPARTMAX], int active[NPARTMAX erase_area(xr, yr - 0.03, 0.4, 0.12, rgb); glColor3f(1.0, 1.0, 1.0); - if (nleft > 1) sprintf(message, "%i particles", nleft); - else sprintf(message, "%i particle", nleft); + if (nleft > 1) sprintf(message, "%i particles", maxnleft); + else sprintf(message, "%i particle", maxnleft); write_text(xl, yl, message); - if (nright > 1) sprintf(message, "%i particles", nright); - else sprintf(message, "%i particle", nright); + if (nright > 1) sprintf(message, "%i particles", maxnright); + else sprintf(message, "%i particle", maxnright); write_text(xr, yr, message); } @@ -770,7 +781,7 @@ void animation() // init_drop_config(-0.5, 0.0, 0.2, 0.4, configs); - init_drop_config(-1.3, 0.2, -0.25*PI, 0.25*PI, configs); + init_drop_config(-1.15, 0.01, 0.00*PI, 0.3*PI, configs); // init_drop_config(-1.3, -0.1, 0.0, DPI, configs); // init_drop_config(1.4, 0.1, 0.0, DPI, configs); @@ -792,7 +803,7 @@ void animation() if (DRAW_BILLIARD) draw_billiard(); if (PRINT_PARTICLE_NUMBER) print_part_number(configs, active, XMIN + 0.1, YMIN + 0.1); else if (PRINT_LEFT_RIGHT_PARTICLE_NUMBER) - print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.1, XMAX - 0.45, YMIN + 0.1, YMAX + MAZE_XSHIFT, YMAX + MAZE_XSHIFT); + print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.05, XMAX - 0.45, YMIN + 0.05, XMIN + MAZE_XSHIFT, XMAX + MAZE_XSHIFT); else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1); glutSwapBuffers(); @@ -850,7 +861,7 @@ void animation() if (DRAW_BILLIARD) draw_billiard(); if (PRINT_PARTICLE_NUMBER) print_part_number(configs, active, XMIN + 0.1, YMIN + 0.1); else if (PRINT_LEFT_RIGHT_PARTICLE_NUMBER) - print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.1, XMAX - 0.45, YMIN + 0.1, YMAX + MAZE_XSHIFT, YMAX + MAZE_XSHIFT); + print_left_right_part_number(configs, active, XMIN + 0.1, YMIN + 0.05, XMAX - 0.45, YMIN + 0.05, YMIN + MAZE_XSHIFT, YMAX + MAZE_XSHIFT); // print_left_right_part_number(configs, XMIN + 0.1, YMIN + 0.1, XMAX - 0.45, YMIN + 0.1, YMIN + MAZE_XSHIFT, YMAX + MAZE_XSHIFT); else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1); diff --git a/rde.c b/rde.c index f32c86e..552b12a 100644 --- a/rde.c +++ b/rde.c @@ -41,6 +41,7 @@ #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 SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ /* General geometrical parameters */ @@ -63,6 +64,8 @@ // // #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 960 /* number of grid points on x axis */ +// // #define NY 540 /* 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 */ @@ -88,7 +91,7 @@ /* Choice of the billiard table */ -#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 197 /* 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 +99,8 @@ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ -#define MU 1.0 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.7 /* parameter controlling the dimensions of domain */ +#define MU 0.15 /* parameter controlling the dimensions of domain */ #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 */ @@ -150,7 +153,13 @@ #define SMOOTHEN_VORTICITY 1 /* set to 1 to smoothen vorticity field in Euler equation */ #define SMOOTHEN_PERIOD 10 /* period between smoothenings */ -#define SMOOTH_FACTOR 0.015 /* factor by which to smoothen */ +// #define SMOOTH_FACTOR 0.05 /* factor by which to smoothen */ +#define SMOOTH_FACTOR 0.03 /* factor by which to smoothen */ +// #define SMOOTH_FACTOR 0.015 /* factor by which to smoothen */ +// #define SMOOTH_FACTOR 0.01 /* factor by which to smoothen */ + +#define ADD_TRACERS 1 /* set to 1 to add tracer particles (for Euler equations) */ +#define N_TRACERS 1000 /* number of tracer particles */ #define T_OUT 2.0 /* outside temperature */ #define T_IN 0.0 /* inside temperature */ @@ -179,12 +188,14 @@ #define B_COND 1 +#define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */ + /* Parameters for length and speed of simulation */ -#define NSTEPS 4000 /* number of frames of movie */ -// #define NSTEPS 2500 /* number of frames of movie */ -#define NVID 50 /* number of iterations between images displayed on screen */ -// #define NVID 100 /* number of iterations between images displayed on screen */ +#define NSTEPS 2250 /* number of frames of movie */ +// #define NSTEPS 500 /* number of frames of movie */ +// #define NVID 90 /* number of iterations between images displayed on screen */ +#define NVID 120 /* 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 */ @@ -208,6 +219,8 @@ #define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ +#define DRAW_PERIODICISED 1 /* set to 1 to repeat wave periodically in x and y directions */ + /* Plot type - color scheme */ #define CPLOT 52 @@ -253,8 +266,8 @@ /* Color schemes, see list in global_pdes.c */ -#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ @@ -264,7 +277,7 @@ #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 0.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_AMPLITUDE 1.5 /* 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) */ @@ -278,18 +291,20 @@ #define HUEMEAN 359.0 /* mean value of hue for color scheme C_HUE */ #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 0.75 /* scaling factor for energy log representation */ +#define FLUX_SCALE 100.0 /* scaling factor for energy representation */ +#define LOG_SCALE 0.5 /* scaling factor for energy log representation */ #define LOG_SHIFT 1.0 +#define LOG_MIN 1.0e-3 /* floor value for log vorticity plot */ #define NXMAZE 7 /* width of maze */ #define NYMAZE 7 /* height of maze */ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ -#define RAND_SHIFT 24 /* seed of random number generator */ +#define RAND_SHIFT 0 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ #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 2.5 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE_B 3.0 /* 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 */ @@ -320,12 +335,12 @@ double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "lig double observer[3] = {8.0, 8.0, 8.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.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 Z_SCALING_FACTOR 0.08 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define XY_SCALING_FACTOR 1.7 /* 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.0 /* overall y shift for REP_PROJ_3D representation */ -#define BORDER_PADDING 2 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */ +#define YSHIFT_3D 0.1 /* overall y shift for REP_PROJ_3D representation */ +#define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */ /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ @@ -514,8 +529,8 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i { nabla_psi = (double *)malloc(2*NX*NY*sizeof(double)); nabla_omega = (double *)malloc(2*NX*NY*sizeof(double)); - compute_gradient_euler(phi_in[0], nabla_psi); - compute_gradient_euler(phi_in[1], nabla_omega); + compute_gradient_euler(phi_in[0], nabla_psi, EULER_GRADIENT_YSHIFT); + compute_gradient_euler(phi_in[1], nabla_omega, 0.0); dx = (XMAX-XMIN)/((double)NX); if (SMOOTHEN_VORTICITY) /* beta: try to reduce formation of ripples */ @@ -622,9 +637,15 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i } case (E_EULER_INCOMP): { - phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*stiffness*(delta_phi[0][i*NY+j] + phi_in[1][i*NY+j]*dx*dx); - phi_out[1][i*NY+j] = phi_in[1][i*NY+j] - intstep*K_EULER*(nabla_omega[i*NY+j]*nabla_psi[NX*NY+i*NY+j]); - phi_out[1][i*NY+j] += intstep*K_EULER*(nabla_omega[NX*NY+i*NY+j]*nabla_psi[i*NY+j]); +// if ((j > 1)&&(j < NY - 1)) + { + phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*stiffness*(delta_phi[0][i*NY+j] + phi_in[1][i*NY+j]*dx*dx); +// phi_out[0][i*NY+j] += intstep*EULER_GRADIENT_YSHIFT; + phi_out[1][i*NY+j] = phi_in[1][i*NY+j] - intstep*K_EULER*(nabla_omega[i*NY+j]*nabla_psi[NX*NY+i*NY+j]); + phi_out[1][i*NY+j] += intstep*K_EULER*(nabla_omega[NX*NY+i*NY+j]*nabla_psi[i*NY+j]); + +// if ((i == 0)&&(j%10 == 0)) printf("j = %i, psi = %.5lg\n", j, phi_out[0][i*NY+j]); + } break; } } @@ -664,7 +685,7 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i } } -void evolve_wave(double *phi[NFIELDS], double *phi_tmp[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY], double vector_potential_field[2*NX*NY]) +void evolve_wave(double *phi[NFIELDS], double *phi_tmp[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY], double vector_potential_field[2*NX*NY]) /* time step of field evolution */ { evolve_wave_half(phi, phi_tmp, xy_in, potential_field, vector_potential_field); @@ -672,6 +693,63 @@ void evolve_wave(double *phi[NFIELDS], double *phi_tmp[NFIELDS], short int xy_in } +void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int time, int nsteps, double step) +/* time steps of tracer particle evolution (for Euler equation) */ +{ + int tracer, i, j, t, ij[2], iplus, jplus; + double x, y, xy[2], vx, vy; + + step = 0.2; + + for (tracer = 0; tracer < N_TRACERS; tracer++) + { + x = tracers[time*2*N_TRACERS + 2*tracer]; + y = tracers[time*2*N_TRACERS + 2*tracer + 1]; + +// printf("Tracer %i position (%.2f, %.2f)\n", tracer, x, y); + + for (t=0; t 0.0)&&(v < 0.1)) +// { +// vx = vx*0.1/v; +// vy = vy*0.1/v; +// } + +// printf("(i, j) = (%i, %i), Tracer %i velocity (%.6f, %.6f)\n", i, j, tracer, vx, vy); + + x += vx*step; + y += vy*step; + } +// printf("Tracer %i velocity (%.2f, %.2f)\n", tracer, vx, vy); + + if (x > XMAX) x += (XMIN - XMAX); + if (x < XMIN) x += (XMAX - XMIN); + if (y > YMAX) y += (YMIN - YMAX); + if (y < YMIN) y += (YMAX - YMIN); + + if (time+1 < NSTEPS) + { + tracers[(time+1)*2*N_TRACERS + 2*tracer] = x; + tracers[(time+1)*2*N_TRACERS + 2*tracer + 1] = y; + } + } +} + + void print_level(int level) { double pos[2]; @@ -852,7 +930,7 @@ void animation() { double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY; - double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field; + double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers; short int *xy_in; int i, j, k, s, nvid, field; static int counter = 0; @@ -886,6 +964,8 @@ void animation() initialize_potential(potential_field); initialize_vector_potential(vector_potential_field); } + + if (ADD_TRACERS) tracers = (double *)malloc(2*NSTEPS*N_TRACERS*sizeof(double)); dx = (XMAX-XMIN)/((double)NX); intstep = DT/(dx*dx); @@ -908,10 +988,10 @@ void animation() // init_fermion_state(-0.5, 0.5, 2.0, 0.0, 0.1, phi, xy_in); // init_boson_state(-0.5, 0.5, 2.0, 0.0, 0.1, phi, xy_in); -// init_vortex_state(0.4, 0.0, 0.1, phi, xy_in); -// add_vortex_state(-0.4, 0.0, 0.1, phi, xy_in); +// init_shear_flow(1.0, 0.02, 0.15, 1, 1, phi, xy_in); +// init_laminar_flow(1.0, 0.1, 0.5, 0.0, phi, xy_in); - init_shear_flow(1.0, 0.02, 0.03, 1, 1, phi, xy_in); + init_shear_flow(-1.0, 0.0, 0.1, 1, 1, 0.0, phi, xy_in); init_cfield_rde(phi, xy_in, CPLOT, rde, 0); if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT, rde, 0); @@ -922,6 +1002,12 @@ void animation() if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT_B, rde, 1); } + if (ADD_TRACERS) for (i=0; i= INITIAL_TIME)&&(DOUBLE_MOVIE)) { draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B); + if (ADD_TRACERS) draw_tracers(phi, tracers, i, 0, 1.0); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); @@ -1062,6 +1158,7 @@ void animation() if (DOUBLE_MOVIE) { draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); + if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 0, 1.0); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); @@ -1072,6 +1169,7 @@ void animation() { fade_value = 1.0 - (double)i/(double)MID_FRAMES; draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0); + if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value); @@ -1079,6 +1177,7 @@ void animation() save_frame_counter(NSTEPS + i + 1); } draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B); + if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 0, 1.0); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); glutSwapBuffers(); @@ -1087,6 +1186,7 @@ void animation() { fade_value = 1.0 - (double)i/(double)END_FRAMES; draw_wave_rde(1, phi, xy_in, rde, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0); + if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value); glutSwapBuffers(); save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i); @@ -1099,6 +1199,7 @@ void animation() { fade_value = 1.0 - (double)i/(double)END_FRAMES; draw_wave_rde(0, phi, xy_in, rde, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0); + if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 1, fade_value); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value); glutSwapBuffers(); save_frame_counter(NSTEPS + 1 + counter + i); @@ -1120,6 +1221,7 @@ void animation() free(potential_field); free(vector_potential_field); } + if (ADD_TRACERS) free(tracers); printf("Time %.5lg\n", time); diff --git a/schrodinger.c b/schrodinger.c index af33c4d..9a1cafa 100644 --- a/schrodinger.c +++ b/schrodinger.c @@ -36,6 +36,7 @@ #include #define MOVIE 0 /* set to 1 to generate movie */ +#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ /* General geometrical parameters */ diff --git a/sub_lj.c b/sub_lj.c index 027411a..d4fa448 100644 --- a/sub_lj.c +++ b/sub_lj.c @@ -955,7 +955,7 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) } } -void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, t_segment segment[NMAXSEGMENTS]) +void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, t_segment segment[NMAXSEGMENTS], int group) /* add four segments forming a rectangle, specified by two adjacent corners and width */ { double tx, ty, ux, uy, norm, x3, y3, x4, y4; @@ -993,7 +993,11 @@ void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, d segment[n+3].x2 = x1; segment[n+3].y2 = y1; - for (i=0; i<4; i++) segment[n+i].concave = 1; + for (i=0; i<4; i++) + { + segment[n+i].concave = 1; + segment[n+i].group = group; + } nsegments += 4; } else printf("Warning: NMAXSEGMENTS too small\n"); @@ -1050,6 +1054,38 @@ void add_rectangle_to_segments(double x1, double y1, double x2, double y2, t_seg } +void add_circle_to_segments(double x, double y, double r, int nsegs, double angle0, t_segment segment[NMAXSEGMENTS], int group) +/* add segments forming a circle to linear obstacle configuration */ +{ + int i, n = nsegments, nplus, nminus; + double angle; + + angle = DPI/(double)nsegs; + + if (nsegments + nsegs < NMAXSEGMENTS) + { + for (i=0; i nsegments - 1) iplus = 0; } angle = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID; angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID; if (angle2 < angle) angle2 += DPI; segment[i].angle1 = angle; segment[i].angle2 = angle2; + + printf("i = %i, iplus = %i, iminus = %i, angle1 = %.0f, angle2 = %.0f\n", i, iplus, iminus, angle*360.0/DPI, angle2*360.0/DPI); } /* make copy of initial values in case of rotation/translation */ @@ -1868,6 +1956,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS]) // printf("Segment %i: (x1, y1) = (%.3lg,%.3lg), (x2, y2) = (%.3lg,%.3lg)\n (nx, ny) = (%.3lg,%.3lg), c = %.3lg, length = %.3lg\n", i, segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2, segment[i].nx, segment[i].ny, segment[i].c, segment[i].length); // if (segment[i].concave) printf("Concave with angles %.3lg Pi, %.3lg Pi\n", segment[i].angle1/PI, segment[i].angle2/PI); // } +// sleep(4); } int in_rocket(double x, double y, int rocket_shape) @@ -2044,6 +2133,13 @@ int in_segment_region(double x, double y) else if (vabs(y) > 0.1*LAMBDA + padding) return(1); else return(0); } + case (S_EXT_CIRCLE_RECT): + { + padding = 0.1; + if ((vabs(x - SEGMENTS_X0) < LAMBDA + padding)&&(vabs(y - SEGMENTS_Y0) < 0.1*LAMBDA + padding)) return(0); + else if (module2(x + SEGMENTS_X0, y - SEGMENTS_Y0) < 0.5*LAMBDA + padding) return(0); + else return(1); + } default: return(1); } } @@ -3192,6 +3288,33 @@ void set_segment_group_color(int group, double lum, double rgb[3]) } +void draw_altitude_lines() +{ + int i, i1, i2; + double x, y; + + glColor3f(0.5, 0.5, 0.5); + glLineWidth(1); + + i1 = (int)(YMIN + ytrack) - 1.0; + i2 = (int)(YMAX + ytrack) + 1.0; + + for (i = i1; i < i2; i++) + { + y = (double)i - ytrack; + draw_line(BCXMIN, y, XMAX - 1.8, y); + } + + i1 = (int)(XMIN + xtrack) - 1.0; + i2 = (int)(XMAX - 1.8 + xtrack) + 1.0; + + for (i = i1; i < i2; i++) + { + x = (double)i - xtrack; + if (x < XMAX - 1.8) draw_line(x, YMIN, x, YMAX); + } +} + void draw_one_triangle(t_particle particle, int same_table[9*HASHMAX], int p, int q, int nsame) { double x, y, dx, dy; @@ -3450,6 +3573,9 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta) } + /* draw "altitude lines" */ + if (ALTITUDE_LINES) draw_altitude_lines(); + /* draw the bonds first */ if ((DRAW_BONDS)||(plot == P_BONDS)) { @@ -4996,7 +5122,7 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr pangle[i] = particle[i].omega; if ((PLOT == P_INITIAL_POS)||(PLOT_B == P_INITIAL_POS)) - particle[i].color_hue = 360.0*(particle[i].yc - YMIN)/(YMAX - YMIN); + particle[i].color_hue = 360.0*(particle[i].yc - INITYMIN)/(INITYMAX - INITYMIN); } /* initialize dummy values in case particles are added */ for (i=ncircles; i < NMAXCIRCLES; i++) @@ -5469,6 +5595,17 @@ void draw_trajectory_plot(t_group_data *group_speeds, int i) first = 0; } + if (ALTITUDE_LINES) + { + glLineWidth(1); + glColor3f(0.5, 0.5, 0.5); + for (j=1; j<(int)scaley + 1; j++) + { + y1 = plot_coord((double)j/scaley, plotymin, plotymax); + if (y1 < plotymax) draw_line(plotxmin, y1, plotxmax + 0.05, y1); + } + } + glLineWidth(2); /* plot trajectories */ @@ -5514,7 +5651,7 @@ void draw_trajectory_plot(t_group_data *group_speeds, int i) // sprintf(message, "%.1f", VMAX_PLOT_SPEEDS); sprintf(message, "y"); - write_text_fixedwidth(plotxmin - 0.1, plotymax - 0.01, message); + write_text_fixedwidth(plotxmin - 0.1, plotymax - 0.005, message); sprintf(message, "x"); write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.1, message); @@ -5532,6 +5669,9 @@ void init_segment_group(t_segment segment[NMAXSEGMENTS], int group, t_group_segm xc += 0.5*(segment[i].x1 + segment[i].x2); yc += 0.5*(segment[i].y1 + segment[i].y2); +// printf("Segment group data %i: z1 = (%.3lg, %.3lg)\n z2 = (%.3lg, %.3lg)\n zc = (%.3lg, %.3lg)\n\n", group, +// segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2, xc, yc); + nseg_group++; } diff --git a/sub_maze.c b/sub_maze.c index 6962d17..d2b441e 100644 --- a/sub_maze.c +++ b/sub_maze.c @@ -12,6 +12,8 @@ typedef struct short int north, east, south, west; /* closed walls */ short int active; /* takes value 1 if currently active in RW path */ short int tested; /* takes value 1 if tested */ + short int connected; /* takes value 1 if connected to exit */ + short int closed; /* takes value 1 if no untested neighbours */ } t_maze; @@ -122,6 +124,8 @@ void init_maze_graph(t_maze maze[NXMAZE*NYMAZE]) n = nmaze(i, j); maze[n].active = 0; maze[n].tested = 0; + maze[n].connected = 0; + maze[n].closed = 0; maze[n].north = 1; maze[n].east = 1; maze[n].south = 1; @@ -129,14 +133,17 @@ void init_maze_graph(t_maze maze[NXMAZE*NYMAZE]) } } -int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0) +int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0, int *path, int *pathlength) /* find a random walk path in the maze */ +/* returns 0 or 1 depending on whether path reaches a tested cell or a deadend */ { - int active_counter = 0, i, n = n0, npaths, inext, nextcell, trial, nnext; + int active_counter = 0, i, n = n0, npaths, inext, nextcell, trial, nnext, deadend = 1, length = 0; int next_table[4]; /* contruct random walk */ npaths = maze[n].nneighb; + path[0] = n0; + // while ((npaths > 0)&&(!maze[n].tested)) while ((npaths > 0)) { @@ -148,7 +155,7 @@ int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0) for (i=0; i pathlength) j = 0; + printf("j = %i\n", j); +// while (deadend) + if (!maze[path[j]].connected) deadend = find_maze_path(maze, path[j], newpath, &newpathlength); + if (!deadend) for (n=0; n pathlength) j = 0; + printf("j = %i\n", j); +// while (deadend) + if (!maze[path[j]].connected) deadend = find_maze_path(maze, path[j], newpath, &newpathlength); + if (!deadend) for (n=0; n0)||(j!=NYMAZE/2))&&(maze[n].west)) + { + polyline[nsides].x1 = x1; + polyline[nsides].y1 = y1; + polyline[nsides].x2 = x1; + polyline[nsides].y2 = y1 + dy; + polyline[nsides].angle = PID; + nsides++; + } + + + if (maze[n].south) + { + polyline[nsides].x1 = x1; + polyline[nsides].y1 = y1; + polyline[nsides].x2 = x1 + dx; + polyline[nsides].y2 = y1; + polyline[nsides].angle = 0.0; + nsides++; + } + + } + + /* 45 degrees parts */ + for (i=0; i 0.0) + if ((conf[0] >= DUMMY_ABSORBING)&&(x2 > XMAX)) +// if (conf[0] >= DUMMY_ABSORBING) { active[i] = 1; if (time > tmaxmax) tmaxmax = time; diff --git a/sub_perco_3d.c b/sub_perco_3d.c new file mode 100644 index 0000000..943b18f --- /dev/null +++ b/sub_perco_3d.c @@ -0,0 +1,185 @@ +/* routines for 3D representation of percolation clusters, taken from sub_wave_3d_rde.c */ + +void xyz_to_xy(double x, double y, double z, double xy_out[2]) +{ + int i; + double s, t, xinter[3]; + static double n2, m2, d, sm2, sn2, v[3], h[2], plane_ratio = 0.5; + static int first = 1; + + if (((first)&&(REPRESENTATION_3D == REP_PROJ_3D))||(reset_view)) + { + m2 = observer[0]*observer[0] + observer[1]*observer[1]; + n2 = m2 + observer[2]*observer[2]; + d = plane_ratio*n2; + sm2 = sqrt(m2); + sn2 = sqrt(n2); + h[0] = observer[1]/sm2; + h[1] = -observer[0]/sm2; + v[0] = -observer[0]*observer[2]/(sn2*sm2); + v[1] = -observer[1]*observer[2]/(sn2*sm2); + v[2] = m2/(sn2*sm2); + first = 0; + reset_view = 0; +// printf("h = (%.3lg, %.3lg)\n", h[0], h[1]); +// printf("v = (%.3lg, %.3lg, %.3lg)\n", v[0], v[1], v[2]); + } + + switch (REPRESENTATION_3D) { + case (REP_AXO_3D): + { + for (i=0; i<2; i++) + xy_out[i] = x*u_3d[i] + y*v_3d[i] + z*w_3d[i]; + break; + } + case (REP_PROJ_3D): + { + if (z > ZMAX_FACTOR*n2) z = ZMAX_FACTOR*n2; + z *= Z_SCALING_FACTOR; + s = observer[0]*x + observer[1]*y + observer[2]*z; + t = (d - s)/(n2 - s); + xinter[0] = t*observer[0] + (1.0-t)*x; + xinter[1] = t*observer[1] + (1.0-t)*y; + xinter[2] = t*observer[2] + (1.0-t)*z; + + xy_out[0] = XSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*h[0] + xinter[1]*h[1]); + xy_out[1] = YSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*v[0] + xinter[1]*v[1] + xinter[2]*v[2]); + break; + } + } +} + +// void draw_vertex_ij(int i, int j) +// { +// double xy[2]; +// +// ij_to_xy(i, j, xy); +// // if (xy[1] > 0.0) printf("(i,j) = (%i,%i), (x,y) = (%.2lg,%.2lg)\n", i, j, xy[0], xy[1]); +// glVertex2d(xy[0], xy[1]); +// } + + +void draw_vertex_xyz(double xy[2], double z) +{ + double xy_screen[2]; + + xyz_to_xy(xy[0], xy[1], z, xy_screen); + glVertex2d(xy_screen[0], xy_screen[1]); +} + +void draw_vertex_x_y_z(double x, double y, double z) +{ + double xy_screen[2]; + + xyz_to_xy(x, y, z, xy_screen); + glVertex2d(xy_screen[0], xy_screen[1]); +} + +void draw_line_3d(double x1, double y1, double z1, double x2, double y2, double z2) +{ + glBegin(GL_LINE_LOOP); + draw_vertex_x_y_z(x1, y1, z1); + draw_vertex_x_y_z(x2, y2, z2); + glEnd(); +} + + +double angle_lum(double cosangle) +{ + double mid_lum = 0.6; + + if (cosangle > 0.0) return(mid_lum + (1.0 - mid_lum)*cosangle); + else return (mid_lum + 0.8*(1.0 - mid_lum)*cosangle); +} + +void draw_cube(double x, double y, double z, double size, double rgb[3]) +/* draw a cube */ +{ + double lum_factor; + + /* front/back face */ + if (observer[0] > 0.0) + { + lum_factor = angle_lum(light[0]); + glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor); + glBegin(GL_TRIANGLE_FAN); + draw_vertex_x_y_z(x + size, y, z); + draw_vertex_x_y_z(x + size, y + size, z); + draw_vertex_x_y_z(x + size, y + size, z + size); + draw_vertex_x_y_z(x + size, y, z + size); + glEnd(); + } + else + { +// cosangle = -light[0]; + lum_factor = angle_lum(-light[0]); + glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor); + glBegin(GL_TRIANGLE_FAN); + draw_vertex_x_y_z(x, y, z); + draw_vertex_x_y_z(x, y + size, z); + draw_vertex_x_y_z(x, y + size, z + size); + draw_vertex_x_y_z(x, y, z + size); + glEnd(); + } + + + /* right/left face */ + if (observer[1] > 0.0) + { +// cosangle = light[1]; + lum_factor = angle_lum(light[1]); + glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor); + glBegin(GL_TRIANGLE_FAN); + draw_vertex_x_y_z(x + size, y + size, z); + draw_vertex_x_y_z(x, y + size, z); + draw_vertex_x_y_z(x, y + size, z + size); + draw_vertex_x_y_z(x + size, y + size, z + size); + glEnd(); + } + else + { +// cosangle = -light[1]; + lum_factor = angle_lum(-light[1]); + glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor); + glBegin(GL_TRIANGLE_FAN); + draw_vertex_x_y_z(x + size, y, z); + draw_vertex_x_y_z(x, y, z); + draw_vertex_x_y_z(x, y, z + size); + draw_vertex_x_y_z(x + size, y, z + size); + glEnd(); + } + + + /* top face */ +// cosangle = light[2]; + lum_factor = angle_lum(light[2]); + glColor3f(rgb[0]*lum_factor, rgb[1]*lum_factor, rgb[2]*lum_factor); + glBegin(GL_TRIANGLE_FAN); + draw_vertex_x_y_z(x + size, y, z + size); + draw_vertex_x_y_z(x + size, y + size, z + size); + draw_vertex_x_y_z(x, y + size, z + size); + draw_vertex_x_y_z(x, y, z + size); + glEnd(); +} + +void viewpoint_schedule(int i) +/* change position of observer */ +{ + int j; + double angle, ca, sa; + static double observer_initial[3]; + static int first = 1; + + if (first) + { + for (j=0; j<3; j++) observer_initial[j] = observer[j]; + first = 0; + } + + angle = (ROTATE_ANGLE*DPI/360.0)*(double)i/(double)NSTEPS; + ca = cos(angle); + sa = sin(angle); + observer[0] = ca*observer_initial[0] - sa*observer_initial[1]; + observer[1] = sa*observer_initial[0] + ca*observer_initial[1]; + printf("Angle %.3lg, Observer position (%.3lg, %.3lg, %.3lg)\n", angle, observer[0], observer[1], observer[2]); +} diff --git a/sub_rde.c b/sub_rde.c index 7ba6dc9..bdb6b86 100644 --- a/sub_rde.c +++ b/sub_rde.c @@ -318,7 +318,7 @@ void init_vortex_state(double x, double y, double scale, double *phi[NFIELDS], s /* phi[0] is stream function, phi[1] is vorticity */ { int i, j; - double xy[2], dist2, module, phase, scale2; + double xy[2], dist2, module, phase, scale2, sign; // scale2 = scale*scale; for (i=0; i 0.0) sign = 1.0; + else sign = -1.0; + + phi[1][i*NY+j] += sign*module; + phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */ } else { @@ -348,7 +351,7 @@ void add_vortex_state(double x, double y, double scale, double *phi[NFIELDS], sh /* phi[0] is stream function, phi[1] is vorticity */ { int i, j; - double xy[2], dist2, module, phase, scale2; + double xy[2], dist2, module, phase, scale2, sign; // scale2 = scale*scale; for (i=0; i 0.0) sign = 1.0; + else sign = -1.0; + + phi[1][i*NY+j] += sign*module; + phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */ } else { @@ -374,7 +380,7 @@ void add_vortex_state(double x, double y, double scale, double *phi[NFIELDS], sh } -void init_shear_flow(double amp, double delta, double rho, int nx, int ny, double *phi[NFIELDS], short int xy_in[NX*NY]) +void init_shear_flow(double amp, double delta, double rho, int nx, int ny, double yshift, double *phi[NFIELDS], short int xy_in[NX*NY]) /* initialise field with a shear flow */ /* phi[0] is stream function, phi[1] is vorticity */ /* amp is global amplitude */ @@ -395,7 +401,7 @@ void init_shear_flow(double amp, double delta, double rho, int nx, int ny, doubl ij_to_xy(i, j, xy); xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]); - y1 = xy[1]*(double)ny/YMAX; + y1 = xy[1]*(double)ny/YMAX - yshift; while (y1 > 1.0) y1 -= 2.0; while (y1 < -1.0) y1 += 2.0; @@ -424,6 +430,39 @@ void init_shear_flow(double amp, double delta, double rho, int nx, int ny, doubl } } +void init_laminar_flow(double amp, double modulation, double period, double yshift, double *phi[NFIELDS], short int xy_in[NX*NY]) +/* initialise field with a laminar flow in x direction */ +/* phi[0] is stream function, phi[1] is vorticity */ +/* amp is global amplitude */ +{ + int i, j; + double xy[2], y1, a, b, f, cplus, cminus; + + a = period*PI/YMAX; +// a = PID/YMAX; + + for (i=0; i 0) + { + z_mid = compute_interpolated_colors_rde(i, j, rde, potential, palette, cplot, + rgb_e, rgb_w, rgb_n, rgb_s, &z_sw, &z_se, &z_nw, &z_ne, + fade, fade_value, movie); + ij_to_xy(i, j, xy_sw); + ij_to_xy(i+1, j, xy_se); + ij_to_xy(i, j+1, xy_nw); + ij_to_xy(i+1, j+1, xy_ne); + + for (k=0; k<2; k++) xy_mid[k] = 0.25*(xy_sw[k] + xy_se[k] + xy_nw[k] + xy_ne[k]); + + if (AMPLITUDE_HIGH_RES == 1) + { + glBegin(GL_TRIANGLE_FAN); + glColor3f(rgb_w[0], rgb_w[1], rgb_w[2]); + draw_vertex_xyz_shift(xy_mid, z_mid, shiftx, shifty); + draw_vertex_xyz_shift(xy_nw, z_nw, shiftx, shifty); + draw_vertex_xyz_shift(xy_sw, z_sw, shiftx, shifty); + + glColor3f(rgb_s[0], rgb_s[1], rgb_s[2]); + draw_vertex_xyz_shift(xy_se, z_se, shiftx, shifty); + + glColor3f(rgb_e[0], rgb_e[1], rgb_e[2]); + draw_vertex_xyz_shift(xy_ne, z_ne, shiftx, shifty); + + glColor3f(rgb_n[0], rgb_n[1], rgb_n[2]); + draw_vertex_xyz_shift(xy_nw, z_nw, shiftx, shifty); + glEnd (); + } + } + else + { + glColor3f(rde[i*NY+j].rgb[0], rde[i*NY+j].rgb[1], rde[i*NY+j].rgb[2]); + + glBegin(GL_TRIANGLE_FAN); + ij_to_xy(i, j, xy); + draw_vertex_xyz_shift(xy, adjust_field(*rde[i*NY+j].p_zfield[movie], potential[i*NY+j]), shiftx, shifty); + ij_to_xy(i+1, j, xy); + draw_vertex_xyz_shift(xy, adjust_field(*rde[(i+1)*NY+j].p_zfield[movie], potential[(i+1)*NY+j]), shiftx, shifty); + ij_to_xy(i+1, j+1, xy); + draw_vertex_xyz_shift(xy, adjust_field(*rde[(i+1)*NY+j+1].p_zfield[movie], potential[(i+1)*NY+j+1]), shiftx, shifty); + ij_to_xy(i, j+1, xy); + draw_vertex_xyz_shift(xy, adjust_field(*rde[i*NY+j+1].p_zfield[movie], potential[i*NY+j+1]), shiftx, shifty); + glEnd (); + } + } +} void draw_wave_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY], @@ -1742,6 +1850,27 @@ void draw_wave_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t } +void draw_periodicised_wave_3d(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], double potential[NX*NY], + int zplot, int cplot, int palette, int fade, double fade_value) +/* a variant where the wave is repeated periodically in x and y directions (beta) */ +{ + int i, j, shiftx, shifty; + double observer_angle; + + blank(); + if (DRAW_BILLIARD) draw_billiard_3d(fade, fade_value); + + if (!ROTATE_VIEW) + { + for (shiftx = -1; shiftx < 2; shiftx++) + for (shifty = -2; shifty < 2; shifty++) + for (i=BORDER_PADDING; i NX-1) ij[0] = NX-1; + if (ij[1] < 0) ij[1] = 0; + if (ij[1] > NY-1) ij[1] = NY-1; +} + + void xy_to_pos(double x, double y, double pos[2]) /* convert (x,y) position to double-valued position in table representing wave */ { @@ -1606,13 +1635,14 @@ int compute_maze_coordinates(t_rectangle polyrect[NMAXPOLY], int closed) { t_maze* maze; int i, j, n, nsides = 0, ropening; - double dx, dy, x1, y1, padding = 0.02, pos[2], width = 0.02; + double dx, dy, x1, y1, x0, padding = 0.02, pos[2], width = 0.02; maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze)); init_maze(maze); /* build walls of maze */ +// x0 = LAMBDA - 1.0; dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE); dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE); @@ -1918,6 +1948,7 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height; int i, j, k, k1, k2, condition = 0, m; static int first = 1, nsides; + static double h, hh, ra, rb; switch (b_domain) { case (D_NOTHING): @@ -1937,6 +1968,28 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ else return(0); break; } + case (D_EXT_ELLIPSE): + { + if (x*x/(LAMBDA*LAMBDA) + y*y/(MU*MU) > 1.0) return(1); + else return(0); + break; + } + case (D_EXT_ELLIPSE_CURVED): + { + y1 = y + 0.4*x*x; + if (x*x/(LAMBDA*LAMBDA) + y1*y1/(MU*MU) > 1.0) return(1); + else return(0); + break; + } + case (D_EXT_ELLIPSE_CURVED_BDRY): + { + if (y > YMAX - 0.05) return(0); + if (y < YMIN + 0.05) return(0); + y1 = y + 0.4*x*x; + if (x*x/(LAMBDA*LAMBDA) + y1*y1/(MU*MU) > 1.0) return(1); + else return(0); + break; + } case (D_STADIUM): { if ((x > -0.5*LAMBDA)&&(x < 0.5*LAMBDA)&&(y > -1.0)&&(y < 1.0)) return(1); @@ -2485,6 +2538,49 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ if ((x > polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0); return(1); } + case (D_CHESSBOARD): + { + i = (int)(vabs(x)/LAMBDA + 0.5); + j = (int)(vabs(y)/LAMBDA + 0.5); + if ((i+j)%2 == 0) return(1); + else return(0); + } + case (D_TRIANGLE_TILES): + { + if (first) + { + h = LAMBDA/(2.0*sqrt(3.0)); + hh = h*3.0; + first = 0; + } + i = (int)((y + h)/hh + 10.0); + y1 = sin(DPI/3.0)*x - 0.5*y; + j = (int)((y1 + h)/hh + 10.0); + y1 = sin(-DPI/3.0)*x -0.5*y; + k = (int)((y1 + h)/hh + 10.0); + if ((i+j+k)%2 == 0) return(1); + else return(0); + } + case (D_HEX_TILES): + { + if (first) + { + ra = -1.0/sqrt(3.0); + rb = -2.0*ra; + first = 0; + } + x1 = (x + ra*y)/LAMBDA + 30.0; + y1 = rb*y/LAMBDA + 30.0; + + x1 = x1 - (double)(3*(int)(x1/3.0)); + y1 = y1 - (double)(3*(int)(y1/3.0)); + + if ((x1 > 2.0)&&(y1 < 1.0)) return(1); + if ((x1 < 1.0)&&(y1 > 2.0)) return(1); + if (x1 + y1 < 1.0) return(1); + if (x1 + y1 > 5.0) return(1); + return(0); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -2675,11 +2771,30 @@ void tvertex_lineto(t_vertex z) } +void hex_transfo(double u, double v, double *x, double *y) +/* linear transformation of plane used for hex tiles */ +{ + static double ra, rb; + static int first = 1; + + if (first) + { + ra = 0.5; + rb = 0.5*sqrt(3.0); + first = 0; + } + + *x = u + ra*v; + *y = rb*v; +} + + void draw_billiard(int fade, double fade_value) /* draws the billiard boundary */ { - double x0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax, height; - int i, j, k, k1, k2, mr2; + double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax, height; + int i, j, k, k1, k2, mr2, ntiles; static int first = 1, nsides; + static double h, hh, sqr3; if (fade) { @@ -2738,6 +2853,67 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound } break; } + case (D_EXT_ELLIPSE): + { + glBegin(GL_LINE_LOOP); + for (i=0; i<=NSEG; i++) + { + phi = (double)i*DPI/(double)NSEG; + x = LAMBDA*cos(phi); + y = MU*sin(phi); + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + glEnd (); + + /* draw foci */ + if (FOCI) + { + if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value); + else glColor3f(0.3, 0.3, 0.3); + x0 = sqrt(LAMBDA*LAMBDA-MU*MU); + + glLineWidth(2); + glEnable(GL_LINE_SMOOTH); + + draw_circle(x0, 0.0, r, NSEG); + draw_circle(-x0, 0.0, r, NSEG); + } + break; + } + case (D_EXT_ELLIPSE_CURVED): + { + glBegin(GL_LINE_LOOP); + for (i=0; i<=NSEG; i++) + { + phi = (double)i*DPI/(double)NSEG; + x = LAMBDA*cos(phi); + y = MU*sin(phi) - 0.4*x*x; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + glEnd (); + + break; + } + case (D_EXT_ELLIPSE_CURVED_BDRY): + { + glBegin(GL_LINE_LOOP); + for (i=0; i<=NSEG; i++) + { + phi = (double)i*DPI/(double)NSEG; + x = LAMBDA*cos(phi); + y = MU*sin(phi) - 0.4*x*x; + xy_to_pos(x, y, pos); + glVertex2d(pos[0], pos[1]); + } + glEnd (); + + draw_line(XMIN, YMAX - 0.05, XMAX, YMAX - 0.05); + draw_line(XMIN, YMIN + 0.05, XMAX, YMIN + 0.05); + + break; + } case (D_STADIUM): { glBegin(GL_LINE_LOOP); @@ -3631,6 +3807,109 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound draw_filled_rectangle(polyrect[i].x1, polyrect[i].y1, polyrect[i].x2, polyrect[i].y2); break; } + case (D_CHESSBOARD): + { + glLineWidth(BOUNDARY_WIDTH); + x = 0.5*LAMBDA; + while (x < XMAX) + { + draw_line(x, YMIN, x, YMAX); + x += LAMBDA; + } + x = -0.5*LAMBDA; + while (x > XMIN) + { + draw_line(x, YMIN, x, YMAX); + x -= LAMBDA; + } + y = 0.5*LAMBDA; + while (y < YMAX) + { + draw_line(XMIN, y, XMAX, y); + y += LAMBDA; + } + y = -0.5*LAMBDA; + while (y > YMIN) + { + draw_line(XMIN, y, XMAX, y); + y -= LAMBDA; + } + + break; + } + case (D_TRIANGLE_TILES): + { + if (first) + { + h = LAMBDA/(2.0*sqrt(3.0)); + hh = h*3.0; + sqr3 = sqrt(3.0); + first = 0; + } + glLineWidth(BOUNDARY_WIDTH); + y = -h; + while (y < YMAX) + { + draw_line(XMIN, y, XMAX, y); + y += hh; + } + y = -h; + while (y > YMIN) + { + draw_line(XMIN, y, XMAX, y); + y -= hh; + } + x = -0.5*LAMBDA; + y = -h; + while (x < 1.5*XMAX) + { + draw_line(x - 10.0, y - 10.0*sqr3, x + 10.0, y + 10.0*sqr3); + draw_line(x - 10.0, y + 10.0*sqr3, x + 10.0, y - 10.0*sqr3); + x += LAMBDA; + } + x = -0.5*LAMBDA; + while (x > 1.5*XMIN) + { + draw_line(x - 10.0, y - 10.0*sqr3, x + 10.0, y + 10.0*sqr3); + draw_line(x - 10.0, y + 10.0*sqr3, x + 10.0, y - 10.0*sqr3); + x -= LAMBDA; + } + break; + } + case (D_HEX_TILES): + { + ntiles = (int)(XMAX/LAMBDA) + 1; + for (i=-ntiles; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, COLOR_PALETTE, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, COLOR_PALETTE, value, 1.0, 1, rgb); +// value = min + 1.0*dy*(double)(j - jmin); +// amp = 0.7*color_amplitude_linear(value, 1.0, 1); +// while (amp > 1.0) amp -= 2.0; +// while (amp < -1.0) amp += 2.0; +// amp_to_rgb(0.5*(1.0 + amp), rgb); + break; + } + case (P_TOTAL_ENERGY_FLUX): + { + value = min + 1.0*dy*(double)(j - jmin); + amp = 0.7*color_amplitude_linear(value, 1.0, 1); + while (amp > 1.0) amp -= 2.0; + while (amp < -1.0) amp += 2.0; + amp_to_rgb(0.5*(1.0 + amp), rgb); + break; + } case (P_PHASE): { value = min + 1.0*dy*(double)(j - jmin); @@ -4042,6 +4342,27 @@ void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int p color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } + case (P_ENERGY_FLUX): + { + value = dy_e*(double)(j - jmin)*100.0/E_SCALE; + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); +// value = min + 1.0*dy*(double)(j - jmin); +// amp = 0.7*color_amplitude_linear(value, 1.0, 1); +// while (amp > 1.0) amp -= 2.0; +// while (amp < -1.0) amp += 2.0; +// amp_to_rgb(0.5*(1.0 + amp), rgb); + break; + } + case (P_TOTAL_ENERGY_FLUX): + { + value = min + 1.0*dy*(double)(j - jmin); + amp = 0.7*color_amplitude_linear(value, 1.0, 1); + while (amp > 1.0) amp -= 2.0; + while (amp < -1.0) amp += 2.0; + amp_to_rgb(0.5*(1.0 + amp), rgb); + break; + } case (P_PHASE): { value = min + 1.0*dy*(double)(j - jmin); @@ -4137,6 +4458,7 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, case (P_LOG_ENERGY): { value = LOG_SHIFT + LOG_SCALE*log(dy_e*(double)(j - jmin)*100.0/E_SCALE); +// printf("value = %.2lg\n", value); // if (value <= 0.0) value = 0.0; color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; @@ -4148,6 +4470,27 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); break; } + case (P_ENERGY_FLUX): + { + value = dy_e*(double)(j - jmin); + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); +// value = min + 1.0*dy*(double)(j - jmin); +// amp = 0.7*color_amplitude_linear(value, 1.0, 1); +// while (amp > 1.0) amp -= 2.0; +// while (amp < -1.0) amp += 2.0; +// amp_to_rgb(0.5*(1.0 + amp), rgb); + break; + } + case (P_TOTAL_ENERGY_FLUX): + { + value = min + 1.0*dy*(double)(j - jmin); + amp = 0.7*color_amplitude_linear(value, 1.0, 1); + while (amp > 1.0) amp -= 2.0; + while (amp < -1.0) amp += 2.0; + amp_to_rgb(0.5*(1.0 + amp), rgb); + break; + } case (P_PHASE): { value = min + 1.0*dy*(double)(j - jmin); diff --git a/sub_wave_3d_rde.c b/sub_wave_3d_rde.c index 3c28e88..8b99b78 100644 --- a/sub_wave_3d_rde.c +++ b/sub_wave_3d_rde.c @@ -82,6 +82,14 @@ void draw_vertex_xyz(double xy[2], double z) glVertex2d(xy_screen[0], xy_screen[1]); } +void draw_vertex_xyz_shift(double xy[2], double z, int shiftx, int shifty) +{ + double xy_screen[2]; + + xyz_to_xy(xy[0] + (double)shiftx*(XMAX - XMIN), xy[1] + (double)shifty*(YMAX - YMIN), z, xy_screen); + glVertex2d(xy_screen[0], xy_screen[1]); +} + void draw_vertex_x_y_z(double x, double y, double z) { double xy_screen[2]; diff --git a/wave_3d.c b/wave_3d.c index 901fdde..a02e9d4 100644 --- a/wave_3d.c +++ b/wave_3d.c @@ -45,44 +45,47 @@ #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 SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ /* General geometrical parameters */ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1000 /* window height */ -#define NX 2000 /* number of grid points on x axis */ -#define NY 2000 /* number of grid points on y axis */ +#define NX 2500 /* number of grid points on x axis */ +#define NY 1250 /* number of grid points on y axis */ // #define NX 2700 /* number of grid points on x axis */ // #define NY 1350 /* number of grid points on y axis */ -// // #define YMIN -1.041666667 -// // #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ - -// #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 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 HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ // #define WINWIDTH 1280 /* window width */ // #define WINHEIGHT 720 /* window height */ -// // // -// #define NX 1500 /* number of grid points on x axis */ -// #define NY 1500 /* 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 2560 /* number of grid points on x axis */ +// #define NY 1440 /* 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 XMIN -1.4 -#define XMAX 1.4 /* x interval */ -#define YMIN -1.4 -#define YMAX 1.4 /* y interval for 9/16 aspect ratio */ -// +// #define XMIN -2.0 +// #define XMAX 2.0 /* x interval */ +// #define YMIN -1.125 +// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ + #define JULIA_SCALE 0.8 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 7 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 57 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */ @@ -98,7 +101,7 @@ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.7 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.25 /* parameter controlling the dimensions of domain */ #define MU 0.0 /* parameter controlling the dimensions of domain */ #define NPOLY 3 /* number of sides of polygon */ #define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ @@ -137,10 +140,10 @@ #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 COURANT 0.06 /* Courant number */ -#define COURANTB 0.15 /* Courant number in medium B */ +#define COURANT 0.05 /* Courant number */ +#define COURANTB 0.1 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMAB 0.015 /* damping factor in wave equation */ +#define GAMMAB 0.0 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ @@ -151,8 +154,9 @@ /* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ -#define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */ -#define OSCILLATING_SOURCE_PERIOD 30 /* period of oscillating source */ +#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ +#define OSCILLATING_SOURCE_PERIOD 100 /* period of oscillating source */ +#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ /* Boundary conditions, see list in global_pdes.c */ @@ -163,9 +167,9 @@ /* Parameters for length and speed of simulation */ -// #define NSTEPS 200 /* number of frames of movie */ -#define NSTEPS 1500 /* number of frames of movie */ -#define NVID 10 /* number of iterations between images displayed on screen */ +#define NSTEPS 1400 /* number of frames of movie */ +// #define NSTEPS 250 /* number of frames of movie */ +#define NVID 5 /* number of iterations between images displayed on screen */ #define NSEG 1000 /* number of segments of boundary */ #define INITIAL_TIME 0 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ @@ -182,12 +186,9 @@ /* Parameters of initial condition */ -// #define INITIAL_AMP 0.75 /* amplitude of initial condition */ -// #define INITIAL_VARIANCE 0.001 /* variance of initial condition */ -// #define INITIAL_WAVELENGTH 0.05 /* wavelength of initial condition */ -#define INITIAL_AMP 0.75 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */ +#define INITIAL_AMP 0.25 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ @@ -195,8 +196,8 @@ #define ZPLOT 103 /* wave height */ #define CPLOT 103 /* color scheme */ -#define ZPLOT_B 105 -#define CPLOT_B 105 /* plot type for second movie */ +#define ZPLOT_B 107 +#define CPLOT_B 107 /* plot type for second movie */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of plot */ #define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ @@ -218,13 +219,13 @@ #define REP_AXO_3D 0 /* linear projection (axonometry) */ #define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */ -#define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */ +#define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ /* Color schemes */ #define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 16 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -237,11 +238,12 @@ #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 ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 100.0 /* scaling factor for energy representation */ -#define LOG_SCALE 0.5 /* scaling factor for energy log representation */ +#define E_SCALE 100.0 /* scaling factor for energy representation */ +#define LOG_SCALE 0.75 /* scaling factor for energy log representation */ #define LOG_SHIFT 0.5 /* shift of colors on log scale */ #define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */ #define LOG_MEAN_ENERGY_SHIFT 1.0 /* additional shift for log of mean energy */ +#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -255,13 +257,11 @@ #define NYMAZE 7 /* height of maze */ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ #define RAND_SHIFT 24 /* seed of random number generator */ -// 0, 9, 18, 20, 22 -// #define RAND_SHIFT 58 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 5.0 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE 3.5 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 3.0 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ @@ -286,7 +286,7 @@ double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "lig double observer[3] = {8.0, 8.0, 8.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.35 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define Z_SCALING_FACTOR 0.3 /* overall scaling factor of z axis for REP_PROJ_3D representation */ #define XY_SCALING_FACTOR 2.4 /* 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 */ @@ -1028,7 +1028,7 @@ void animation() // init_circular_wave_mod(polyline[85].x, polyline[85].y, phi, psi, xy_in); - init_circular_wave_mod(-0.85, 0.0, phi, psi, xy_in); + init_circular_wave_mod(0.0, 0.0, phi, psi, xy_in); // init_wave_flat_mod(phi, psi, xy_in); // add_circular_wave_mod(1.0, 1.0, 0.0, phi, psi, xy_in); @@ -1108,15 +1108,16 @@ void animation() // add_circular_wave_mod(1.0, -1.0, 0.0, phi, psi, xy_in); // add_circular_wave(1.0, -1.5*LAMBDA, 0.0, phi, psi, xy_in); // add_circular_wave(-1.0, 0.6*cos((double)(period)*DPI/3.0), 0.6*sin((double)(period)*DPI/3.0), phi, psi, xy_in); - yshift = (double)period*a + (double)(period*period)*b; - add_circular_wave_mod(sign, -1.5 + yshift, 0.0, phi, psi, xy_in); +// yshift = (double)period*a + (double)(period*period)*b; + if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign; + add_circular_wave_mod(sign, 0.0, 0.0, phi, psi, xy_in); // speed = (a + 2.0*(double)(period)*b)/((double)(NVID)); // speed = 0.55*(a + 2.0*(double)(period)*b)/((double)(NVID*OSCILLATING_SOURCE_PERIOD)); - speed = (a + 2.0*(double)(period)*b)/((double)(3*NVID*OSCILLATING_SOURCE_PERIOD)); - printf("v = %.3lg, c = %.3lg\n", speed, c); - speed = speed/c; - sign = -sign; - period++; +// speed = (a + 2.0*(double)(period)*b)/((double)(3*NVID*OSCILLATING_SOURCE_PERIOD)); +// printf("v = %.3lg, c = %.3lg\n", speed, c); +// speed = speed/c; +// sign = -sign; +// period++; } if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0); diff --git a/wave_billiard.c b/wave_billiard.c index 975a8cd..0ce855e 100644 --- a/wave_billiard.c +++ b/wave_billiard.c @@ -45,6 +45,7 @@ #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 SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ /* General geometrical parameters */ @@ -81,7 +82,7 @@ /* Choice of the billiard table */ -#define B_DOMAIN 53 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 57 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 1 /* pattern of circles or polygons, see list in global_pdes.c */ @@ -93,8 +94,8 @@ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ -#define MU 0.3 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.25 /* parameter controlling the dimensions of domain */ +#define MU 0.0 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 6 /* depth of computation of Menger gasket */ @@ -130,10 +131,10 @@ #define AMPLITUDE 0.8 /* amplitude of periodic excitation */ #define ACHIRP 0.25 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.1 /* Courant number */ -#define COURANTB 0.05 /* Courant number in medium B */ +#define COURANT 0.05 /* Courant number */ +#define COURANTB 0.1 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMAB 5.0e-3 /* damping factor in wave equation */ +#define GAMMAB 0.0 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ #define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ #define KAPPA 0.0 /* "elasticity" term enforcing oscillations */ @@ -145,7 +146,8 @@ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ #define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ -#define OSCILLATING_SOURCE_PERIOD 50 /* period of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 50 /* period of oscillating source */ +#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ /* Boundary conditions, see list in global_pdes.c */ @@ -153,8 +155,8 @@ /* Parameters for length and speed of simulation */ -// #define NSTEPS 1000 /* number of frames of movie */ -#define NSTEPS 2400 /* number of frames of movie */ +#define NSTEPS 1700 /* number of frames of movie */ +// #define NSTEPS 3500 /* number of frames of movie */ #define NVID 12 /* number of iterations between images displayed on screen */ #define NSEG 1000 /* number of segments of boundary */ #define INITIAL_TIME 0 /* time after which to start saving frames */ @@ -172,8 +174,7 @@ /* Parameters of initial condition */ -// #define INITIAL_AMP 0.75 /* amplitude of initial condition */ -#define INITIAL_AMP 1.0 /* amplitude of initial condition */ +#define INITIAL_AMP 0.35 /* amplitude of initial condition */ #define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ #define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */ @@ -186,8 +187,8 @@ /* Color schemes */ -#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 13 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 18 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -198,9 +199,10 @@ #define PHASE_FACTOR 1.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 ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 200.0 /* scaling factor for energy representation */ -#define LOG_SCALE 1.0 /* scaling factor for energy log representation */ -#define LOG_SHIFT 1.0 /* shift of colors on log scale */ +#define E_SCALE 100.0 /* scaling factor for energy representation */ +#define LOG_SCALE 0.25 /* scaling factor for energy log representation */ +#define LOG_SHIFT 0.0 /* shift of colors on log scale */ +#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -211,8 +213,8 @@ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 1.5 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 1.5 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 7.0 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ @@ -234,6 +236,8 @@ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 10.0 /* max value of wave amplitude */ +#define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX) + #include "global_pdes.c" /* constants and global variables */ #include "sub_maze.c" /* support for generating mazes */ #include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ @@ -568,7 +572,7 @@ 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 = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c; - double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX]; + double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX], *total_flux; short int *xy_in[NX]; int i, j, s, sample_left[2], sample_right[2], period = 0, fade; static int counter = 0; @@ -591,6 +595,8 @@ void animation() color_scale[i] = (double *)malloc(NY*sizeof(double)); } + if (MEAN_FLUX) total_flux = (double *)malloc(4*NX*NY*sizeof(double)); + /* initialise positions and radii of circles */ 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); @@ -627,6 +633,10 @@ void animation() for (i=0; i DPI) *arg -= DPI; + + if ((xy_in[i][j])||(TWOSPEEDS)) + { + *module = velocity*module2(gradientx, gradienty); + *gx = velocity*gradientx; + *gy = velocity*gradienty; + } + else + { + *module = 0.0; + *gx = 0.0; + *gy = 0.0; + } + + // if (xy_in[i][j]) return(E_SCALE*E_SCALE*(velocity*COURANT*module2(gradientx,gradienty))); +// else if (TWOSPEEDS) return(E_SCALE*E_SCALE*(velocity*COURANTB*module2(gradientx,gradienty))); + +} + + + double compute_variance(double *phi[NX], double *psi[NX], short int *xy_in[NX]) /* compute the variance of the field, to adjust color scheme */ { @@ -682,12 +721,12 @@ void draw_wave_highres_diss(int size, double *phi[NX], double *psi[NX], double * } -void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[NX], double *color_scale[NX], short int *xy_in[NX], - double scale, int time, int plot, int palette, int fade, double fade_value) +void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, double *color_scale[NX], + short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value) /* same as draw_wave_e, but with color scheme specification */ { int i, j, k, iplus, iminus, jplus, jminus; - double rgb[3], xy[2], x1, y1, x2, y2, velocity, field_value, energy, gradientx2, gradienty2, r2; + double rgb[3], xy[2], x1, y1, x2, y2, velocity, field_value, energy, gradientx2, gradienty2, r2, arg, mod, flux_factor, gx, gy, mgx, mgy, ffactor; static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX); glBegin(GL_QUADS); @@ -755,6 +794,50 @@ void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[N color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb); break; } + case (P_ENERGY_FLUX): + { + compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); +// color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb); +// flux_factor = tanh(mod*E_SCALE); +// for (k=0; k<3; k++) rgb[k] *= flux_factor; + + color_scheme_asym_palette(COLOR_SCHEME, palette, mod*FLUX_SCALE, scale, time, rgb); + break; + } + case (P_TOTAL_ENERGY_FLUX): + { +// ffactor = 1.0; + compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); + total_flux[2*NX*NY + 2*j*NX + 2*i] += gx; + total_flux[2*NX*NY + 2*j*NX + 2*i + 1] += gy; + total_flux[2*j*NX + 2*i] += total_flux[2*NX*NY + 2*j*NX + 2*i]; + total_flux[2*j*NX + 2*i + 1] += total_flux[2*NX*NY + 2*j*NX + 2*i + 1]; +// total_flux[2*j*NX + 2*i] *= 1.0 + 1.0/(double)(time+1); +// total_flux[2*j*NX + 2*i + 1] *= 1.0 + 1.0/(double)(time+1); +// total_flux[2*j*NX + 2*i] *= ffactor; +// total_flux[2*j*NX + 2*i + 1] *= ffactor; +// total_flux[2*j*NX + 2*i] += gx; +// total_flux[2*j*NX + 2*i + 1] += gy; +// total_flux[2*j*NX + 2*i] *= 1.0/ffactor; +// total_flux[2*j*NX + 2*i + 1] *= 1.0/ffactor; +// total_flux[2*j*NX + 2*i] += 0.1*gx; +// total_flux[2*j*NX + 2*i + 1] += 0.1*gy; + mgx = total_flux[2*j*NX + 2*i]; + mgy = total_flux[2*j*NX + 2*i + 1]; +// mgx = total_flux[2*j*NX + 2*i]/sqrt((double)(time+1)); +// mgy = total_flux[2*j*NX + 2*i + 1]/sqrt((double)(time+1)); +// mgx = total_flux[2*j*NX + 2*i]/(1.0 + 0.1*log((double)(time+2))); +// mgy = total_flux[2*j*NX + 2*i + 1]/(1.0 + 0.1*log((double)(time+2))); + mod = module2(mgx, mgy); + arg = argument(mgx, mgy); + if (arg < 0.0) arg += DPI; + color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb); + flux_factor = tanh(mod*FLUX_SCALE); + for (k=0; k<3; k++) rgb[k] *= flux_factor; + +// color_scheme_asym_palette(COLOR_SCHEME, palette, mod, scale, time, rgb); + break; + } } if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; glColor3f(rgb[0], rgb[1], rgb[2]); @@ -770,11 +853,11 @@ void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[N } -void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value) +void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value) /* same as draw_wave_highres, but with color scheme option */ { int i, j, k, iplus, iminus, jplus, jminus; - double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2; + double rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx2, gradienty2, arg, mod, flux_factor, gx, gy, mgx, mgy; static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX); glBegin(GL_QUADS); @@ -834,6 +917,35 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl color_scheme_palette(COLOR_SCHEME, palette, LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb); break; } + case (P_ENERGY_FLUX): + { + compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); + color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb); + flux_factor = tanh(mod*E_SCALE); + for (k=0; k<3; k++) rgb[k] *= flux_factor; + break; + } + case (P_TOTAL_ENERGY_FLUX): + { + compute_energy_flux(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); + total_flux[2*j*NX + 2*i] *= 0.99; + total_flux[2*j*NX + 2*i + 1] *= 0.99; + total_flux[2*j*NX + 2*i] += gx; + total_flux[2*j*NX + 2*i + 1] += gy; +// mgx = total_flux[2*j*NX + 2*i]/(double)(time+1); +// mgy = total_flux[2*j*NX + 2*i + 1]/(double)(time+1); + mgx = total_flux[2*j*NX + 2*i]; + mgy = total_flux[2*j*NX + 2*i + 1]; +// mgx = total_flux[2*j*NX + 2*i]/log((double)(time+2)); +// mgy = total_flux[2*j*NX + 2*i + 1]/log((double)(time+2)); + mod = module2(mgx, mgy); + arg = argument(mgx, mgy); + if (arg < 0.0) arg += DPI; + color_scheme_palette(C_ONEDIM_LINEAR, palette, arg/DPI, 1.0, 1, rgb); + flux_factor = tanh(mod*E_SCALE); + for (k=0; k<3; k++) rgb[k] *= flux_factor; + break; + } } if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; diff --git a/wave_comparison.c b/wave_comparison.c index f6c1a70..bd0b724 100644 --- a/wave_comparison.c +++ b/wave_comparison.c @@ -44,6 +44,7 @@ #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 SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ #define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */ #define TIME_LAPSE_FACTOR 4 /* factor of time-lapse movie */ @@ -200,6 +201,7 @@ #define E_SCALE 200.0 /* scaling factor for energy representation */ #define LOG_SCALE 1.5 /* scaling factor for energy log representation */ #define LOG_SHIFT 1.0 /* shift of colors on log scale */ +#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ diff --git a/wave_energy.c b/wave_energy.c index dad945f..5bddf6c 100644 --- a/wave_energy.c +++ b/wave_energy.c @@ -43,6 +43,7 @@ #include #define MOVIE 0 /* set to 1 to generate movie */ +#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ @@ -178,6 +179,7 @@ #define E_SCALE 500.0 /* scaling factor for energy representation */ #define LOG_SCALE 1.5 /* scaling factor for energy log representation */ #define LOG_SHIFT 1.0 /* shift of colors on log scale */ +#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */