From 46a381dcf370412b39fc9716bc70a57415e5afde Mon Sep 17 00:00:00 2001 From: Nils Berglund <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Tue, 18 Oct 2022 23:28:20 +0200 Subject: [PATCH] Add files via upload --- drop_billiard.c | 9 + global_3d.c | 13 + global_ljones.c | 34 +- global_particles.c | 1 + global_pdes.c | 12 + global_perc.c | 10 + heat.c | 14 +- lennardjones.c | 404 +++++++++++++++++---- mangrove.c | 12 + particle_billiard.c | 83 ++++- particle_pinball.c | 7 + percolation.c | 113 +++--- rde.c | 308 +++++++++++----- schrodinger.c | 13 + sub_lj.c | 838 ++++++++++++++++++++++++++++++++++++++------ sub_maze.c | 231 ++++++++++++ sub_part_billiard.c | 181 +++++++++- sub_perco.c | 663 ++++++++++++++++++++++++++++++----- sub_rde.c | 350 +++++++++++++++++- sub_wave.c | 397 ++++++++++++++++++++- sub_wave_3d.c | 22 +- sub_wave_3d_rde.c | 24 +- wave_3d.c | 151 ++++---- wave_billiard.c | 112 +++--- wave_comparison.c | 13 + wave_energy.c | 12 + 26 files changed, 3484 insertions(+), 543 deletions(-) create mode 100644 sub_maze.c diff --git a/drop_billiard.c b/drop_billiard.c index 12a2364..830940a 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -122,12 +122,21 @@ #define SLEEP2 100 /* final sleeping time */ #define END_FRAMES 25 /* number of frames at end of movie */ +#define NXMAZE 8 /* width of maze */ +#define NYMAZE 8 /* 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 MAZE_XSHIFT 0.0 /* horizontal shift of maze */ + + + #define PI 3.141592654 #define DPI 6.283185307 #define PID 1.570796327 #include "global_particles.c" +#include "sub_maze.c" #include "sub_part_billiard.c" diff --git a/global_3d.c b/global_3d.c index 07a069e..68de88e 100644 --- a/global_3d.c +++ b/global_3d.c @@ -23,6 +23,7 @@ #define E_RPS 4 /* rock-paper-scissors equation */ #define E_RPSLZ 41 /* rock-paper-scissors-lizard-Spock equation */ #define E_SCHRODINGER 5 /* Schrodinger equation */ +#define E_EULER_INCOMP 6 /* incompressigle Euler equation */ /* Choice of potential */ @@ -32,6 +33,12 @@ #define POT_DOUBLE_COULOMB 4 /* sum of Coulomb potentials located at focal points of ellipse */ #define POT_FERMIONS 5 /* two interacting 1D fermions */ #define POT_FERMIONS_PERIODIC 6 /* two interacting 1D fermions on the circle */ +#define POT_MAZE 7 /* higher potential on walls of a maze */ + +/* Choice of vector potential */ + +#define VPOT_CONSTANT_FIELD 100 /* constant magnetic field */ +#define VPOT_AHARONOV_BOHM 101 /* single flux line for Aharonov-Bohm effect */ /* plot types used by rde */ @@ -56,6 +63,11 @@ #define Z_THETA_RPSLZ 41 /* polar angle */ #define Z_NORM_GRADIENT_RPSLZ 42 /* gradient of polar angle */ +/* for Euler equation */ +#define Z_EULER_VORTICITY 50 /* vorticity of velocity */ +#define Z_EULER_LOG_VORTICITY 51 /* log of vorticity of velocity */ +#define Z_EULER_VORTICITY_ASYM 52 /* vorticity of velocity */ + /* macros to avoid unnecessary computations in 3D plots */ #define COMPUTE_THETA ((cplot == Z_POLAR)||(cplot == Z_NORM_GRADIENT)||(cplot == Z_ANGLE_GRADIENT)||(cplot == Z_NORM_GRADIENT_INTENSITY)||(cplot == Z_VORTICITY)||(cplot == Z_VORTICITY_ABS)) @@ -104,6 +116,7 @@ typedef struct double field_arg; /* argument of field or gradient */ double curl; /* curl of field */ double cos_angle; /* cos of angle between normal vector and direction of light */ + double log_vorticity; /* logarithm of vorticity (for Euler equation) */ double rgb[3]; /* RGB color code */ double *p_zfield[2]; /* pointers to z field (second pointer for option DOUBLE_MOVIE) */ double *p_cfield[2]; /* pointers to color field (second pointer for option DOUBLE_MOVIE) */ diff --git a/global_ljones.c b/global_ljones.c index 03e3d2e..2b0b909 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -14,7 +14,8 @@ #define NMAXCIRCLES 20000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define MAXNEIGH 20 /* max number of neighbours kept in memory */ #define NMAXOBSTACLES 100 /* max number of obstacles */ -#define NMAXSEGMENTS 1000 /* max number of repelling segments */ +#define NMAXSEGMENTS 1000 /* max number of repelling segments */ +#define NMAXGROUPS 50 /* max number of groups of segments */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -60,6 +61,9 @@ #define S_DAM 12 /* segments forming a dam that can break */ #define S_DAM_WITH_HOLE 13 /* segments forming a dam in which a hole can open */ #define S_DAM_WITH_HOLE_AND_RAMP 14 /* segments forming a dam in which a hole can open */ +#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 */ /* particle interaction */ @@ -102,6 +106,12 @@ #define G_INCREASE_RELEASE 1 /* slow increase and instant release */ #define G_INCREASE_DECREASE 2 /* slow increase an decrease */ +/* Rocket shapes */ + +#define RCK_DISC 0 /* disc-shaped rocket */ +#define RCK_RECT 1 /* rectangular rocket */ +#define RCK_RECT_HAT 2 /* rectangular rocket with a hat */ + /* Nozzle shapes */ #define NZ_STRAIGHT 0 /* straight nozzle */ @@ -109,6 +119,7 @@ #define NZ_GLAS 2 /* glas-shaped nozzle */ #define NZ_CONE 3 /* cone-shaped nozzle */ #define NZ_TRUMPET 4 /* trumpet-shaped nozzle */ +#define NZ_BROAD 5 /* broad straight nozzle */ #define NZ_NONE 99 /* no nozzle */ /* Plot types */ @@ -218,6 +229,7 @@ typedef struct typedef struct { double x1, x2, y1, y2; /* extremities of segment */ + double xc, yc; /* mid-point of segment */ double nx, ny; /* normal vector */ double c; /* constant term in cartesian eq nx*x + ny*y = c */ double length; /* length of segment */ @@ -230,6 +242,7 @@ typedef struct double nx0, ny0; /* initial normal vector */ double angle01, angle02; /* initial values of angles in which concave corners repel */ double fx, fy; /* x and y-components of force on segment */ + double torque; /* torque on segment with respect to its center */ short int inactivate; /* set to 1 for segment to become inactive at time SEGMENT_DEACTIVATION_TIME */ } t_segment; @@ -238,6 +251,23 @@ typedef struct double xc, yc; /* center of circle */ } t_tracer; +typedef struct +{ + double xc, yc; /* center of mass of obstacle */ + double angle; /* orientation of obstacle */ + double vx, vy; /* velocity of center of mass */ + double omega; /* angular velocity */ + double mass; /* mass of obstacle */ + double moment_inertia; /* moment of inertia */ +} t_group_segments; -int ncircles, nobstacles, nsegments, counter = 0; +typedef struct +{ + double xc, yc; /* coordinates of centers of mass */ + double vx, vy; /* velocities */ + double omega; /* angular velocity */ +} t_group_data; + + +int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0; diff --git a/global_particles.c b/global_particles.c index 3e011ca..16cf1fe 100644 --- a/global_particles.c +++ b/global_particles.c @@ -78,6 +78,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define P_TOKA_PRIME 6 /* Tokarsky room made of 86 triangles */ #define P_TREE 7 /* pine tree */ #define P_TOKA_NONSELF 8 /* Tokarsky non-self-unilluminable room */ +#define P_MAZE 10 /* maze */ /* Color palettes */ diff --git a/global_pdes.c b/global_pdes.c index 70f83a0..0f5d99d 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -61,6 +61,10 @@ #define D_GROOVE 48 /* groove array supposed to induce polaritons */ #define D_FABRY_PEROT 49 /* Fabry-Perrot cavity (in fact simply a vertical slab) */ #define D_LSHAPE 50 /* L-shaped billiard (surface of genus 2) */ +#define D_WAVEGUIDE 51 /* wave guide */ +#define D_WAVEGUIDE_W 52 /* W-shaped wave guide */ +#define D_MAZE 53 /* maze */ +#define D_MAZE_CLOSED 54 /* closed maze */ #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) */ @@ -192,6 +196,12 @@ typedef struct double posi, posj; /* (i,j) coordinates of vertex */ } t_vertex; +typedef struct +{ + double x1, y1, x2, y2; /* (x,y) coordinates of vertices */ + double posi1, posj1, posi2, posj2; /* (i,j) coordinates of vertices */ +} t_rectangle; + typedef struct { int nneighb; /* number of neighbours to compute Laplacian */ @@ -201,10 +211,12 @@ typedef struct int ncircles = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */ int npolyline = NMAXPOLY; /* actual length of polyline */ +int npolyrect = NMAXPOLY; /* actual number of polyrect */ t_circle circles[NMAXCIRCLES]; /* circular scatterers */ t_polygon polygons[NMAXCIRCLES]; /* polygonal scatterers */ t_vertex polyline[NMAXPOLY]; /* vertices of polygonal line */ +t_rectangle polyrect[NMAXPOLY]; /* vertices of rectangles */ /* the same for comparisons between different domains */ int ncircles_b = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */ diff --git a/global_perc.c b/global_perc.c index 33db13c..e2f0e7e 100644 --- a/global_perc.c +++ b/global_perc.c @@ -16,6 +16,8 @@ #define BC_TRIANGLE_SITE_DIRICHLET 20 /* triangular lattice, site percolation, Dirichlet b.c. */ #define BC_POISSON_DISC 50 /* Poisson disc (blue noise) lattice */ +#define BC_CUBIC_DIRICHLET 100 /* cubic lattice */ + /* Plot types */ #define PLOT_SQUARES 0 /* plot squares */ @@ -25,6 +27,14 @@ #define PLOT_HEX_BONDS 4 /* plot edges of hexagonal lattice */ #define PLOT_POISSON_DISC 5 /* plot Poisson disc process */ +#define PLOT_CUBES 10 /* plot cubes */ + +/* 3D representation */ + +#define REP_AXO_3D 0 /* linear projection (axonometry) */ +#define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */ + + /* Color schemes */ #define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */ diff --git a/heat.c b/heat.c index 935c5b5..e13107e 100644 --- a/heat.c +++ b/heat.c @@ -35,7 +35,7 @@ #include /* Sam Leffler's libtiff library. */ #include -#define MOVIE 1 /* set to 1 to generate movie */ +#define MOVIE 0 /* set to 1 to generate movie */ /* General geometrical parameters */ @@ -188,7 +188,19 @@ #define AMPLITUDE 0.8 /* amplitude of periodic excitation */ /* end of constants only used by wave_billiard and wave_3d */ +/* for compatibility with sub_wave and sub_maze */ +#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 MAZE_XSHIFT 0.0 /* horizontal shift of maze */ +#define ADD_POTENTIAL 0 +#define POT_MAZE 7 +#define POTENTIAL 0 +/* end of constants only used by sub_wave and sub_maze */ + #include "global_pdes.c" +#include "sub_maze.c" #include "sub_wave.c" diff --git a/lennardjones.c b/lennardjones.c index 21c36e3..e2d1a36 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -49,36 +49,48 @@ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ -#define XMIN -2.0 -#define XMAX 2.0 /* x interval */ -#define YMIN -1.125 -#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +// #define XMIN -2.3 +// #define XMAX 3.7 /* x interval */ +// #define YMIN -1.6875 +// #define YMAX 1.6875 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -2.0 -#define INITXMAX -0.04 /* x interval for initial condition */ -#define INITYMIN -1.125 -#define INITYMAX 0.65 /* y interval for initial condition */ +#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 BCXMIN -2.05 -#define BCXMAX 2.0 /* x interval for boundary condition */ -#define BCYMIN -1.125 -#define BCYMAX 1.25 /* y interval for boundary 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 OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ -#define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_ljones.c */ +#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_ljones.c */ #define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */ #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 14 /* pattern of repelling segments, see list in global_ljones.c */ +#define SEGMENT_PATTERN 102 /* 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 */ -#define NOZZLE_SHAPE_B 2 /* shape of nozzle for second rocket, see list in global_ljones.c */ +#define NOZZLE_SHAPE_B 0 /* shape of nozzle for second rocket, see list in global_ljones.c */ #define TWO_TYPES 0 /* set to 1 to have two types of particles */ -#define TPYE_PROPORTION 0.7 /* proportion of particles of first type */ +#define TPYE_PROPORTION 0.5 /* proportion of particles of first type */ #define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ #define CENTER_PX 0 /* set to 1 to center horizontal momentum */ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ @@ -91,12 +103,15 @@ #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.5 /* minimal distance in Poisson disc process, controls density of particles */ -#define PDISC_CANDIDATES 50 /* number of candidates in construction of Poisson disc process */ +#define PDISC_DISTANCE 2.7 /* 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.0075 /* parameter controlling radius of particles */ +// #define MU 0.02 /* parameter controlling radius of particles */ +// #define MU 0.015 /* parameter controlling radius of particles */ +#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 */ @@ -119,11 +134,13 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 3500 /* number of frames of movie */ -#define NVID 60 /* number of iterations between images displayed on screen */ +#define NSTEPS 3300 /* number of frames of movie */ +// #define NSTEPS 2000 /* 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 0 /* time after which to start saving frames */ -#define OBSTACLE_INITIAL_TIME 10 /* time after which to start moving obstacle */ +#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 BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ #define CONTAINER_WIDTH 4 /* width of container boundary */ @@ -137,16 +154,17 @@ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 0 +#define BOUNDARY_COND 20 /* Plot type, see list in global_ljones.c */ -#define PLOT 11 -#define PLOT_B 0 /* plot type for second movie */ +#define PLOT 0 +#define PLOT_B 8 /* plot type for second movie */ -#define DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */ +#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 COLOR_SEG_GROUPS 1 /* set to 1 to collor segment groups differently */ /* Color schemes */ @@ -173,36 +191,36 @@ #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 1.0e3 /* energy of particle with hottest color */ -#define HUE_TYPE0 280.0 /* hue of particles of type 0 */ -#define HUE_TYPE1 70.0 /* hue of particles of type 1 */ -#define HUE_TYPE2 70.0 /* hue of particles of type 2 */ +#define PARTICLE_EMAX 2.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 */ #define HUE_TYPE3 210.0 /* hue of particles of type 3 */ #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 3.5 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST 4.5 /* 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 10.0 /* damping coefficient of particles */ -#define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ +#define DAMPING 5.0 /* damping coefficient of particles */ +#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 */ #define PARTICLE_INERTIA_MOMENT_B 0.02 /* moment of inertia of second type of particle */ #define V_INITIAL 0.0 /* initial velocity range */ #define OMEGA_INITIAL 10.0 /* initial angular velocity range */ -#define THERMOSTAT 0 /* set to 1 to switch on thermostat */ +#define THERMOSTAT 1 /* 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 */ +#define BETA 0.02 /* initial inverse temperature */ #define MU_XI 0.01 /* friction constant in thermostat */ -#define KSPRING_BOUNDARY 1.0e11 /* confining harmonic potential outside simulation region */ +#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.0 /* radius in which to count neighbours */ -#define GRAVITY 2000.0 /* gravity acting on all particles */ -#define GRAVITY_X 0.0 /* gravity acting on all particles */ +#define NBH_DIST_FACTOR 7.5 /* radius in which to count neighbours */ +#define GRAVITY 15.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 */ #define GRAVITY_FACTOR 100.0 /* factor by which to increase gravity */ @@ -222,13 +240,13 @@ #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 0 /* set to 1 to increase BETA during simulation */ -#define BETA_FACTOR 0.01 /* factor by which to change BETA during simulation */ +#define INCREASE_BETA 1 /* 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 */ -#define MIDDLE_CONSTANT_PHASE 370 /* final phase in which temperature is constant */ -#define FINAL_DECREASE_PHASE 350 /* final phase in which temperature decreases */ -#define FINAL_CONSTANT_PHASE 1180 /* final phase in which temperature is constant */ +#define MIDDLE_CONSTANT_PHASE 1600 /* final phase in which temperature is constant */ +#define FINAL_DECREASE_PHASE 1500 /* final phase in which temperature decreases */ +#define FINAL_CONSTANT_PHASE -1 /* final phase in which temperature is constant */ #define DECREASE_CONTAINER_SIZE 0 /* set to 1 to decrease size of container */ #define SYMMETRIC_DECREASE 0 /* set tp 1 to decrease container symmetrically */ @@ -249,7 +267,7 @@ #define N_P_AVERAGE 100 /* size of pressure averaging window */ #define N_T_AVERAGE 50 /* size of temperature averaging window */ #define MAX_PRESSURE 3.0e10 /* pressure shown in "hottest" color */ -#define PARTIAL_THERMO_COUPLING 1 /* set to 1 to couple only particles to the right of obstacle to thermostat */ +#define PARTIAL_THERMO_COUPLING 0 /* set to 1 to couple only particles to the right of obstacle to thermostat */ #define PARTIAL_THERMO_REGION 1 /* region for partial thermostat coupling (see list in global_ljones.c) */ #define PARTIAL_THERMO_SHIFT 0.2 /* distance from obstacle at the right of which particles are coupled to thermostat */ #define PARTIAL_THERMO_WIDTH 0.5 /* vertical size of partial thermostat coupling */ @@ -283,21 +301,33 @@ #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 0 /* set to 1 to print velocity of moving segments */ +#define PRINT_SEGMENTS_SPEEDS 1 /* 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 500 /* time at which to deactivate last segment */ -#define RELEASE_ROCKET_AT_DEACTIVATION 0 /* set to 1 to limit segments velocity before segment release */ -#define SEGMENTS_X0 0.0 /* initial position of segments */ -#define SEGMENTS_Y0 1.5 /* initial position of segments */ +#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_VX0 0.0 /* initial velocity of segments */ -#define SEGMENTS_VY0 -4.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_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_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_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 */ -#define POSITION_Y_DEPENDENCE 0 /* set to 1 for the separation between particles to be vertical */ +#define POSITION_Y_DEPENDENCE 0 /* set to 1 for the separation between particles to be horizontal */ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ #define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */ @@ -305,7 +335,10 @@ #define REACTION_PROB 0.0045 /* probability controlling reaction term */ #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ -#define PRINT_LEFT 1 /* set to 1 to print certain parameters at the top left instead of right */ +#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 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) */ @@ -315,15 +348,21 @@ #define WALL_FRICTION 0.0 /* friction on wall for BC_RECTANGLE_WALL b.c. */ #define WALL_WIDTH 0.1 /* width of wall for BC_RECTANGLE_WALL b.c. */ #define WALL_VMAX 100.0 /* max speed of wall */ -#define WALL_TIME 500 /* time during which to keep wall */ +#define WALL_TIME 0 /* time during which to keep wall */ + +#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 200 /* seed of random number generator */ +#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ #define FLOOR_FORCE 1 /* set to 1 to limit force on particle to FMAX */ #define FMAX 1.0e12 /* maximal force */ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 160 /* size of hashgrid in x direction */ -#define HASHY 80 /* size of hashgrid in y direction */ +#define HASHX 120 /* size of hashgrid in x direction */ +#define HASHY 450 /* 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 */ @@ -343,6 +382,8 @@ double vylid = 0.0; /* y speed coordinate of lid (for BC_RECTANGLE_LID b.c double xwall = 0.0; /* x coordinate of wall (for BC_RECTANGLE_WALL b.c.) */ double vxwall = 0.0; /* x speed of wall (for BC_RECTANGLE_WALL b.c.) */ double angular_speed = 0.0; /* angular speed of rotating segments */ +double xtrack = 0.0; /* traking coordinate */ +double ytrack = 0.0; /* traking coordinate */ double xsegments[2] = {SEGMENTS_X0, -SEGMENTS_X0}; /* x coordinate of segments (for option MOVE_BOUNDARY) */ double ysegments[2] = {SEGMENTS_Y0, SEGMENTS_Y0}; /* y coordinate of segments (for option MOVE_BOUNDARY) */ double vxsegments[2] = {SEGMENTS_VX0, SEGMENTS_VX0}; /* vx coordinate of segments (for option MOVE_BOUNDARY) */ @@ -352,6 +393,7 @@ int thermostat_on = 1; /* thermostat switch used when VARY_THERMOSTAT is on * #define THERMOSTAT_ON ((THERMOSTAT)&&((!VARY_THERMOSTAT)||(thermostat_on))) #include "global_ljones.c" +#include "sub_maze.c" #include "sub_lj.c" #include "sub_hashgrid.c" @@ -836,32 +878,195 @@ void evolve_segments(t_segment segment[NMAXSEGMENTS], int time) } +void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_segments segment_group[NMAXGROUPS]) +/* new version of evolve_segments that takes the group structure into account */ +{ + double fx[NMAXGROUPS], fy[NMAXGROUPS], torque[NMAXGROUPS], dx[NMAXGROUPS], dy[NMAXGROUPS], dalpha[NMAXGROUPS]; + double x, y, dx0, dy0, padding, proj, distance, f, xx[2], yy[2], xmean = 0.0, ymean = 0.0; + int i, j, k, group = 0; + static double maxdepth, saturation_depth; + + maxdepth = 0.5*GROUP_WIDTH; + saturation_depth = 0.1*GROUP_WIDTH; + + for (group=0; group 0)) + { + group = segment[i].group; + + fx[group] += segment[i].fx; + fy[group] += segment[i].fy; + torque[group] += segment[i].torque; + + dx0 = segment[i].xc - segment_group[group].xc; + dy0 = segment[i].yc - segment_group[group].yc; + torque[group] += dx0*segment[i].fy - dy0*segment[i].fx; + + if (BOUNDARY_COND == BC_SCREEN) /* add force from simulation boundary */ + { + x = 0.5*(segment[i].x1 + segment[i].x2); + y = 0.5*(segment[i].y1 + segment[i].y2); + if (x < XMIN + padding) fx[group] += KSPRING_BOUNDARY*(XMIN + padding - x); + else if (x > XMAX - padding) fx[group] -= KSPRING_BOUNDARY*(x - XMAX + padding); + if (y < YMIN + padding) fy[group] += KSPRING_BOUNDARY*(YMIN + padding - y); + else if (y > YMAX - padding) fy[group] -= KSPRING_BOUNDARY*(y - YMAX + padding); + } + else if (BOUNDARY_COND == BC_REFLECT_ABS) /* add force from simulation boundary */ + { + y = 0.5*(segment[i].y1 + segment[i].y2); + if (y < YMIN) fy[group] += KSPRING_BOUNDARY*(YMIN - y); + } + + /* repulsion between different groups */ + if (GROUP_REPULSION) for (j=0; j 0.0)&&(proj < 1.0)) + { + distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; + if ((distance > -maxdepth)&&(distance < 0.0)) + { + if (distance < -saturation_depth) distance = -saturation_depth; + f = KSPRING_GROUPS*(-distance); + segment[j].fx += f*segment[i].nx; + segment[j].fy += f*segment[i].ny; + segment[j].torque += (x - segment[i].xc)*f*segment[i].ny - (y - segment[i].yc)*f*segment[i].nx; + + fx[group] -= f*segment[i].nx; + fy[group] -= f*segment[i].ny; + torque[group] -= (x - segment[i].xc)*f*segment[i].ny - (y - segment[i].yc)*f*segment[i].nx; + } + } + } + } + } + + if (GROUP_G_REPEL) for (i=0; i FMAX) fx[group] = FMAX; + else if (fx[group] < -FMAX) fx[group] = -FMAX; + if (fy[group] > FMAX) fy[group] = FMAX; + else if (fy[group] < -FMAX) fy[group] = -FMAX; + } + + for (group=1; group 0)) + { + group = segment[i].group; + + translate_one_segment(segment, i, dx[group], dy[group]); + rotate_one_segment(segment, i, dalpha[group], segment_group[group].xc, segment_group[group].yc); + } + + if (TRACK_SEGMENT_GROUPS) + { + /* compute mean position */ + for (group=1; group ytrack) ytrack = ymean; + if (xmean > XMAX - TRACK_X_PADDING) + xtrack = xmean - XMAX + TRACK_X_PADDING; + else if (xmean < XMIN + TRACK_X_PADDING) + xtrack = xmean - XMIN - TRACK_X_PADDING; + } +} + + void animation() { double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, a, b, length, fx, fy, force[2], totalenergy = 0.0, krepel = KREPEL, pos[2], prop, vx, beta = BETA, xi = 0.0, xmincontainer = BCXMIN, xmaxcontainer = BCXMAX, torque, torque_ij, - fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY; - double *qx, *qy, *px, *py, *qangle, *pangle, *pressure; + fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio; + double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds; int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, min_nb, max_nb, close, wrapx = 0, wrapy = 0, nactive = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, - tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0; + tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0, + group, gshift; static int imin, imax; static short int first = 1; t_particle *particle; t_obstacle *obstacle; t_segment *segment; + t_group_segments *segment_group; t_tracer *trajectory; + t_group_data *group_speeds; t_hashgrid *hashgrid; char message[100]; particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* circular obstacles */ - if (ADD_FIXED_SEGMENTS) segment = (t_segment *)malloc(NMAXSEGMENTS*sizeof(t_segment)); /* linear obstacles */ + if (ADD_FIXED_SEGMENTS) + { + segment = (t_segment *)malloc(NMAXSEGMENTS*sizeof(t_segment)); /* linear obstacles */ + segment_group = (t_group_segments *)malloc(NMAXGROUPS*sizeof(t_group_segments)); + } if (TRACER_PARTICLE) trajectory = (t_tracer *)malloc(TRAJECTORY_LENGTH*N_TRACER_PARTICLES*sizeof(t_tracer)); - + hashgrid = (t_hashgrid *)malloc(HASHX*HASHY*sizeof(t_hashgrid)); /* hashgrid */ qx = (double *)malloc(NMAXCIRCLES*sizeof(double)); @@ -872,17 +1077,27 @@ void animation() pangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); pressure = (double *)malloc(N_PRESSURES*sizeof(double)); + /* initialise positions and radii of circles */ init_particle_config(particle); init_hashgrid(hashgrid); xshift = OBSTACLE_XMIN; + speed_ratio = (double)(25*NVID)*DT_PARTICLE; if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle); if (ADD_FIXED_SEGMENTS) init_segment_config(segment); + if (MOVE_SEGMENT_GROUPS) + { + for (i=0; i OBSTACLE_INITIAL_TIME)) evolve_segments(segment, i); + + if ((MOVE_SEGMENT_GROUPS)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segment_groups(segment, i, segment_group); } /* end of for (n=0; nN_T_AVERAGE)) @@ -1091,13 +1332,19 @@ void animation() print_entropy(entropy); } + if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i); + if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i); + if (PRINT_OMEGA) print_omega(angular_speed); else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); - else if (PRINT_SEGMENTS_SPEEDS) print_segments_speeds(vxsegments, vysegments); + else if (PRINT_SEGMENTS_SPEEDS) + { + if (MOVE_BOUNDARY) print_segments_speeds(vxsegments, vysegments); + else print_segment_group_speeds(segment_group); + } glutSwapBuffers(); - if (MOVIE) { if (i >= INITIAL_TIME) @@ -1129,11 +1376,14 @@ void animation() draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity); + if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i); + if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i); if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles); if (PRINT_OMEGA) print_omega(angular_speed); else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); - else if (PRINT_SEGMENTS_SPEEDS) print_segments_speeds(vxsegments, vysegments); + else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group); +// print_segments_speeds(vxsegments, vysegments); glutSwapBuffers(); save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); counter++; @@ -1160,11 +1410,14 @@ void animation() draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); print_parameters(beta, mean_energy, krepel, xmaxcontainer - xmincontainer, fboundary/(double)(ncircles*NVID), PRINT_LEFT, pressure, gravity); + if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i); + if (PLOT_TRAJECTORIES) draw_trajectory_plot(group_speeds, i); if (BOUNDARY_COND == BC_EHRENFEST) print_ehrenfest_parameters(particle, pleft, pright); else if (PRINT_PARTICLE_NUMBER) print_particle_number(ncircles); if (PRINT_OMEGA) print_omega(angular_speed); else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); - else if (PRINT_SEGMENTS_SPEEDS) print_segments_speeds(vxsegments, vysegments); + else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group); +// print_segments_speeds(vxsegments, vysegments); glutSwapBuffers(); } for (i=0; i xmax) nright++; + } + + hsl_to_rgb(0.0, 0.0, 0.0, rgb); + + erase_area(xl, yl - 0.03, 0.5, 0.12, rgb); + 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); + write_text(xl, yl, message); + if (nright > 1) sprintf(message, "%i particles", nright); + else sprintf(message, "%i particle", nright); + write_text(xr, yr, message); +} + void print_collision_number(int ncollisions, double x, double y) { char message[50]; @@ -730,7 +770,7 @@ void animation() // init_drop_config(-0.5, 0.0, 0.2, 0.4, configs); - init_drop_config(0.0, 0.0, 0.0, DPI, configs); + init_drop_config(-1.3, 0.2, -0.25*PI, 0.25*PI, configs); // init_drop_config(-1.3, -0.1, 0.0, DPI, configs); // init_drop_config(1.4, 0.1, 0.0, DPI, configs); @@ -746,11 +786,13 @@ void animation() // init_line_config(-1.0, -0.3, -1.0, 0.3, 0.0, configs); // init_line_config(-0.7, -0.45, -0.7, 0.45, 0.0, configs); // init_line_config(-1.5, 0.1, -0.1, 1.0, -0.5*PID, configs); - + if (!SHOWTRAILS) blank(); glColor3f(0.0, 0.0, 0.0); 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); else if (PRINT_COLLISION_NUMBER) print_collision_number(ncollisions, XMIN + 0.1, YMIN + 0.1); glutSwapBuffers(); @@ -787,6 +829,8 @@ void animation() color[i] = (i*NCOLORS)/NPART; newcolor[i] = (i*NCOLORS)/NPART; } + + if (TEST_INITIAL_COND) nparticles = test_initial_condition(configs, active, newcolor); sleep(SLEEP1); @@ -805,6 +849,9 @@ void animation() // draw_config(newcolor, configs, active); 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, 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); for (j=0; j #include -#define MOVIE 1 /* set to 1 to generate movie */ +#define MOVIE 0 /* set to 1 to generate movie */ /* General geometrical parameters */ -// #define WINWIDTH 1920 /* window width */ -// #define WINHEIGHT 1000 /* window height */ +#define WINWIDTH 1920 /* window width */ +#define WINHEIGHT 1000 /* window height */ // #define NX 1920 /* number of grid points on x axis */ // #define NY 992 /* number of grid points on y axis */ -// -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -// #define YMIN -1.041666667 -// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ - -#define HIGHRES 0 /* 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 1280 /* number of grid points on x axis */ -#define NY 720 /* number of grid points on y axis */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ -#define YMIN -1.125 -#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +#define YMIN -1.041666667 +#define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ + +#define HIGHRES 0 /* 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 256 /* number of grid points on x axis */ +#define NY 256 /* number of grid points on y axis */ +#define NZ 256 /* number of grid points on z axis, for 3D percolation */ + +// #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 */ /* Boundary conditions, see list in global_pdes.c */ -#define LATTICE 2 +#define LATTICE 100 -#define FLOOD_LEFT_BOUNDARY 0 /* set to 1 to flood cells on left boundary */ -#define FIND_ALL_CLUSTERS 1 /* set to 1 to find all open clusters */ +#define FLOOD_LEFT_BOUNDARY 0 /* set to 1 to flood cells on left boundary */ +#define FLOOD_BOTTOM_BOUNDARY 0 /* set to 1 to flood cells on bottom boundary */ +#define FIND_ALL_CLUSTERS 1 /* set to 1 to find all open clusters */ +#define PLOT_ONLY_FLOODED_CELLS 0 /* set to 1 to plot only flooded cells */ #define PLOT_CLUSTER_SIZE 0 /* set to 1 to add a plot for the size of the percolation cluster */ #define PLOT_CLUSTER_NUMBER 0 /* set to 1 to add a graph of the number of clusters */ #define PLOT_CLUSTER_HISTOGRAM 1 /* set to 1 to add a histogram of the number of clusters */ -#define PRINT_LARGEST_CLUSTER_SIZE 1 /* set to 1 to print size of largest cluster */ +#define PRINT_LARGEST_CLUSTER_SIZE 0 /* set to 1 to print size of largest cluster */ +#define HISTO_X_LOG_SCALE 1 /* set to 1 to use a log scale on cluster sizes */ + +#define P_SCHEDULE_POWER 4 /* power controlling slowing down near pc - 2 is standard, higher values mean slower passage */ #define MAX_CLUSTER_NUMBER 6 /* vertical scale of the cluster number plot */ #define HISTO_BINS 30 /* number of bins in histogram */ #define NSTEPS 100 /* number of frames of movie */ -// #define NSTEPS 700 /* number of frames of movie */ -// #define NSTEPS 830 /* number of frames of movie */ +// #define NSTEPS 760 /* number of frames of movie */ #define PAUSE 200 /* number of frames after which to pause */ #define PSLEEP 2 /* sleep time during pause */ @@ -95,11 +100,13 @@ /* Color schemes */ -#define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ #define COLOR_CLUSTERS_BY_SIZE 1 /* set to 1 to link cluster color to their size */ +#define COLOR_CELLS_BY_XCOORD 0 /* set to 1 to color cells according to their x-coordinate */ +#define COLOR_CELLS_BY_ZCOORD 0 /* set to 1 to color cells according to their z-coordinate */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ @@ -107,10 +114,11 @@ #define HUE_CLOSED 350.0 /* color hue of closed cells */ #define HUE_OPEN 200.0 /* color hue of open (dry) cells */ -#define HUE_FLOODED 45.0 /* color hue of open flooded cells */ -#define HUE_GRAPH 250.0 /* color hue in graph of cluster size */ +#define HUE_FLOODED 300.0 /* color hue of open flooded cells */ +#define HUE_GRAPH_SIZE 250.0 /* color hue in graph of cluster size */ +#define HUE_GRAPH_CLUSTERS 150.0 /* color hue in graph of cluster size */ -#define CLUSTER_HUEMIN 60.0 /* minimal color hue of clusters */ +#define CLUSTER_HUEMIN 10.0 /* minimal color hue of clusters */ #define CLUSTER_HUEMAX 300.0 /* maximal color hue of clusters */ #define N_CLUSTER_COLORS 20 /* number of different colors of clusters */ @@ -121,6 +129,24 @@ #define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ +/* parameters of 3D representation */ +double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_3D representation */ +double v_3d[2] = {-0.75, -0.45}; +double w_3d[2] = {0.0, 0.015}; +double light[3] = {0.816496581, 0.33333333, 0.4714045}; /* vector of "light" direction for P_3D_ANGLE color scheme */ +double observer[3] = {8.0, 11.0, 10.0}; /* location of observer for REP_PROJ_3D representation */ +int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ + +#define REPRESENTATION_3D 1 /* choice of 3D representation */ +#define Z_SCALING_FACTOR 1.0 /* 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.4 /* overall x shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */ + +#define ROTATE_VIEW 0 /* set to 1 to rotate viewpoint */ +#define ROTATE_ANGLE 360.0 /* total angle of viewpoint rotation */ + /* debugging options */ #define VERBOSE 0 /* set to 1 to print more messages in shell */ #define DEBUG 0 /* set to 1 for some debugging features */ @@ -130,24 +156,25 @@ #define ADD_PLOT ((PLOT_CLUSTER_SIZE)||(PLOT_CLUSTER_NUMBER)||(PLOT_CLUSTER_HISTOGRAM)) #define FIND_CLUSTER_SIZES ((COLOR_CLUSTERS_BY_SIZE)||(PLOT_CLUSTER_HISTOGRAM)) +#define PLOT_3D (LATTICE == BC_CUBIC_DIRICHLET) #include "global_perc.c" /* constants and global variables */ +#include "sub_perco_3d.c" /* support for 3D graphics */ #include "sub_perco.c" - void animation(int size) { - int i, j, k, s, nx, ny, nmaxcells, maxsize, nopen, nflooded, nstack, nclusters, maxclustersize = 0, maxclusterlabel; + int i, j, k, s, nx, ny, nz, nmaxcells, maxsize, nopen, nflooded, nstack, nclusters, maxclustersize = 0, maxclusterlabel; int *plot_cluster_number, *cluster_sizes; int ncells; double p, *plot_cluster_size; t_perco *cell; t_perco **pstack; - compute_nxny(size, &nx, &ny); + compute_nxnynz(size, &nx, &ny, &nz); - nmaxcells = cell_number(NX, NY); + nmaxcells = cell_number(NX, NY, NZ); cell = (t_perco *)malloc(nmaxcells*sizeof(t_perco)); if (PLOT_CLUSTER_SIZE) plot_cluster_size = (double *)malloc(NSTEPS*sizeof(double)); @@ -155,8 +182,9 @@ void animation(int size) // if (FIND_CLUSTER_SIZES) cluster_sizes = (int *)malloc(2*nmaxcells*sizeof(int)); - ncells = init_cell_lattice(cell, nx, ny); - printf("nx = %i, ny = %i, ncells = %i, maxcells = %i\n", nx, ny, ncells, nmaxcells); + ncells = init_cell_lattice(cell, nx, ny, nz); + + printf("nx = %i, ny = %i, nz = %i, ncells = %i, maxcells = %i\n", nx, ny, nz, ncells, nmaxcells); pstack = (t_perco* *)malloc(ncells*sizeof(struct t_perco *)); @@ -169,14 +197,21 @@ void animation(int size) p = p_schedule(i); printf("\ni = %i, p = %.4lg\n", i, p); + if (ROTATE_VIEW) + { + viewpoint_schedule(i); + reset_view = 1; + } + init_cell_state(cell, p, ncells, (i == 0)); - if (FLOOD_LEFT_BOUNDARY) nstack = init_flooded_cells(cell, ncells, nx, ny, pstack); + if (FLOOD_LEFT_BOUNDARY) nstack = init_flooded_cells(cell, ncells, nx, ny, nz, 0, pstack); + if (FLOOD_BOTTOM_BOUNDARY) nstack = init_flooded_cells(cell, ncells, nx, ny, nz, 1, pstack); nopen = count_open_cells(cell, ncells); printf("Flooded cells, %i open cells, nstack = %i\n", nopen, nstack); - if (FLOOD_LEFT_BOUNDARY) + if ((FLOOD_LEFT_BOUNDARY)||(FLOOD_BOTTOM_BOUNDARY)) { nflooded = find_percolation_cluster(cell, ncells, pstack, nstack); printf("Found percolation cluster with %i flooded cells\n", nflooded); @@ -196,7 +231,7 @@ void animation(int size) // print_cluster_sizes(cell, ncells, cluster_sizes); - draw_configuration(cell, cluster_sizes, ncells, nx, ny, size, ncells); + draw_configuration(cell, cluster_sizes, ncells, nx, ny, nz, size, ncells); print_p(p); if (PRINT_LARGEST_CLUSTER_SIZE) print_largest_cluster_size(maxclustersize); @@ -265,8 +300,8 @@ void display(void) // animation(64); // animation(32); // animation(16); -// animation(8); - animation(4); + animation(8); +// animation(4); // animation(2); // animation(1); diff --git a/rde.c b/rde.c index e1e8914..f32c86e 100644 --- a/rde.c +++ b/rde.c @@ -39,58 +39,48 @@ #include #include -#define MOVIE 1 /* set to 1 to generate movie */ +#define MOVIE 0 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ /* General geometrical parameters */ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1000 /* window height */ -// // #define NX 640 /* number of grid points on x axis */ -// // #define NY 360 /* number of grid points on y axis */ -// #define NX 600 /* number of grid points on x axis */ -// #define NY 300 /* number of grid points on y axis */ -// // #define NX 480 /* number of grid points on x axis */ -// // #define NY 240 /* number of grid points on y axis */ -// // #define NX 1920 /* number of grid points on x axis */ -// // #define NY 1000 /* number of grid points on y axis */ -// -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -// #define YMIN -1.041666667 -// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ +#define NX 960 /* number of grid points on x axis */ +#define NY 500 /* number of grid points on y axis */ +// #define NX 480 /* number of grid points on x axis */ +// #define NY 250 /* number of grid points on y axis */ + +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ +#define YMIN -1.041666667 +#define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ // #define WINWIDTH 1280 /* window width */ // #define WINHEIGHT 720 /* window height */ - -// #define NX 200 /* number of grid points on x axis */ -// #define NY 200 /* number of grid points on y axis */ -#define NX 500 /* number of grid points on x axis */ -#define NY 500 /* number of grid points on y axis */ - +// +// // #define NX 320 /* number of grid points on x axis */ +// // #define NY 180 /* number of grid points on y axis */ // #define NX 640 /* number of grid points on x axis */ // #define NY 360 /* number of grid points on y axis */ - -// #define NX 1280 /* number of grid points on x axis */ -// #define NY 720 /* number of grid points on y axis */ - -#define XMIN -1.8 -#define XMAX 1.8 /* x interval */ -#define YMIN -1.8 -#define YMAX 1.8 /* y interval for 9/16 aspect ratio */ +// +// // #define NX 1280 /* number of grid points on x axis */ +// // #define NY 720 /* number of grid points on y axis */ +// +// #define XMIN -2.0 +// #define XMAX 2.0 /* x interval */ +// #define YMIN -1.125 +// #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ /* Choice of simulated equation */ -#define RDE_EQUATION 5 /* choice of reaction term, see list in global_3d.c */ +#define RDE_EQUATION 6 /* choice of reaction term, see list in global_3d.c */ #define NFIELDS 2 /* number of fields in reaction-diffusion equation */ -#define NLAPLACIANS 2 /* number of fields for which to compute Laplacian */ -// #define RDE_EQUATION 4 /* choice of reaction term, see list in global_3d.c */ -// #define NFIELDS 3 /* number of fields in reaction-diffusion equation */ -// #define NLAPLACIANS 3 /* number of fields for which to compute Laplacian */ +#define NLAPLACIANS 1 /* number of fields for which to compute Laplacian */ -#define ADD_POTENTIAL 1 /* set to 1 to add a potential (for Schrodinger equation) */ -#define POTENTIAL 1 /* type of potential, see list in global_3d.c */ -#define ADD_MAGNETIC_FIELD 1 /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */ +#define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */ +#define ADD_MAGNETIC_FIELD 0 /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */ +#define POTENTIAL 7 /* type of potential or vector potential, see list in global_3d.c */ #define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */ @@ -134,10 +124,11 @@ /* Physical patameters of wave equation */ -#define DT 0.00000002 +// #define DT 0.00000002 // #define DT 0.00000003 // #define DT 0.000000011 -// #define DT 0.00000001 +#define DT 0.0000012 +// #define DT 0.000001 #define VISCOSITY 2.0 @@ -148,10 +139,18 @@ #define DELTA 0.1 /* time scale separation */ #define FHNA 1.0 /* parameter in FHN equation */ #define FHNC -0.01 /* parameter in FHN equation */ -#define K_HARMONIC 0.5 /* spring constant of harmonic potential */ +#define K_HARMONIC 1.0 /* spring constant of harmonic potential */ #define K_COULOMB 0.5 /* constant in Coulomb potential */ +#define V_MAZE 0.4 /* potential in walls of maze */ #define BZQ 0.0008 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */ +#define B_FIELD 10.0 /* magnetic field */ +#define AB_RADIUS 0.2 /* radius of region with magnetic field for Aharonov-Bohm effect */ +#define K_EULER 50.0 /* constant in stream function integration of Euler equation */ + +#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 T_OUT 2.0 /* outside temperature */ #define T_IN 0.0 /* inside temperature */ @@ -182,9 +181,10 @@ /* Parameters for length and speed of simulation */ -// #define NSTEPS 500 /* number of frames of movie */ -#define NSTEPS 1100 /* number of frames of movie */ -#define NVID 500 /* number of iterations between images displayed on screen */ +#define 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 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 */ @@ -203,22 +203,21 @@ /* Visualisation */ -#define PLOT_3D 1 /* controls whether plot is 2D or 3D */ +#define PLOT_3D 0 /* controls whether plot is 2D or 3D */ #define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ /* Plot type - color scheme */ -#define CPLOT 32 -// #define CPLOT 32 -#define CPLOT_B 31 +#define CPLOT 52 +#define CPLOT_B 51 /* Plot type - height of 3D plot */ -#define ZPLOT 32 /* z coordinate in 3D plot */ +#define ZPLOT 52 /* z coordinate in 3D plot */ // #define ZPLOT 32 /* z coordinate in 3D plot */ -#define ZPLOT_B 30 /* z coordinate in second 3D plot */ +#define ZPLOT_B 51 /* z coordinate in second 3D plot */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */ #define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ @@ -226,8 +225,8 @@ #define WRAP_ANGLE 1 /* experimental: wrap angle to [0, 2Pi) for interpolation in angle schemes */ #define FADE_IN_OBSTACLE 0 /* set to 1 to fade color inside obstacles */ #define DRAW_OUTSIDE_GRAY 0 /* experimental - draw outside of billiard in gray */ -#define ADD_POTENTIAL_TO_Z 0 /* set to 1 to add the external potential to z-coordinate of plot */ -#define ADD_POT_CONSTANT 1.0 /* constant in front of added potential */ +#define ADD_POTENTIAL_TO_Z 1 /* set to 1 to add the external potential to z-coordinate of plot */ +#define ADD_POT_CONSTANT 0.35 /* constant in front of added potential */ #define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */ @@ -255,7 +254,7 @@ /* Color schemes, see list in global_pdes.c */ #define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 0 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ @@ -265,7 +264,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 30.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_AMPLITUDE 0.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) */ @@ -277,13 +276,19 @@ #define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */ #define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ #define HUEMEAN 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 1.0 /* scaling factor for energy log representation */ -#define LOG_SHIFT 0.0 +#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 LOG_SHIFT 1.0 + +#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 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.5 /* scale of color scheme bar */ +#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 ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ @@ -312,38 +317,54 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {8.0, 8.0, 12.0}; /* location of observer for REP_PROJ_3D representation */ +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 1.25 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#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 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 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 */ /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ -#define VMAX 2.0 /* max value of wave amplitude */ +#define VMAX 10.0 /* max value of wave amplitude */ +#define TEST_GRADIENT 0 /* print norm squared of gradient */ #define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT) /* to save computing time, to be improved */ #define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX))) #define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_PROBABILITIES)||(PRINT_NOISE)) #include "global_pdes.c" +#include "global_3d.c" /* constants and global variables */ + +#include "sub_maze.c" #include "sub_wave.c" #include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ -#include "global_3d.c" /* constants and global variables */ #include "sub_wave_3d_rde.c" /* should be later replaced by sub_wave_rde.c */ #include "sub_rde.c" -double potential(int i, int j) -/* compute potential (e.g. for Schrödinger equation) */ +double f_aharonov_bohm(double r2) +/* radial part of Aharonov-Bohm vector potential */ { - double x, y, xy[2], r, small = 2.0e-1, kx, ky, lx = XMAX - XMIN, r1, r2, r3; + double r02 = AB_RADIUS*AB_RADIUS; + + if (r2 > r02) return(-0.25*r02/r2); + else return(0.25*(r2 - 2.0*r02)/r02); + +// if (r2 > r02) return(1.0/r2); +// else return((2.0*r02 - r2)/(r02*r02)); +} + +double potential(int i, int j) +/* compute potential (e.g. for Schrödinger equation), or potential part if there is a magnetic field */ +{ + double x, y, xy[2], r, small = 1.0e-1, kx, ky, lx = XMAX - XMIN, r1, r2, r3, f; + int rect; ij_to_xy(i, j, xy); x = xy[0]; @@ -380,6 +401,23 @@ double potential(int i, int j) // r = r/3.0; return (-0.5*K_COULOMB*(1.0/r1 + 1.0/r2 + 1.0/r3)); } + case (VPOT_CONSTANT_FIELD): + { + return (K_HARMONIC*(x*x + y*y)); /* magnetic field strength b is chosen such that b^2 = K_HARMONIC */ + } + case (VPOT_AHARONOV_BOHM): + { + r2 = x*x + y*y; + f = f_aharonov_bohm(r2); + return (B_FIELD*B_FIELD*f*f*r2); /* magnetic field strength b is chosen such that b^2 = K_HARMONIC */ +// return (K_HARMONIC*f); /* magnetic field strength b is chosen such that b^2 = K_HARMONIC */ + } + case (POT_MAZE): + { + for (rect=0; rect= SMOOTHEN_PERIOD) smooth = 0; + } + } + + if (TEST_GRADIENT) { + for (i=0; i<2*NX*NY; i++){ + test += nabla_omega[i]*nabla_omega[i]; + test += nabla_psi[i]*nabla_psi[i]; + } + printf("nabla square = %.5lg\n", test/((double)NX*NY)); + } + + #pragma omp parallel for private(i,j,k,x,y,z,deltax,deltay,deltaz,rho) for (i=0; i0)||(j!=NYMAZE/2))&&(maze[n].west)) add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0); + if (maze[n].south) add_rectangle_to_segments(x1, y1, x1 + dx, y1 - width, segment, 0); + } + + /* top side of maze */ + add_rectangle_to_segments(YMIN + padding + MAZE_XSHIFT, YMAX - padding, YMAX - padding + MAZE_XSHIFT, YMAX - padding - width, segment, 0); + + /* right side of maze */ + y1 = YMIN + padding + dy*((double)NYMAZE/2); + x1 = YMAX - padding + MAZE_XSHIFT; + add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, y1 - dy, segment, 0); + add_rectangle_to_segments(x1, y1, x1 - width, YMAX + 1.0, segment, 0); + + /* left side of maze */ + x1 = YMIN + padding + MAZE_XSHIFT; + add_rectangle_to_segments(x1, YMIN - 1.0, x1 - width, YMIN + padding, segment, 0); + add_rectangle_to_segments(x1, YMAX - padding, x1 - width, YMAX + 1.0, segment, 0); + + free(maze); +} + + void init_segment_config(t_segment segment[NMAXSEGMENTS]) /* initialise linear obstacle configuration */ { - int i, j, cycle = 0, iminus, iplus, nsides, n; - double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy; + int i, j, cycle = 0, iminus, iplus, nsides, n, concave = 1; + double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, dy; switch (SEGMENT_PATTERN) { case (S_RECTANGLE): @@ -1317,13 +1415,13 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS]) width = MU; lpocket = 0.1; - add_rectangle_to_segments(BCXMIN + lpocket, BCYMIN, xmid - lpocket, BCYMIN - width, segment); - add_rectangle_to_segments(xmid + lpocket, BCYMIN, BCXMAX - lpocket, BCYMIN - width, segment); - add_rectangle_to_segments(BCXMAX + width, BCYMIN + lpocket, BCXMAX, BCYMAX - lpocket, segment); + add_rectangle_to_segments(BCXMIN + lpocket, BCYMIN, xmid - lpocket, BCYMIN - width, segment, 0); + add_rectangle_to_segments(xmid + lpocket, BCYMIN, BCXMAX - lpocket, BCYMIN - width, segment, 0); + add_rectangle_to_segments(BCXMAX + width, BCYMIN + lpocket, BCXMAX, BCYMAX - lpocket, segment, 0); - add_rectangle_to_segments(BCXMAX - lpocket, BCYMAX, xmid + lpocket, BCYMAX + width, segment); - add_rectangle_to_segments(xmid - lpocket, BCYMAX, BCXMIN + lpocket, BCYMAX + width, segment); - add_rectangle_to_segments(BCXMIN - width, BCYMAX - lpocket, BCXMIN, BCYMIN + lpocket, segment); + add_rectangle_to_segments(BCXMAX - lpocket, BCYMAX, xmid + lpocket, BCYMAX + width, segment, 0); + add_rectangle_to_segments(xmid - lpocket, BCYMAX, BCXMIN + lpocket, BCYMAX + width, segment, 0); + add_rectangle_to_segments(BCXMIN - width, BCYMAX - lpocket, BCXMIN, BCYMIN + lpocket, segment, 0); cycle = 0; for (i=0; i 0.05)&&(y1 < b - 0.05)); + } + case (RCK_RECT_HAT) : + { + a = 0.5*LAMBDA; + b = 0.49*PI*LAMBDA; + y1 = y - YMIN - LAMBDA; + if (vabs(x) > 0.95*a) return(0); + if (y1 < 0.05) return(0); + if (y1 < b - 0.05) return(1); + return(y1 < a + b - 0.05 - vabs(x)); + } + } +} int in_segment_region(double x, double y) /* returns 1 if (x,y) is inside region delimited by obstacle segments */ { - double angle, dx, height, width, theta, lx, ly, x1, y1, x2, y2; + double angle, dx, height, width, theta, lx, ly, x1, y1, x2, y2, padding; if (x >= BCXMAX) return(0); if (x <= BCXMIN) return(0); @@ -1804,14 +2008,14 @@ int in_segment_region(double x, double y) y1 = y - ysegments[0]; if (y1 < YMIN + LAMBDA) return(0); else if (y1 > YMIN + 2.4*LAMBDA) return(0); - else - { - ly = 0.7*LAMBDA; - lx = 0.7*LAMBDA; - x1 = x - xsegments[0]; - y1 -= YMIN + 1.7*LAMBDA; - if (x1*x1/(lx*lx) + y1*y1/(ly*ly) + MU*MU < 0.925) return(1); - } + else if (in_rocket(x - xsegments[0], y1, ROCKET_SHAPE)) return(1); +// { +// ly = 0.7*LAMBDA; +// lx = 0.7*LAMBDA; +// x1 = x - xsegments[0]; +// y1 -= YMIN + 1.7*LAMBDA; +// if (x1*x1/(lx*lx) + y1*y1/(ly*ly) + MU*MU < 0.925) return(1); +// } return(0); } case (S_TWO_ROCKETS): @@ -1819,17 +2023,11 @@ int in_segment_region(double x, double y) y1 = y - ysegments[0]; y2 = y - ysegments[1]; if ((y1 < YMIN + LAMBDA)&&(y2 < YMIN + LAMBDA)) return(0); - else if ((y1 > YMIN + 2.4*LAMBDA)&&(y2 > YMIN + 2.4*LAMBDA)) return(0); + else if ((y1 > YMIN + 3.5*LAMBDA)&&(y2 > YMIN + 3.5*LAMBDA)) return(0); else { - ly = 0.7*LAMBDA; - lx = 0.7*LAMBDA; - x1 = x - xsegments[0]; - y1 -= YMIN + 1.7*LAMBDA; - if (x1*x1/(lx*lx) + y1*y1/(ly*ly) + MU*MU < 0.925) return(1); - x2 = x - xsegments[1]; - y2 -= YMIN + 1.7*LAMBDA; - if (x2*x2/(lx*lx) + y2*y2/(ly*ly) + MU*MU < 0.925) return(1); + if (in_rocket(x - xsegments[0], y1, ROCKET_SHAPE)) return(1); + if (in_rocket(x - xsegments[1], y2, ROCKET_SHAPE_B)) return(1); } return(0); } @@ -1839,6 +2037,13 @@ int in_segment_region(double x, double y) else if (y > LAMBDA) return(1); else return(0); } + case (S_EXT_RECTANGLE): + { + padding = 0.1; + if (vabs(x) > LAMBDA + padding) return(1); + else if (vabs(y) > 0.1*LAMBDA + padding) return(1); + else return(0); + } default: return(1); } } @@ -1901,6 +2106,69 @@ void translate_segments(t_segment segment[NMAXSEGMENTS], double deltax[2], doubl } } +void translate_one_segment(t_segment segment[NMAXSEGMENTS], int i, double deltax, double deltay) +/* translates the repelling segment by given vector */ +{ + segment[i].x1 += deltax; + segment[i].x2 += deltax; + + segment[i].y1 += deltay; + segment[i].y2 += deltay; + segment[i].c = segment[i].nx*segment[i].x1 + segment[i].ny*segment[i].y1; + + segment[i].xc += deltax; + segment[i].yc += deltay; + +} + +void rotate_one_segment(t_segment segment[NMAXSEGMENTS], int i, double dalpha, double xc, double yc) +/* rotates the repelling segment by given angle around (xc, yc) */ +{ + double ca, sa, x, y, nx, ny; + + ca = cos(dalpha); + sa = sin(dalpha); + + x = segment[i].x1 - xc; + y = segment[i].y1 - yc; + + segment[i].x1 = xc + x*ca - y*sa; + segment[i].y1 = yc + x*sa + y*ca; + + x = segment[i].x2 - xc; + y = segment[i].y2 - yc; + + segment[i].x2 = xc + x*ca - y*sa; + segment[i].y2 = yc + x*sa + y*ca; + + segment[i].xc = 0.5*(segment[i].x1 + segment[i].x2); + segment[i].yc = 0.5*(segment[i].y1 + segment[i].y2); + + nx = segment[i].nx; + ny = segment[i].ny; + + segment[i].nx = ca*nx - sa*ny; + segment[i].ny = sa*nx + ca*ny; + + segment[i].c = segment[i].nx*segment[i].x1 + segment[i].ny*segment[i].y1; + + if (segment[i].concave) + { + segment[i].angle1 += dalpha; + segment[i].angle2 += dalpha; + while (segment[i].angle1 > DPI) + { + segment[i].angle1 -= DPI; + segment[i].angle2 -= DPI; + } + while (segment[i].angle2 < 0.0) + { + segment[i].angle1 += DPI; + segment[i].angle2 += DPI; + } + } +} + /* Computation of interaction force */ double lennard_jones_force(double r, t_particle particle) @@ -2899,6 +3167,31 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl } +void set_segment_group_color(int group, double lum, double rgb[3]) +{ + switch (group) { + case (1): + { + hsl_to_rgb_palette(270.0, 0.9, 0.5, rgb, COLOR_PALETTE); + break; + } + case (2): + { + hsl_to_rgb_palette(90.0, 0.9, 0.5, rgb, COLOR_PALETTE); + break; + } + default: + { + rgb[0] = 1.0; + rgb[1] = 1.0; + rgb[2] = 1.0; + } + } + + glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); +} + + void draw_one_triangle(t_particle particle, int same_table[9*HASHMAX], int p, int q, int nsame) { double x, y, dx, dy; @@ -2906,6 +3199,12 @@ void draw_one_triangle(t_particle particle, int same_table[9*HASHMAX], int p, in x = particle.xc + (double)p*(BCXMAX - BCXMIN); y = particle.yc + (double)q*(BCYMAX - BCYMIN); + + if (TRACK_SEGMENT_GROUPS) + { + x -= xtrack; + y -= ytrack; + } glBegin(GL_TRIANGLE_FAN); glVertex2d(x, y); @@ -3004,15 +3303,20 @@ void draw_one_particle_links(t_particle particle) void draw_one_particle(t_particle particle, double xc, double yc, double radius, double angle, int nsides, double width, double rgb[3]) /* draw one of the particles */ { - double ca, sa, x1, x2, y1, y2, xc1, wangle; + double ca, sa, x1, x2, y1, y2, xc1, yc1, wangle; int wsign; if (CENTER_VIEW_ON_OBSTACLE) xc1 = xc - xshift; else xc1 = xc; + if (TRACK_SEGMENT_GROUPS) + { + xc1 -= xtrack; + yc1 = yc - ytrack; + } glColor3f(rgb[0], rgb[1], rgb[2]); if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)) - draw_colored_rhombus(xc1, yc, radius, angle + APOLY*PID, rgb); - else draw_colored_polygon(xc1, yc, radius, nsides, angle + APOLY*PID, rgb); + draw_colored_rhombus(xc1, yc1, radius, angle + APOLY*PID, rgb); + else draw_colored_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID, rgb); /* draw crosses on particles of second type */ if ((TWO_TYPES)&&(DRAW_CROSS)) @@ -3025,28 +3329,28 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, glLineWidth(3); glColor3f(0.0, 0.0, 0.0); x1 = xc1 - MU_B*ca; - y1 = yc - MU_B*sa; + y1 = yc1 - MU_B*sa; x2 = xc1 + MU_B*ca; - y2 = yc + MU_B*sa; + y2 = yc1 + MU_B*sa; draw_line(x1, y1, x2, y2); x1 = xc1 - MU_B*sa; - y1 = yc + MU_B*ca; + y1 = yc1 + MU_B*ca; x2 = xc1 + MU_B*sa; - y2 = yc - MU_B*ca; + y2 = yc1 - MU_B*ca; draw_line(x1, y1, x2, y2); } glLineWidth(width); glColor3f(1.0, 1.0, 1.0); if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)) - draw_rhombus(xc1, yc, radius, angle + APOLY*PID); - else draw_polygon(xc1, yc, radius, nsides, angle + APOLY*PID); + draw_rhombus(xc1, yc1, radius, angle + APOLY*PID); + else draw_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID); if (particle.interaction == I_LJ_WATER) for (wsign = -1; wsign <= 1; wsign+=2) { wangle = particle.angle + (double)wsign*DPI/3.0; x1 = xc1 + particle.radius*cos(wangle); - y1 = yc + particle.radius*sin(wangle); + y1 = yc1 + particle.radius*sin(wangle); draw_colored_polygon(x1, y1, 0.5*radius, nsides, angle + APOLY*PID, rgb); glColor3f(1.0, 1.0, 1.0); draw_polygon(x1, y1, 0.5*radius, nsides, angle + APOLY*PID); @@ -3310,17 +3614,20 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] glLineWidth(CONTAINER_WIDTH); hsl_to_rgb(300.0, 0.1, 0.5, rgb); for (i = 0; i < nobstacles; i++) - draw_colored_circle(obstacle[i].xc, obstacle[i].yc, obstacle[i].radius, NSEG, rgb); + draw_colored_circle(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius, NSEG, rgb); glColor3f(1.0, 1.0, 1.0); for (i = 0; i < nobstacles; i++) - draw_circle(obstacle[i].xc, obstacle[i].yc, obstacle[i].radius, NSEG); + draw_circle(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius, NSEG); } if (ADD_FIXED_SEGMENTS) { glLineWidth(CONTAINER_WIDTH); glColor3f(1.0, 1.0, 1.0); for (i = 0; i < nsegments; i++) if (segment[i].active) - draw_line(segment[i].x1, segment[i].y1, segment[i].x2, segment[i].y2); + { + if (COLOR_SEG_GROUPS) set_segment_group_color(segment[i].group, 1.0, rgb); + draw_line(segment[i].x1 - xtrack, segment[i].y1 - ytrack, segment[i].x2 - xtrack, segment[i].y2 - ytrack); + } } switch (BOUNDARY_COND) { @@ -3547,7 +3854,7 @@ void print_parameters(double beta, double temperature, double krepel, double len xxtext = XMIN + 0.08*scale; } xmid = 0.5*(XMIN + XMAX) - 0.1*scale; - xmidtext = xmid - 0.24*scale; + xmidtext = xmid - 0.3*scale; for (i=0; i=t_window) position[group] = 0; + for (i=0; i segment[i].angle2) angle -= DPI; + r = 1.5*particle[j].radius; if ((distance < r)&&(angle > segment[i].angle1)&&(angle < segment[i].angle2)) @@ -4041,12 +4424,11 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl f = KSPRING_OBSTACLE*(r - distance); particle[j].fx += f*cos(angle); particle[j].fy += f*sin(angle); - if (MOVE_BOUNDARY) + if ((MOVE_BOUNDARY)||(MOVE_SEGMENT_GROUPS)) { segment[i].fx -= f*cos(angle); segment[i].fy -= f*sin(angle); -// segment[i].fx += f*cos(angle); -// segment[i].fy += f*sin(angle); + segment[i].torque -= (x - segment[i].xc)*f*sin(angle) - (y - segment[i].yc)*f*cos(angle); } } } @@ -4941,4 +5323,228 @@ void update_types(t_particle particle[NMAXCIRCLES]) } } } - \ No newline at end of file + + +double plot_coord(double x, double xmin, double xmax) +{ + return(xmin + x*(xmax - xmin)); +} + + +void draw_speed_plot(t_group_data *group_speeds, int i) +/* draw plot of obstacle speeds as a function of time */ +{ + int j, group; + char message[100]; + static double xmin, xmax, ymin, ymax, xmid, ymid, dx, dy, plotxmin, plotxmax, plotymin, plotymax; + double pos[2], x1, y1, x2, y2, rgb[3]; + static int first = 1, gshift = INITIAL_TIME + NSTEPS; + + if (first) + { +// xmin = XMAX - 1.35; +// xmax = XMAX - 0.05; +// ymin = YMAX - 1.8; +// ymax = YMAX - 0.5; + + xmin = XMAX - 1.8; + xmax = XMAX - 0.066; + ymin = YMAX - 2.5; + ymax = YMAX - 0.76; + + xmid = 0.5*(xmin + xmax); + ymid = 0.5*(ymin + ymax); + + dx = 0.5*(xmax - xmin); + dy = 0.5*(ymax - ymin); + + plotxmin = xmin + 0.05; + plotxmax = xmax - 0.1; + plotymin = ymin + 0.07; + plotymax = ymax - 0.15; + + first = 0; + } + +// rgb[0] = 1.0; rgb[1] = 1.0; rgb[2] = 1.0; + + glLineWidth(2); + + /* plot angular speed */ + for (group=1; group scaley) scaley = y0; + + x2 = plot_coord(0.5 + x0/scalex, plotxmin, plotxmax); + y2 = plot_coord(y0/scaley + 0.05, plotymin, plotymax); + +// printf("yinitial = %.3lg, (x0, y0) = (%.3lg, %.3lg), (x2, y2) = (%.3lg, %.3lg)\n", yinitial, x0, y0, x2, y2); + + if ((j>0)&&(module2(x2-x1, y2-y1) < 0.1)) draw_line(x1, y1, x2, y2); + x1 = x2; + y1 = y2; + } + if (i>0) draw_colored_circle(x1, y1, 0.015, NSEG, rgb); + } + + glColor3f(1.0, 1.0, 1.0); + + /* axes and labels */ + draw_line(plotxmin, plotymin, plotxmax + 0.05, plotymin); + draw_line(plotxmin, plotymin, plotxmin, plotymax + 0.1); + + for (j=1; j<(int)scaley; j++) + { + y1 = plot_coord((double)j/scaley, plotymin, plotymax); + draw_line(plotxmin - 0.02, y1, plotxmin + 0.02, y1); + } + + sprintf(message, "%i", (int)scaley - 1); + write_text_fixedwidth(plotxmin - 0.15, y1 - 0.025, message); + +// sprintf(message, "%.1f", VMAX_PLOT_SPEEDS); + sprintf(message, "y"); + write_text_fixedwidth(plotxmin - 0.1, plotymax - 0.01, message); + + sprintf(message, "x"); + write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.1, message); +} + + +void init_segment_group(t_segment segment[NMAXSEGMENTS], int group, t_group_segments segment_group[NMAXGROUPS]) +/* initialize center of mass and similar data of grouped segments */ +{ + int i, nseg_group = 0; + double xc = 0.0, yc = 0.0; + + for (i=0; i 0)&&(!maze[n].tested)) + while ((npaths > 0)) + { + maze[n].active = 1; + + printf("Cell (%i, %i) ", n%NXMAZE, n/NXMAZE); + + nnext = 0; + for (i=0; i= 12) + { + free(image); + mem_counter = 0; + } + return 0; } @@ -198,6 +208,20 @@ void save_frame() } +void save_frame_counter(int counter) +{ + char *name="part.", n2[100]; + char format[6]=".%05i"; + + strcpy(n2, name); + sprintf(strstr(n2,"."), format, counter); + strcat(n2, ".tif"); + printf(" saving frame %s \n",n2); + writetiff(n2, "Billiard in an ellipse", 0, 0, + WINWIDTH, WINHEIGHT, COMPRESSION_LZW); + +} + void write_text_fixedwidth( double x, double y, char *st) { @@ -5353,7 +5377,8 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ case D_POLYLINE: { /* not easy to implement for non-convex polygons */ - return(1); + if (POLYLINE_PATTERN == P_MAZE) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); + else return(1); break; } default: @@ -5722,7 +5747,9 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) { int i, j, k, l, n, z, ii, jj, terni[SDEPTH], ternj[SDEPTH], quater[SDEPTH], cond; short int vkoch[NMAXCIRCLES], turnright; - double ratio, omega, angle, sw, length, dist, x, y, ta, tb, a, b; + double ratio, omega, angle, sw, length, dist, x, y, ta, tb, a, b, + x1, y1, x2, y2, dx, dy, padding = 0.02, width = 0.01; + t_maze* maze; switch (POLYLINE_PATTERN) { case (P_RECTANGLE): @@ -6148,6 +6175,96 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) } break; } + case (P_MAZE): + { + maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze)); + + init_maze(maze); + + /* build walls of maze */ + dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE); + dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE); + + nsides = 0; + ncircles = 0; + + for (i=0; i0)||(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; + +// add_rectangle_to_segments(x1, y1, x1 - width, y1 + dy, segment, 0); + } + 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; + +// add_rectangle_to_segments(x1, y1, x1 + dx, y1 - width, segment, 0); + } + + nsides++; + } + + /* top side of maze */ + polyline[nsides].x1 = YMIN + padding + MAZE_XSHIFT; + polyline[nsides].y1 = YMAX - padding; + polyline[nsides].x2 = YMAX - padding + MAZE_XSHIFT; + polyline[nsides].y2 = YMAX - padding; + polyline[nsides].angle = 0.0; + nsides++; + + /* right side of maze */ + y1 = YMIN + padding + dy*((double)NYMAZE/2); + x1 = YMAX - padding + MAZE_XSHIFT; + polyline[nsides].x1 = x1; + polyline[nsides].y1 = YMIN - 1.0; + polyline[nsides].x2 = x1; + polyline[nsides].y2 = y1 - dy; + polyline[nsides].angle = PID; + nsides++; + + polyline[nsides].x1 = x1; + polyline[nsides].y1 = y1; + polyline[nsides].x2 = x1; + polyline[nsides].y2 = YMAX + 1.0; + polyline[nsides].angle = PID; + nsides++; + + /* left side of maze */ + x1 = YMIN + padding + MAZE_XSHIFT; + polyline[nsides].x1 = x1; + polyline[nsides].y1 = YMIN - 1.0; + polyline[nsides].x2 = x1; + polyline[nsides].y2 = YMIN + padding; + polyline[nsides].angle = PID; + nsides++; + + polyline[nsides].x1 = x1; + polyline[nsides].y1 = YMAX - padding; + polyline[nsides].x2 = x1; + polyline[nsides].y2 = YMAX + 1.0; + polyline[nsides].angle = PID; + nsides++; + + free(maze); + break; + } } } @@ -6249,3 +6366,63 @@ void init_polyline_depth(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRC } +int test_initial_condition(double *configs[NPARTMAX], int active[NPARTMAX], int color[NPARTMAX]) +/* apply a test to initial conditions - so far, whether particle has left maze on the right */ +{ + int i, j, time, nactive = 0, counter = 0, tmax[NPARTMAX], tmaxmax = 0, tmaxmin = 1000; + double conf[8], newconf[8], cosphi, x2, pcolor; + + for (i=0; i 0.0) + { + active[i] = 1; + if (time > tmaxmax) tmaxmax = time; + else if (time < tmaxmin) tmaxmin = time; + } + else active[i] = 0; + + nactive += active[i]; + } + + printf("%i particles, %i active particles\n", nparticles, nactive); + printf("tmin = %i, tmax = %i\n", tmaxmin, tmaxmax); +// sleep(5); + + /* reorder particles */ + for (i=0; i 0.5) return(pstar + factor*(1.0 - pstar)*ipow(time - 0.5, power)); - else return(pstar - factor*pstar*ipow(0.5 - time, power)); + if (time > 0.5) return(pstar + factor*(1.0 - pstar)*ipow(time - 0.5, P_SCHEDULE_POWER)); + else return(pstar - factor*pstar*ipow(0.5 - time, P_SCHEDULE_POWER)); } int in_plot_box(double x, double y) @@ -496,12 +531,17 @@ int in_plot_box_screencoord(double x, double y) double size_ratio_color(int clustersize, int ncells) /* color of cell as function of the size of its cluster */ { - double ratio, minratio = 1.0e-2, x; + double ratio, minratio = 1.0e-2, x, p = 0.1; + + minratio = 1.0/(double)ncells; ratio = (double)clustersize/(double)ncells; if (ratio > 1.0) ratio = 1.0; else if (ratio < minratio) ratio = minratio; - x = log(ratio/minratio)/log(1.0/minratio); +// x = log(ratio/minratio)/log(1.0/minratio); +// x = log(1.0 + log(ratio/minratio))/log(1.0 - log(minratio)); +// x = pow(log(ratio/minratio)/(-log(minratio)), 0.1); + x = (pow(ratio, p) - pow(minratio, p))/(1.0 - pow(minratio, p)); return(CLUSTER_HUEMIN*x + CLUSTER_HUEMAX*(1.0 -x)); /* other attempts that seem to bug */ @@ -513,17 +553,30 @@ double size_ratio_color(int clustersize, int ncells) // return(CLUSTER_HUEMAX*ratio + CLUSTER_HUEMIN*(1.0 -ratio)); } -void set_cell_color(t_perco cell, int *cluster_sizes, int fade, int max_cluster_size) +void compute_cell_color(t_perco cell, int *cluster_sizes, int fade, int max_cluster_size, int kx, int nx, int kz, int nz, double rgb[3]) /* compute color of cell */ { int k, color, csize; - double rgb[3], fade_factor = 0.15, hue; + double fade_factor = 0.15, hue; if (!cell.open) hsl_to_rgb_palette(HUE_CLOSED, 0.9, 0.5, rgb, COLOR_PALETTE); else { if (!cell.flooded) hsl_to_rgb_palette(HUE_OPEN, 0.9, 0.5, rgb, COLOR_PALETTE); - else hsl_to_rgb_palette(HUE_FLOODED, 0.9, 0.5, rgb, COLOR_PALETTE); + else + { + if (COLOR_CELLS_BY_XCOORD) + { + hue = CLUSTER_HUEMIN + (CLUSTER_HUEMAX - CLUSTER_HUEMIN)*(double)kx/(double)nx; + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE); + } + else if (COLOR_CELLS_BY_ZCOORD) + { + hue = CLUSTER_HUEMIN + (CLUSTER_HUEMAX - CLUSTER_HUEMIN)*(double)kz/(double)nz; + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE); + } + else hsl_to_rgb_palette(HUE_FLOODED, 0.9, 0.5, rgb, COLOR_PALETTE); + } if ((FIND_ALL_CLUSTERS)&&(COLOR_CLUSTERS_BY_SIZE)) { @@ -537,13 +590,22 @@ void set_cell_color(t_perco cell, int *cluster_sizes, int fade, int max_cluster_ hue = CLUSTER_HUEMIN + (CLUSTER_HUEMAX - CLUSTER_HUEMIN)*(double)color/(double)N_CLUSTER_COLORS; hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE); } + +// if ((FLOOD_LEFT_BOUNDARY)&&(cell.flooded == 1)) hsl_to_rgb_palette(HUE_FLOODED, 0.9, 0.5, rgb, COLOR_PALETTE); } if (fade) for (k=0; k<3; k++) rgb[k] = 1.0 - fade_factor + fade_factor*rgb[k]; - - glColor3f(rgb[0], rgb[1], rgb[2]); + } +void set_cell_color(t_perco cell, int *cluster_sizes, int fade, int max_cluster_size, int i, int nx, int k, int nz) +/* set color of cell */ +{ + double rgb[3]; + + compute_cell_color(cell, cluster_sizes, fade, max_cluster_size, i, nx, k, nz, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); +} double plot_coord(double x, double xmin, double xmax) { @@ -598,15 +660,18 @@ void draw_size_plot(double plot_cluster_size[NSTEPS], int i, double pcrit) sprintf(message, "p"); write_text_fixedwidth(pos[0], pos[1], message); - xy_to_pos(plotxmax - 0.015, plotymin - 0.06, pos); - sprintf(message, "1"); - write_text_fixedwidth(pos[0], pos[1], message); - xy_to_pos(x - 0.02, plotymin - 0.04, pos); sprintf(message, "pc"); write_text_fixedwidth(pos[0], pos[1], message); - xy_to_pos(plotxmin + 0.02, plotymax + 0.05, pos); + hsl_to_rgb_palette(HUE_GRAPH_SIZE, 0.9, 0.5, rgb, COLOR_PALETTE); + glColor3f(rgb[0], rgb[1], rgb[2]); + + xy_to_pos(plotxmax - 0.015, plotymin - 0.06, pos); + sprintf(message, "1"); + write_text_fixedwidth(pos[0], pos[1], message); + + xy_to_pos(plotxmin + 0.02, plotymax + 0.1, pos); sprintf(message, "nflooded/nopen"); write_text_fixedwidth(pos[0], pos[1], message); @@ -615,8 +680,6 @@ void draw_size_plot(double plot_cluster_size[NSTEPS], int i, double pcrit) write_text_fixedwidth(pos[0], pos[1], message); /* plot */ - hsl_to_rgb_palette(HUE_GRAPH, 0.9, 0.5, rgb, COLOR_PALETTE); - glColor3f(rgb[0], rgb[1], rgb[2]); x1 = plotxmin; y1 = plotymin; for (j=0; j 0) { - bin = (cluster_sizes[i]-1)/binwidth; + if (HISTO_X_LOG_SCALE) + { + bin = (int)(log(logfactor*(double)cluster_sizes[i])/logbinwidth) - 1; + if (bin >= nbins) bin = nbins - 1; + else if (bin < 0) bin = 0; + } + else bin = (cluster_sizes[i]-1)/binwidth; +// printf("cluster size = %i, bin = %i\n", cluster_sizes[i], bin); histo[bin]++; } for (bin=0; bin maxheight) maxheight = histo[bin]; @@ -772,20 +849,25 @@ void draw_cluster_histogram(int ncells, int *cluster_sizes, int maxclustersize, x1 = plotxmin; y1 = plotymin; - max_x_axis = maxclustersize; + if (HISTO_X_LOG_SCALE) max_x_axis = log((double)maxclustersize); + else max_x_axis = maxclustersize; if (max_x_axis < HISTO_BINS) max_x_axis = HISTO_BINS; for (bin=0; bin < nbins; bin++) { - csize = bin*binwidth + binwidth/2; +// csize = bin*binwidth + binwidth/2; + if (HISTO_X_LOG_SCALE) csize = (int)(exp((double)bin*logbinwidth)/logfactor); + else csize = bin*binwidth; + if (csize >= ncells) csize = ncells-1; hue = size_ratio_color(csize, ncells); - x2 = plot_coord((double)((bin+1)*binwidth)/(double)(max_x_axis+1), plotxmin, plotxmax); + x2 = plot_coord((double)((bin+1)*binwidth)/(double)(max_x_axis), plotxmin, plotxmax); // x2 = plot_coord((double)(bin+1)/(double)nbins, plotxmin, plotxmax); y = log((double)(histo[bin]+1))/log((double)(maxheight+1)); y2 = plot_coord(y, plotymin, plotymax); +// printf("x1 = %.2f, x2 %.2f\n", x1, x2); draw_colored_rectangle(x1, y1, x2, y2, hue); x1 = x2; } @@ -808,33 +890,54 @@ void draw_cluster_histogram(int ncells, int *cluster_sizes, int maxclustersize, } /* graduation of x axis */ - x = log((double)(max_x_axis+1))/log(10.0); - n = ipowi(10, (int)x); - y = (double)n/10.0; - - delta = plot_coord((double)n/((double)(max_x_axis+1)), plotxmin, plotxmax) - plotxmin; - if (delta > 0.13 + 0.01*x) di = 1; - else if (delta > 0.08 + 0.01*x) di = 2; - else di = 5; - for (i=di; i<10; i+=di) + if (HISTO_X_LOG_SCALE) { - x1 = plot_coord(y*(double)i*10.0/(double)(max_x_axis+1), plotxmin, plotxmax); - if (i*n < max_x_axis*11/10) + y = log(logfactor*(double)(maxclustersize+1))/log(10.0); + printf("y = %.3lg\n", y); + for (i=1; i<(int)y + 1; i++) { - sprintf(message, "%i", i*n); - xy_to_pos(x1 + 0.005 - 0.012*x, plotymin - 0.07, pos); + n = ipowi(10, i); + x = log((double)(n+1))/(log(logfactor*(double)(maxclustersize+1))); +// printf("n = %i, x = %.3lg\n", n, x); + x1 = plot_coord(x, plotxmin, plotxmax); + xy_to_pos(x1, plotymin - 0.1, pos); + draw_line(x1, plotymin - 0.02, x1, plotymin + 0.02); + if (n <= 1000) sprintf(message, "%i", n); + else sprintf(message, "1e%i", i); + xy_to_pos(x1 - 0.015, plotymin - 0.05, pos); write_text_fixedwidth(pos[0], pos[1], message); } } - - delta = plot_coord((double)n/(10.0*(double)(max_x_axis+1)), plotxmin, plotxmax) - plotxmin; - for (i=0; i<100; i++) if (i*n < 11*max_x_axis) + else { - y = (double)(i*n)/10.0; - x1 = plot_coord(y/(double)(max_x_axis+1), plotxmin, plotxmax); - xy_to_pos(x1, plotymin , pos); - if (i%10 == 0) draw_line(x1, plotymin - 0.02, x1, plotymin + 0.02); - else if (delta > 0.02) draw_line(x1, plotymin - 0.01, x1, plotymin + 0.01); + x = log((double)(max_x_axis+1))/log(10.0); + n = ipowi(10, (int)x); + y = (double)n/10.0; + + delta = plot_coord((double)n/((double)(max_x_axis+1)), plotxmin, plotxmax) - plotxmin; + if (delta > 0.13 + 0.01*x) di = 1; + else if (delta > 0.08 + 0.01*x) di = 2; + else di = 5; + for (i=di; i<10; i+=di) + { + x1 = plot_coord(y*(double)i*10.0/(double)(max_x_axis+1), plotxmin, plotxmax); + if (i*n < max_x_axis*11/10) + { + sprintf(message, "%i", i*n); + xy_to_pos(x1 + 0.005 - 0.012*x, plotymin - 0.07, pos); + write_text_fixedwidth(pos[0], pos[1], message); + } + } + + delta = plot_coord((double)n/(10.0*(double)(max_x_axis+1)), plotxmin, plotxmax) - plotxmin; + for (i=0; i<100; i++) if (i*n < 11*max_x_axis) + { + y = (double)(i*n)/10.0; + x1 = plot_coord(y/(double)(max_x_axis+1), plotxmin, plotxmax); + xy_to_pos(x1, plotymin , pos); + if (i%10 == 0) draw_line(x1, plotymin - 0.02, x1, plotymin + 0.02); + else if (delta > 0.02) draw_line(x1, plotymin - 0.01, x1, plotymin + 0.01); + } } /* for debugging */ @@ -1010,12 +1113,33 @@ void test_neighbours(int start, t_perco *cell, int nx, int ny, int size, int nce } +void draw_cube_ijk(int i, int j, int k, t_perco *cell, int *cluster_sizes, int nx, int ny, int nz, int size, int max_cluster_size) +/* draw one cube of 3d configuration */ +{ + double dx, x, y, z, rgb[3]; + + dx = 1.0/(double)nx; + + compute_cell_color(cell[k*nx*ny+j*nx+i], cluster_sizes, 0, max_cluster_size, i, nx, k, nz, rgb); + + x = (double)i*dx - 0.5; + y = (double)j*dx - 0.5; + z = (double)k*dx - 0.5; + draw_cube(x, y, z, dx, rgb); +} + +int plot_cube(t_perco cell) +/* returns 1 if cube is plotted */ +{ + if (PLOT_ONLY_FLOODED_CELLS) return(cell.flooded); + else return(cell.open); +} -void draw_configuration(t_perco *cell, int *cluster_sizes, int ncells, int nx, int ny, int size, int max_cluster_size) +void draw_configuration(t_perco *cell, int *cluster_sizes, int ncells, int nx, int ny, int nz, int size, int max_cluster_size) /* draw the configuration */ { - int i, j, ishift, k, n; - double rgb[3], x, y, alpha, r, fade = 0, dsize, radius, x2, y2; + int i, j, ishift, k, n, sector; + double rgb[3], x, y, z, alpha, r, fade = 0, dsize, radius, x2, y2, dx, dy, dz, observer_angle; static double h, h1; static int first = 1; @@ -1034,7 +1158,7 @@ void draw_configuration(t_perco *cell, int *cluster_sizes, int ncells, int nx, i if ((ADD_PLOT)&&(in_plot_box((double)(i+1)*dsize, (double)(j)*dsize))) fade = 1; else fade = 0; - set_cell_color(cell[i*ny+j], cluster_sizes, fade, max_cluster_size); + set_cell_color(cell[i*ny+j], cluster_sizes, fade, max_cluster_size, 0, 1, 0, 1); glVertex2i(i*size, j*size); glVertex2i((i+1)*size, j*size); @@ -1070,7 +1194,7 @@ void draw_configuration(t_perco *cell, int *cluster_sizes, int ncells, int nx, i if ((ADD_PLOT)&&(in_plot_box((double)(i+1)*dsize, (double)(j)*dsize))) fade = 1; else fade = 0; - set_cell_color(cell[i+nx*j], cluster_sizes, fade, max_cluster_size); + set_cell_color(cell[i+nx*j], cluster_sizes, fade, max_cluster_size, 0, 1, 0, 1); glVertex2i(i*size, j*size); glVertex2i((i+1)*size, j*size); @@ -1083,7 +1207,7 @@ void draw_configuration(t_perco *cell, int *cluster_sizes, int ncells, int nx, i if ((ADD_PLOT)&&(in_plot_box((double)(i+1)*dsize, (double)(j+1)*dsize))) fade = 1; else fade = 0; - set_cell_color(cell[ishift+ny*i+j], cluster_sizes, fade, max_cluster_size); + set_cell_color(cell[ishift+ny*i+j], cluster_sizes, fade, max_cluster_size, 0, 1, 0, 1); glVertex2i(i*size, j*size); glVertex2i(i*size, (j+1)*size); @@ -1115,7 +1239,7 @@ void draw_configuration(t_perco *cell, int *cluster_sizes, int ncells, int nx, i if ((ADD_PLOT)&&(in_plot_box(x, y))) fade = 1; else fade = 0; - set_cell_color(cell[j*nx+i], cluster_sizes, fade, max_cluster_size); + set_cell_color(cell[j*nx+i], cluster_sizes, fade, max_cluster_size, 0, 1, 0, 1); glBegin(GL_TRIANGLE_FAN); @@ -1165,7 +1289,7 @@ void draw_configuration(t_perco *cell, int *cluster_sizes, int ncells, int nx, i /* vertical bonds */ if (j=0; i--) if (plot_cube(cell[k*nx*ny+j*nx+i])) + draw_cube_ijk(i, j, k, cell, cluster_sizes, nx, ny, nz, size, max_cluster_size); + break; + } + case (3): + { + for (i=nx-1; i>= 0; i--) + for (j=0; j= 0; i--) + for (j=ny-1; j>=0; j--) if (plot_cube(cell[k*nx*ny+j*nx+i])) + draw_cube_ijk(i, j, k, cell, cluster_sizes, nx, ny, nz, size, max_cluster_size); + break; + } + case (5): + { + for (j=ny-1; j>=0; j--) + for (i=nx-1; i>=0; i--) if (plot_cube(cell[k*nx*ny+j*nx+i])) + draw_cube_ijk(i, j, k, cell, cluster_sizes, nx, ny, nz, size, max_cluster_size); + break; + } + case (6): + { + for (j=ny-1; j>=0; j--) + for (i=0; i=0; j--) if (plot_cube(cell[k*nx*ny+j*nx+i])) + draw_cube_ijk(i, j, k, cell, cluster_sizes, nx, ny, nz, size, max_cluster_size); + break; + } + } + } +// else +// { +// for (i=0; i < nx; i++) +// for (j=0; j 1.0) y1 -= 2.0; + while (y1 < -1.0) y1 += 2.0; + + if (xy_in[i*NY+j]) + { + f = delta*cos(a*xy[0]); + cplus = cosh((y1 + b)/rho); + cminus = cosh((y1 - b)/rho); + + if (y1 > 0.0) + { + phi[1][i*NY+j] = amp*(f + 1.0/(rho*cminus*cminus)); + phi[0][i*NY+j] = amp*(f/(a*a) + rho/(cminus*cminus)); + } + else + { + phi[1][i*NY+j] = amp*(f - 1.0/(rho*cplus*cplus)); + phi[0][i*NY+j] = amp*(f/(a*a) - rho/(cplus*cplus)); + } + } + else + { + phi[0][i*NY+j] = 0.0; + phi[1][i*NY+j] = 0.0; + } + } +} + + /*********************/ /* animation part */ /*********************/ @@ -417,7 +561,7 @@ void compute_gradient_xy(double phi[NX*NY], double gradient[2*NX*NY]) /* compute the gradient of the field */ { int i, j, iplus, iminus, jplus, jminus, padding = 0; - double deltaphi; + double deltaphi, maxgradient = 1.0e10; double dx = (XMAX-XMIN)/((double)NX); dx = (XMAX-XMIN)/((double)NX); @@ -426,12 +570,17 @@ void compute_gradient_xy(double phi[NX*NY], double gradient[2*NX*NY]) for (i=1; i0)||(j!=ropening))&&(maze[n].west)) + { + polyrect[nsides].x1 = x1 - width; + polyrect[nsides].y1 = y1 - width; + polyrect[nsides].x2 = x1 + width; + polyrect[nsides].y2 = y1 + width + dy; + } + nsides++; + + if (maze[n].south) + { + polyrect[nsides].x1 = x1 - width; + polyrect[nsides].y1 = y1 - width; + polyrect[nsides].x2 = x1 + width + dx; + polyrect[nsides].y2 = y1 + width; + } + nsides++; + } + + /* top side of maze */ + polyrect[nsides].x1 = YMIN + padding + MAZE_XSHIFT; + polyrect[nsides].y1 = YMAX - padding - width; + polyrect[nsides].x2 = YMAX - padding + MAZE_XSHIFT; + polyrect[nsides].y2 = YMAX - padding + width; + nsides++; + + /* right side of maze */ + y1 = YMIN + padding + dy*((double)ropening); + x1 = YMAX - padding + MAZE_XSHIFT; + polyrect[nsides].x1 = x1 - width; + polyrect[nsides].y1 = YMIN - 1.0; + polyrect[nsides].x2 = x1 + width; + polyrect[nsides].y2 = y1 - dy; + nsides++; + + polyrect[nsides].x1 = x1 - width; + polyrect[nsides].y1 = y1; + polyrect[nsides].x2 = x1 + width; + polyrect[nsides].y2 = YMAX + 1.0; + nsides++; + + /* left side of maze */ + x1 = YMIN + padding + MAZE_XSHIFT; + polyrect[nsides].x1 = x1 - width; + polyrect[nsides].y1 = YMIN - 1.0; + polyrect[nsides].x2 = x1 + width; + polyrect[nsides].y2 = YMIN + padding; + nsides++; + + polyrect[nsides].x1 = x1 - width; + polyrect[nsides].y1 = YMAX - padding; + polyrect[nsides].x2 = x1 + width; + polyrect[nsides].y2 = YMAX + 1.0; + nsides++; + + if (closed) + { + polyrect[nsides].x1 = XMIN - 0.5*width; + polyrect[nsides].y1 = YMIN - 0.5*width; + polyrect[nsides].x2 = XMIN + 0.5*width; + polyrect[nsides].y2 = YMAX + 0.5*width; + nsides++; + + polyrect[nsides].x1 = XMIN - 0.5*width; + polyrect[nsides].y1 = YMIN - 0.5*width; + polyrect[nsides].x2 = x1 + 0.5*width; + polyrect[nsides].y2 = YMIN + 0.5*width; + nsides++; + + polyrect[nsides].x1 = XMIN - 0.5*width; + polyrect[nsides].y1 = YMAX - 0.5*width; + polyrect[nsides].x2 = x1 + 0.5*width; + polyrect[nsides].y2 = YMAX + 0.5*width; + nsides++; + } + + for (i=0; i x2) return(0); + if (y < y1) return(0); + if (y > y2) return(0); + return(1); +} + + +int ij_in_polyrect(double i, double j, t_rectangle rectangle) +/* returns 1 if (x,y) is in rectangle */ +{ + int i1, i2, j1, j2; + + if (rectangle.posi1 < rectangle.posi2) + { + i1 = rectangle.posi1; + i2 = rectangle.posi2; + } + else + { + i1 = rectangle.posi2; + i2 = rectangle.posi1; + } + if (rectangle.posj1 < rectangle.posj2) + { + j1 = rectangle.posj1; + j2 = rectangle.posj2; + } + else + { + j1 = rectangle.posj2; + j2 = rectangle.posj1; + } + if (i < i1) return(0); + if (i > i2) return(0); + if (j < j1) return(0); + if (j > j2) return(0); + return(1); +} + + int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_circle *circles) /* returns 1 if (x,y) represents a point in the billiard */ { - double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width, alpha, s, a; + 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; @@ -2198,6 +2434,57 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ else if ((x > 0.0)&&(y > 0.0)) return(0); else return(1); } + case (D_WAVEGUIDE): + { + x1 = XMIN + MU; + x2 = XMAX - 2.0*MU - 1.5*LAMBDA; + y1 = 0.5*LAMBDA; + y2 = 1.5*LAMBDA; + if (x < x1) return(0); + if (x > x2 + 1.5*LAMBDA) return(0); + if (vabs(y) > y2) return(0); + if (x < x2) return(vabs(y) >= y1); + r = module2(x-x2, y); + if (r < 0.5*LAMBDA) return(0); + if (r > 1.5*LAMBDA) return(0); + return(1); + } + case (D_WAVEGUIDE_W): + { + x1 = vabs(x); + width = LAMBDA - 2.0*MU; + height = 0.5*MU; + if (x1 > 2.0*LAMBDA - MU) return(0); + if (vabs(y) > MU + width + height) return(0); + if (y >= height) + { + r = module2(x1, y-height); + if ((r > MU)&&(r < MU + width)) return(1); + if (x1 > LAMBDA + MU) return(1); + return(0); + } + if (y <= -height) + { + r = module2(x1-LAMBDA, y+height); + if ((r > MU)&&(r < MU + width)) return(1); + return(0); + } + if (x1 > LAMBDA + MU) return(1); + if ((x1 > MU)&&(x1 < MU + width)) return(1); + return(0); + } + case (D_MAZE): + { + for (i=0; i polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0); + return(1); + } + case (D_MAZE_CLOSED): + { + for (i=0; i polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0); + return(1); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -2390,7 +2677,7 @@ void tvertex_lineto(t_vertex z) void draw_billiard(int fade, double fade_value) /* draws the billiard boundary */ { - double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax; + 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; static int first = 1, nsides; @@ -3209,6 +3496,76 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound draw_line(-LAMBDA + 0.5*MU, YMIN, LAMBDA + 0.5*MU, YMAX); break; } + case (D_WAVEGUIDE): + { + x1 = XMIN + MU; + x2 = XMAX - 2.0*MU - 1.5*LAMBDA; + y1 = 0.5*LAMBDA; + y2 = 1.5*LAMBDA; + + glLineWidth(BOUNDARY_WIDTH); + draw_line(x1, y1, x2, y1); + draw_line(x1, y1, x1, y2); + draw_line(x1, y2, x2, y2); + draw_line(x1, -y1, x2, -y1); + draw_line(x1, -y1, x1, -y2); + draw_line(x1, -y2, x2, -y2); + + dphi = PI/(double)NSEG; + glBegin(GL_LINE_STRIP); + for (i=0; i #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 */ /* 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 1920 /* number of grid points on x axis */ -// // #define NY 1000 /* number of grid points on y axis */ -// // // #define YMIN -1.041666667 -// // // #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ -// +#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 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 HIGHRES 0 /* set to 1 if resolution of grid is double that of displayed image */ +#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 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 640 /* number of grid points on x axis */ -// #define NY 360 /* number of grid points on y axis */ - -#define XMIN -2.0 -#define XMAX 2.0 /* x interval */ -#define YMIN -1.125 -#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ - +// +// // #define 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 JULIA_SCALE 0.8 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 49 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 7 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 201 /* pattern of circles or polygons, see list in global_pdes.c */ +#define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */ #define COMPARISON 0 /* set to 1 to compare two different patterns */ #define B_DOMAIN_B 20 /* second domain shape, for comparisons */ #define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */ #define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ -#define IOR 1 /* choice of index of refraction, see list in global_pdes.c */ +#define IOR 2 /* choice of index of refraction, see list in global_pdes.c */ #define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA -0.5 /* parameter controlling the dimensions of domain */ -#define MU 0.25 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.7 /* 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 */ #define MDEPTH 7 /* depth of computation of Menger gasket */ @@ -110,8 +107,8 @@ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 36 /* number of grid point for grid of disks */ -#define NGRIDY 6 /* number of grid point for grid of disks */ +#define NGRIDX 12 /* number of grid point for grid of disks */ +#define NGRIDY 16 /* number of grid point for grid of disks */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -131,18 +128,19 @@ /* Physical parameters of wave equation */ #define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ -#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */ -#define OSCILLATE_TOPBOT 1 /* set to 1 to enforce a planar wave on top and bottom boundary */ +#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ +#define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ #define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */ #define OMEGA 0.001 /* frequency of periodic excitation */ -#define AMPLITUDE 0.8 /* amplitude of periodic excitation */ +// #define OMEGA 0.0005 /* frequency of periodic excitation */ +#define AMPLITUDE 0.8 /* amplitude of periodic excitation */ #define ACHIRP 0.2 /* acceleration coefficient in chirp */ -#define COURANT 0.1 /* Courant number */ -#define COURANTB 0.025 /* Courant number in medium B */ -#define DAMPING 0.0 /* damping of periodic excitation */ -#define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMAB 0.0 /* damping factor in wave equation */ +#define DAMPING 0.0 /* damping of periodic excitation */ +#define COURANT 0.06 /* Courant number */ +#define COURANTB 0.15 /* Courant number in medium B */ +#define GAMMA 0.0 /* damping factor in wave equation */ +#define GAMMAB 0.015 /* 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 */ @@ -158,14 +156,16 @@ /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 3 +// #define B_COND 1 +#define B_COND 2 #define PRECOMPUTE_BC 0 /* set to 1 to compute neighbours for Laplacian in advance */ /* Parameters for length and speed of simulation */ -#define NSTEPS 3000 /* number of frames of movie */ -#define NVID 8 /* number of iterations between images displayed on screen */ +// #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 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,9 +182,13 @@ /* 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.0005 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ +#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */ + /* Plot type, see list in global_pdes.c */ @@ -199,13 +203,13 @@ #define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */ #define FLOOR_ZCOORD 1 /* set to 1 to draw only facets with z not too negative */ #define DRAW_BILLIARD 1 /* set to 1 to draw boundary */ -#define DRAW_BILLIARD_FRONT 0 /* set to 1 to draw front of boundary after drawing wave */ +#define DRAW_BILLIARD_FRONT 1 /* set to 1 to draw front of boundary after drawing wave */ #define DRAW_CONSTRUCTION_LINES 0 /* set to 1 to draw construction lines of certain domains */ #define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */ #define DRAW_OUTSIDE_GRAY 0 /* experimental, draw outside of billiard in gray */ -#define PLOT_SCALE_ENERGY 0.5 /* vertical scaling in energy plot */ -#define PLOT_SCALE_LOG_ENERGY 1.0 /* vertical scaling in log energy plot */ +#define PLOT_SCALE_ENERGY 0.01 /* vertical scaling in energy plot */ +#define PLOT_SCALE_LOG_ENERGY 0.25 /* vertical scaling in log energy plot */ /* 3D representation */ @@ -219,24 +223,24 @@ /* Color schemes */ -#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ +#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 BLACK 1 /* background */ #define COLOR_SCHEME 3 /* choice of color scheme, see list in global_pdes.c */ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ -#define SLOPE 0.5 /* sensitivity of color on wave amplitude */ +#define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ -#define VSCALE_ENERGY 60.0 /* additional scaling factor for color scheme P_3D_ENERGY */ +#define VSCALE_ENERGY 100.0 /* additional scaling factor for color scheme P_3D_ENERGY */ #define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 250.0 /* scaling factor for energy representation */ +#define E_SCALE 100.0 /* scaling factor for energy representation */ #define LOG_SCALE 0.5 /* scaling factor for energy log representation */ -#define LOG_SHIFT 0.0 /* shift of colors on log scale */ -#define LOG_ENERGY_FLOOR -100.0 /* floor value for log of (total) energy */ +#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 RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ @@ -247,13 +251,28 @@ #define HUEMEAN 240.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -200.0 /* amplitude of variation of hue for color scheme C_HUE */ +#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 */ +// 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 3.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 15.0 /* scale of color scheme bar for 2nd part */ +#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 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 */ +/* for compatibility with sub_wave and sub_maze */ +#define ADD_POTENTIAL 0 +// #define POT_MAZE 7 +#define POTENTIAL 0 +/* end of constants only used by sub_wave and sub_maze */ + + /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 10.0 /* max value of wave amplitude */ @@ -264,21 +283,22 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {8.0, 8.0, 6.0}; /* location of observer for REP_PROJ_3D representation */ +double observer[3] = {8.0, 8.0, 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.2 /* overall scaling factor of z axis for REP_PROJ_3D representation */ -#define XY_SCALING_FACTOR 2.1 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define Z_SCALING_FACTOR 0.35 /* 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 */ #define YSHIFT_3D 0.1 /* overall y shift for REP_PROJ_3D representation */ #include "global_pdes.c" /* constants and global variables */ +#include "global_3d.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 */ #include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ -#include "global_3d.c" /* constants and global variables */ #include "sub_wave_3d.c" /* graphical functions specific to wave_3d */ FILE *time_series_left, *time_series_right; @@ -961,6 +981,9 @@ void animation() npolyline = init_polyline(MDEPTH, polyline); for (i=0; i #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 */ /* General geometrical parameters */ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1000 /* window height */ // #define NX 1920 /* number of grid points on x axis */ -// #define NY 1000 /* number of grid points on y axis */ +// #define NY 1000 /* number of grid points on x axis */ #define NX 3840 /* number of grid points on x axis */ #define NY 2000 /* number of grid points on y axis */ @@ -65,6 +65,8 @@ // #define WINWIDTH 1280 /* window width */ // #define WINHEIGHT 720 /* window height */ // +// // #define NX 640 /* number of grid points on x axis */ +// // #define NY 360 /* number of grid points on y axis */ // // #define NX 1280 /* number of grid points on x axis */ // // #define NY 720 /* number of grid points on y axis */ // #define NX 2560 /* number of grid points on x axis */ @@ -79,11 +81,11 @@ /* Choice of the billiard table */ -#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 53 /* 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 */ -#define COMPARISON 0 /* set to 1 to compare two different patterns */ +#define COMPARISON 0 /* set to 1 to compare two different patterns (beta) */ #define B_DOMAIN_B 20 /* second domain shape, for comparisons */ #define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */ @@ -91,10 +93,10 @@ #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 MU 0.028 /* parameter controlling the dimensions of domain */ -#define NPOLY 3 /* number of sides of polygon */ -#define APOLY 0.3333333333333333 /* angle by which to turn polygon, in units of Pi/2 */ +#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ +#define MU 0.3 /* 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 */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ @@ -120,18 +122,18 @@ /* Physical parameters of wave equation */ #define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ -#define OSCILLATE_LEFT 1 /* set to 1 to add oscilating boundary condition on the left */ +#define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ -#define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */ +#define OSCILLATION_SCHEDULE 1 /* oscillation schedule, see list in global_pdes.c */ #define OMEGA 0.0005 /* frequency of periodic excitation */ #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.05 /* Courant number */ -#define COURANTB 0.0125 /* Courant number in medium B */ +#define COURANT 0.1 /* Courant number */ +#define COURANTB 0.05 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ -#define GAMMAB 0.0 /* damping factor in wave equation */ +#define GAMMAB 5.0e-3 /* 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 */ @@ -142,24 +144,26 @@ /* 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 6 /* period of oscillating source */ +#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ +#define OSCILLATING_SOURCE_PERIOD 50 /* period of oscillating source */ /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 3 +#define B_COND 2 /* Parameters for length and speed of simulation */ -#define NSTEPS 2700 /* number of frames of movie */ -#define NVID 15 /* number of iterations between images displayed on screen */ +// #define NSTEPS 1000 /* number of frames of movie */ +#define NSTEPS 2400 /* 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 50 /* time after which to start saving frames */ -#define BOUNDARY_WIDTH 1 /* width of billiard boundary */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ +// #define INITIAL_TIME 400 /* time after which to start saving frames */ +#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PRINT_SPEED 0 /* print speed of moving source */ #define PAUSE 200 /* number of frames after which to pause */ -#define PSLEEP 2 /* sleep time during pause */ +#define PSLEEP 1 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ #define MID_FRAMES 20 /* number of still frames between parts of two-part movie */ @@ -168,20 +172,22 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 0.75 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.00005 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.004 /* wavelength of initial condition */ +// #define INITIAL_AMP 0.75 /* amplitude of initial condition */ +#define INITIAL_AMP 1.0 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -#define PLOT 1 +#define PLOT 0 +// #define PLOT 1 -#define PLOT_B 0 /* plot type for second movie */ +#define PLOT_B 5 /* plot type for second movie */ /* Color schemes */ -#define COLOR_PALETTE 13 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */ +#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 BLACK 1 /* background */ @@ -192,7 +198,7 @@ #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 250.0 /* scaling factor for energy representation */ +#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 RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ @@ -205,17 +211,31 @@ #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 5.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 2.0 /* scale of color scheme bar for 2nd part */ +#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 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 */ +#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 MAZE_XSHIFT 0.0 /* horizontal shift of maze */ + +/* for compatibility with sub_wave and sub_maze */ +#define ADD_POTENTIAL 0 +#define POT_MAZE 7 +#define POTENTIAL 0 +/* end of constants only used by sub_wave and sub_maze */ + + /* For debugging purposes only */ #define FLOOR 0 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 10.0 /* max value of wave amplitude */ #include "global_pdes.c" /* constants and global variables */ +#include "sub_maze.c" /* support for generating mazes */ #include "sub_wave.c" /* common functions for wave_billiard, heat and schrodinger */ #include "wave_common.c" /* common functions for wave_billiard, wave_comparison, etc */ @@ -548,7 +568,6 @@ 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], *phi_tmp[NX], *psi_tmp[NX], *total_energy[NX], *color_scale[NX]; double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX]; short int *xy_in[NX]; int i, j, s, sample_left[2], sample_right[2], period = 0, fade; @@ -566,8 +585,6 @@ void animation() { phi[i] = (double *)malloc(NY*sizeof(double)); psi[i] = (double *)malloc(NY*sizeof(double)); -// phi_tmp[i] = (double *)malloc(NY*sizeof(double)); -// psi_tmp[i] = (double *)malloc(NY*sizeof(double)); tmp[i] = (double *)malloc(NY*sizeof(double)); total_energy[i] = (double *)malloc(NY*sizeof(double)); xy_in[i] = (short int *)malloc(NY*sizeof(short int)); @@ -582,6 +599,9 @@ void animation() /* initialise polyline for von Koch and similar domains */ npolyline = init_polyline(MDEPTH, polyline); for (i=0; i