diff --git a/drop_billiard.c b/drop_billiard.c index 830940a..53ba3ad 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -29,6 +29,7 @@ #include /* Sam Leffler's libtiff library. */ #define MOVIE 0 /* set to 1 to generate movie */ +#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ @@ -104,6 +105,8 @@ #define NCOLORS 16 /* number of colors */ #define COLORSHIFT 3 /* hue of initial color */ +#define COLOR_HUEMIN 0 /* minimal color hue */ +#define COLOR_HUEMAX 360 /* maximal color hue */ #define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ #define NSEG 100 /* number of segments of boundary */ #define BILLIARD_WIDTH 4 /* width of billiard */ @@ -127,6 +130,8 @@ #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 MAZE_RANDOM_FACTOR 0.1 /* randomization factor for S_MAZE_RANDOM */ +#define MAZE_CORNER_RADIUS 0.5 /* radius of tounded corners in maze */ diff --git a/global_3d.c b/global_3d.c index 68de88e..859b3bf 100644 --- a/global_3d.c +++ b/global_3d.c @@ -3,7 +3,7 @@ /* plot types used by wave_3d */ #define P_3D_AMPLITUDE 101 /* height/color depends on amplitude - DEPRECATED, instead use set SHADE_3D to 0 */ -#define P_3D_ANGLE 102 /* height/color depends on angle with fixed direction - TO DO */ +#define P_3D_ANGLE 102 /* height/color depends on angle with fixed direction - TODO */ #define P_3D_AMP_ANGLE 103 /* height/color depends on amplitude, luminosity depends on angle */ #define P_3D_ENERGY 104 /* height/color depends on energy, luminosity depends on angle */ #define P_3D_LOG_ENERGY 105 /* height/color depends on logarithm of energy, luminosity depends on angle */ @@ -13,6 +13,8 @@ #define P_3D_LOG_MEAN_ENERGY 109 /* height/color depends on log on energy averaged over time, luminosity depends on angle */ #define P_3D_PHASE 111 /* phase of wave */ +#define P_3D_FLUX_INTENSITY 112 /* energy flux intensity */ +#define P_3D_FLUX_DIRECTION 113 /* energy flux direction */ /* Choice of simulated reaction-diffusion equation in rde.c */ @@ -23,7 +25,8 @@ #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 */ +#define E_EULER_INCOMP 6 /* incompressible Euler equation */ +#define E_EULER_COMP 7 /* compressible Euler equation */ /* Choice of potential */ @@ -40,40 +43,17 @@ #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 */ +/* Choice of force field in compressible Euler equation */ -#define Z_AMPLITUDE 0 /* amplitude of first field */ -#define Z_RGB 20 /* RGB plot */ -#define Z_POLAR 21 /* polar angle associated with RBG plot */ -#define Z_NORM_GRADIENT 22 /* gradient of polar angle */ -#define Z_ANGLE_GRADIENT 221 /* direction of polar angle */ -#define Z_NORM_GRADIENTX 23 /* norm of gradient of u */ -#define Z_ANGLE_GRADIENTX 231 /* direction of gradient of u */ -#define Z_NORM_GRADIENT_INTENSITY 24 /* gradient and intensity of polar angle */ -#define Z_VORTICITY 25 /* curl of polar angle */ -#define Z_VORTICITY_ABS 251 /* absolute value of curl of polar angle */ - -/* for Schrodinger equation */ -#define Z_MODULE 30 /* module squared of first two fields */ -#define Z_ARGUMENT 31 /* argument of first two fields, with luminosity depending on module */ -#define Z_REALPART 32 /* first field, with luminosity depending on module */ - -/* for RPSLZ equation */ -#define Z_MAXTYPE_RPSLZ 40 /* color of type with maximal density */ -#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 */ +#define GF_VERTICAL 0 /* gravity acting vertically */ +#define GF_CIRCLE 1 /* repelling circle */ /* 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)) #define COMPUTE_THETAZ ((zplot == Z_POLAR)||(zplot == Z_NORM_GRADIENT)||(zplot == Z_ANGLE_GRADIENT)||(zplot == Z_NORM_GRADIENT_INTENSITY)||(zplot == Z_VORTICITY)||(zplot == Z_VORTICITY_ABS)) -#define COMPUTE_ENERGY ((zplot == P_3D_ENERGY)||(cplot == P_3D_ENERGY)||(zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_ENERGY)||(zplot == P_3D_TOTAL_ENERGY)||(cplot == P_3D_TOTAL_ENERGY)||(zplot == P_3D_LOG_TOTAL_ENERGY)||(cplot == P_3D_LOG_TOTAL_ENERGY)||(zplot == P_3D_MEAN_ENERGY)||(cplot == P_3D_MEAN_ENERGY)||(zplot == P_3D_LOG_MEAN_ENERGY)||(cplot == P_3D_LOG_MEAN_ENERGY)) +#define COMPUTE_ENERGY ((zplot == P_3D_ENERGY)||(cplot == P_3D_ENERGY)||(zplot == P_3D_LOG_ENERGY)||(cplot == P_3D_LOG_ENERGY)||(zplot == P_3D_TOTAL_ENERGY)||(cplot == P_3D_TOTAL_ENERGY)||(zplot == P_3D_LOG_TOTAL_ENERGY)||(cplot == P_3D_LOG_TOTAL_ENERGY)||(zplot == P_3D_MEAN_ENERGY)||(cplot == P_3D_MEAN_ENERGY)||(zplot == P_3D_LOG_MEAN_ENERGY)||(cplot == P_3D_LOG_MEAN_ENERGY)||(ZPLOT == P_3D_FLUX_INTENSITY)||(CPLOT == P_3D_FLUX_INTENSITY)||(ZPLOT_B == P_3D_FLUX_INTENSITY)||(CPLOT_B == P_3D_FLUX_INTENSITY)||(ZPLOT == P_3D_FLUX_DIRECTION)||(CPLOT == P_3D_FLUX_DIRECTION)||(ZPLOT_B == P_3D_FLUX_DIRECTION)||(CPLOT_B == P_3D_FLUX_DIRECTION)) #define COMPUTE_LOG_TOTAL_ENERGY ((ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)) @@ -83,6 +63,8 @@ #define COMPUTE_MEAN_ENERGY ((ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY)) +#define COMPUTE_ENERGY_FLUX ((ZPLOT == P_3D_FLUX_INTENSITY)||(CPLOT == P_3D_FLUX_INTENSITY)||(ZPLOT_B == P_3D_FLUX_INTENSITY)||(CPLOT_B == P_3D_FLUX_INTENSITY)||(ZPLOT == P_3D_FLUX_DIRECTION)||(CPLOT == P_3D_FLUX_DIRECTION)||(ZPLOT_B == P_3D_FLUX_DIRECTION)||(CPLOT_B == P_3D_FLUX_DIRECTION)) + #define COMPUTE_TOTAL_ENERGY ((ZPLOT == P_3D_TOTAL_ENERGY)||(CPLOT == P_3D_TOTAL_ENERGY)||(ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_TOTAL_ENERGY)||(CPLOT_B == P_3D_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY)) @@ -101,9 +83,14 @@ typedef struct double mean_energy; /* energy averaged since beginning of simulation */ double log_mean_energy; /* log of energy averaged since beginning of simulation */ double cos_angle; /* cos of angle between normal vector and direction of light */ + double flux_intensity; /* intensity of energy flux */ + double flux_direction; /* direction of energy flux */ + double flux_int_table[FLUX_WINDOW]; /* table of energy flux intensities (for averaging) */ + short int flux_counter; /* counter for averaging of energy flux */ 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) */ + double *p_cfield[4]; /* pointers to color field (second pointer for option DOUBLE_MOVIE) */ + /* third and fourth pointer for color luminosity (for energy flux) */ } t_wave; @@ -117,6 +104,8 @@ typedef struct 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 Lpressure; /* Laplacian of pressure (for Euler equation) */ + double dxu, dyu, dxv, dyv; /* gradient of velocity field (for compressible 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 09d71e0..8366ffc 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -42,6 +42,7 @@ #define O_GALTON_BOARD 1 /* Galton board pattern */ #define O_GENUS_TWO 2 /* obstacles in corners of L-shape domeain (for genus 2 b.c.) */ #define O_POOL_TABLE 3 /* obstacles around pockets of pool table */ +#define O_HLINE_HOLE_SPOKES 181 /* tips of spokes for S_HLINE_HOLE_SPOKES segment pattern */ /* pattern of additional repelling segments */ #define S_RECTANGLE 0 /* segments forming a rectangle */ @@ -65,7 +66,9 @@ #define S_EXT_RECTANGLE 16 /* particles outside a rectangle */ #define S_DAM_BRICKS 17 /* dam made of several bricks */ #define S_HLINE_HOLE 18 /* horizontal line with a hole in the bottom */ +#define S_HLINE_HOLE_SPOKES 181 /* horizontal line with a hole in the bottom and extra spokes */ #define S_EXT_CIRCLE_RECT 19 /* particles outside a circle and a rectangle */ +#define S_BIN_OPENING 20 /* bin containing particles opening at deactivation time */ /* particle interaction */ @@ -124,6 +127,27 @@ #define NZ_BROAD 5 /* broad straight nozzle */ #define NZ_NONE 99 /* no nozzle */ +/* Types of chemical reactions */ + +#define CHEM_RPS 0 /* rock-paper-scissors reaction */ +#define CHEM_AAB 1 /* reaction A + A -> B */ +#define CHEM_ABC 2 /* reaction A + B -> C */ +#define CHEM_A2BC 3 /* reaction 2A + B -> C */ +#define CHEM_CATALYSIS 4 /* reaction 2A + C -> B + C */ +#define CHEM_BAA 5 /* reaction B -> A + A (dissociation) */ +#define CHEM_AABAA 6 /* reaction A + A <-> B (reversible) */ +#define CHEM_POLYMER 7 /* reaction A + B -> C, A + C -> D, etc */ +#define CHEM_POLYMER_DISS 8 /* polimerisation with dissociation */ + +/* Initial conditions for chemical reactions */ + +#define IC_UNIFORM 0 /* all particles have type 1 */ +#define IC_UNIFORM2 20 /* all particles have type 2 */ +#define IC_RANDOM_UNIF 1 /* particle type chosen uniformly at random */ +#define IC_RANDOM_TWO 2 /* particle type chosen randomly between 1 and 2, with TYPE_PROPORTION */ +#define IC_CIRCLE 3 /* type 1 in a disc */ +#define IC_CATALYSIS 4 /* mix of 1 and 2 in left half, only 1 in right half */ + /* Plot types */ #define P_KINETIC 0 /* colors represent kinetic energy of particles */ @@ -194,6 +218,7 @@ typedef struct double spin_range; /* range of spin-spin interaction */ double spin_freq; /* angular frequency of spin-spin interaction */ double color_hue; /* color hue of particle, for P_INITIAL_POS plot type */ + int color_rgb[3]; /* RGB colors code of particle, for use in ljones_movie.c */ } t_particle; typedef struct @@ -271,5 +296,12 @@ typedef struct } t_group_data; +typedef struct +{ + double x, y; /* location of collision */ + int time; /* time since collision */ + int color; /* color hue in case of different collisions */ +} t_collision; + int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0; diff --git a/global_particles.c b/global_particles.c index 642a416..010340c 100644 --- a/global_particles.c +++ b/global_particles.c @@ -4,6 +4,7 @@ // int circlecolor[NMAXCIRCLES]; /* color of circular scatterer */ int ncircles = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */ int nsides = NMAXPOLY; /* actual number of sides of polygonal line */ +int narcs = NMAXCIRCLES; /* actual number of arcs */ typedef struct { @@ -18,8 +19,20 @@ typedef struct int color; } t_segment; +typedef struct +{ + double xc, yc, radius, angle1, dangle; + int color; +} t_arc; + +typedef struct +{ + short int left, right; +} t_exit; + t_circle circles[NMAXCIRCLES]; t_segment polyline[NMAXPOLY]; +t_arc arcs[NMAXCIRCLES]; double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; /* shooter and target positions for "laser in room of mirrors" simulations, with default values for square domain */ @@ -67,7 +80,8 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define C_LASER 11 /* laser fight in a room of mirrors */ #define C_LASER_GENUSN 12 /* laser fight in a translation surface */ -#define D_POLYLINE 30 /* general polygon */ +#define D_POLYLINE 30 /* polygonal line */ +#define D_POLYLINE_ARCS 31 /* polygonal line and circular arcs */ #define P_RECTANGLE 0 /* rectangle (for test purposes) */ #define P_TOKARSKY 1 /* Tokarsky unilluminable room */ @@ -80,6 +94,10 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define P_TOKA_NONSELF 8 /* Tokarsky non-self-unilluminable room */ #define P_MAZE 10 /* maze */ #define P_MAZE_DIAG 11 /* maze with 45 degrees angles */ +#define P_MAZE_RANDOM 12 /* maze with randomized wall positions */ +#define P_MAZE_CIRCULAR 13 /* circular maze */ +#define P_MAZE_CIRC_SCATTERER 14 /* circular maze with scatterers */ +#define P_MAZE_HEX 15 /* hexagonal maze */ /* Color palettes */ diff --git a/global_pdes.c b/global_pdes.c index a453062..fde056a 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -68,13 +68,17 @@ #define D_WAVEGUIDE_W 52 /* W-shaped wave guide */ #define D_MAZE 53 /* maze */ #define D_MAZE_CLOSED 54 /* closed maze */ +#define D_MAZE_CHANNELS 541 /* maze with two channels attached */ #define D_CHESSBOARD 55 /* chess board configuration */ #define D_TRIANGLE_TILES 56 /* triangular tiling */ #define D_HEX_TILES 57 /* honeycomb tiling */ +#define D_FUNNELS 58 /* two funnels */ +#define D_ONE_FUNNEL 581 /* one funnel */ +#define D_LENSES_RING 59 /* several lenses forming a ring */ +#define D_MAZE_CIRCULAR 60 /* circular 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) */ -// #define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -114,6 +118,7 @@ #define IOR_MANDELBROT 1 /* index of refraction depends on escape time in Mandelbrot set (log) */ #define IOR_MANDELBROT_LIN 100 /* index of refraction depends on escape time in Mandelbrot set (linear) */ #define IOR_EARTH 2 /* index of refraction models speed of seismic waves */ +#define IOR_EXPLO_LENSING 3 /* explosive lensing */ /* Boundary conditions */ @@ -185,6 +190,47 @@ #define NPWIDTH 0.02 /* width of noise panel separation */ +/* plot types used by rde */ + +#define Z_AMPLITUDE 0 /* amplitude of first field */ +#define Z_RGB 20 /* RGB plot */ +#define Z_POLAR 21 /* polar angle associated with RBG plot */ +#define Z_NORM_GRADIENT 22 /* gradient of polar angle */ +#define Z_ANGLE_GRADIENT 221 /* direction of polar angle */ +#define Z_NORM_GRADIENTX 23 /* norm of gradient of u */ +#define Z_ANGLE_GRADIENTX 231 /* direction of gradient of u */ +#define Z_NORM_GRADIENT_INTENSITY 24 /* gradient and intensity of polar angle */ +#define Z_VORTICITY 25 /* curl of polar angle */ +#define Z_VORTICITY_ABS 251 /* absolute value of curl of polar angle */ + +/* for Schrodinger equation */ +#define Z_MODULE 30 /* module squared of first two fields */ +#define Z_ARGUMENT 31 /* argument of first two fields, with luminosity depending on module */ +#define Z_REALPART 32 /* first field, with luminosity depending on module */ + +/* for RPSLZ equation */ +#define Z_MAXTYPE_RPSLZ 40 /* color of type with maximal density */ +#define Z_THETA_RPSLZ 41 /* polar angle */ +#define Z_NORM_GRADIENT_RPSLZ 42 /* gradient of polar angle */ + +/* for Euler incompressible 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 */ +#define Z_EULER_LPRESSURE 53 /* Laplacian of pressure */ +#define Z_EULER_PRESSURE 54 /* pressure */ + +/* for Euler compressible Euler equation */ +#define Z_EULER_DENSITY 60 /* density */ +#define Z_EULER_SPEED 61 /* norm of velocity */ +#define Z_EULERC_VORTICITY 62 /* vorticity of velocity */ + +/* special boundary conditions for Euler equation */ +#define BCE_TOPBOTTOM 1 /* laminar flow at top and bottom */ +#define BCE_TOPBOTTOMLEFT 2 /* laminar flow at top, bottom and left side */ +#define BCE_CHANNELS 3 /* laminar flow in channels at left and right */ +#define BCE_MIDDLE_STRIP 4 /* laminar flow in horizontal strip in the middle */ + typedef struct { double xc, yc, radius; /* center and radius of circle */ @@ -210,6 +256,23 @@ typedef struct double posi1, posj1, posi2, posj2; /* (i,j) coordinates of vertices */ } t_rectangle; +typedef struct +{ + double x1, y1, x2, y2; /* (x,y) coordinates of long symmetry axis */ + double width; /* width of rectangle */ + double posi1, posj1, posi2, posj2; /* (i,j) coordinates of vertices */ + double posi3, posj3, posi4, posj4; /* (i,j) coordinates of vertices */ +} t_rect_rotated; + +typedef struct +{ + double xc, yc; /* (x,y) coordinates of center */ + double r, width; /* radius and width of arc */ + double angle1, dangle; /* start angle and angular width */ + double posi1, posj1; /* (i,j) coordinates of center */ +// double posi3, posj3, posi4, posj4; /* (i,j) coordinates of vertices */ +} t_arc; + typedef struct { int nneighb; /* number of neighbours to compute Laplacian */ @@ -220,11 +283,15 @@ 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 */ +int npolyrect_rot = NMAXPOLY; /* actual number of rotated polyrect */ +int npolyarc = NMAXPOLY; /* actual number of arcs */ 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 */ +t_rect_rotated polyrectrot[NMAXPOLY]; /* data of rotated rectangles */ +t_arc polyarc[NMAXPOLY]; /* data of arcs */ /* 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/heat.c b/heat.c index 299cd6c..dd81957 100644 --- a/heat.c +++ b/heat.c @@ -76,7 +76,7 @@ #define MU 0.8 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ -#define MDEPTH 6 /* depth of computation of Menger gasket */ +#define MDEPTH 1 /* depth of computation of Menger gasket */ #define MRATIO 5 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ diff --git a/lennardjones.c b/lennardjones.c index 57b0e49..b4aaa4d 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -37,12 +37,16 @@ #include #define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define SAVE_MEMORY 1 /* set to 1 to save memory while saving frames */ +#define NO_EXTRA_BUFFER_SWAP 0 /* some OS require one less buffer swap when recording images */ -#define TIME_LAPSE 0 /* set to 1 to add a time-lapse movie at the end */ +#define TIME_LAPSE 1 /* set to 1 to add a time-lapse movie at the end */ /* so far incompatible with double movie */ #define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */ -#define TIME_LAPSE_FIRST 0 /* set to 1 to show time-lapse version first */ +#define TIME_LAPSE_FIRST 1 /* set to 1 to show time-lapse version first */ + +#define SAVE_TIME_SERIES 0 /* set to 1 to save time series of particle positions */ /* General geometrical parameters */ @@ -54,10 +58,10 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.95 -#define INITXMAX 1.95 /* x interval for initial condition */ -#define INITYMIN -1.1 -#define INITYMAX 0.15 /* y interval for initial condition */ +#define INITXMIN -1.92 +#define INITXMAX 1.92 /* x interval for initial condition */ +#define INITYMIN -1.05 +#define INITYMAX 1.05 /* y interval for initial condition */ #define BCXMIN -2.0 #define BCXMAX 2.0 /* x interval for boundary condition */ @@ -67,20 +71,20 @@ #define OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ -#define CIRCLE_PATTERN 8 /* 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 OBSTACLE_PATTERN 181 /* 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 19 /* pattern of repelling segments, see list in global_ljones.c */ +#define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */ +#define SEGMENT_PATTERN 181 /* 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 0 /* shape of nozzle for second rocket, see list in global_ljones.c */ +#define NOZZLE_SHAPE 2 /* shape of nozzle, see list in global_ljones.c */ +#define NOZZLE_SHAPE_B 4 /* 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.5 /* proportion of particles of first type */ +#define TYPE_PROPORTION 0.66 /* 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 */ @@ -93,16 +97,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 3.0 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 4.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.3 /* parameter controlling the dimensions of domain */ -#define MU 0.009 /* parameter controlling radius of particles */ -#define MU_B 0.018 /* parameter controlling radius of particles of second type */ -// #define NPOLY 4 /* number of sides of polygon */ -#define NPOLY 3 /* number of sides of polygon */ -#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ +#define MU 0.008 /* parameter controlling radius of particles */ +#define MU_B 0.012 /* parameter controlling radius of particles of second type */ +#define NPOLY 25 /* number of sides of polygon */ +#define APOLY 0.666666666 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ @@ -122,11 +125,11 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1500 /* number of frames of movie */ -#define NVID 100 /* number of iterations between images displayed on screen */ +#define NSTEPS 1000 /* number of frames of movie */ +#define NVID 150 /* number of iterations between images displayed on screen */ #define NSEG 250 /* number of segments of boundary */ -#define INITIAL_TIME 250 /* time after which to start saving frames */ -#define OBSTACLE_INITIAL_TIME 350 /* time after which to start moving obstacle */ +#define INITIAL_TIME 25 /* time after which to start saving frames */ +#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 */ @@ -140,18 +143,18 @@ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 0 +#define BOUNDARY_COND 3 /* Plot type, see list in global_ljones.c */ -#define PLOT 0 -#define PLOT_B 8 /* plot type for second movie */ +#define PLOT 5 +#define PLOT_B 0 /* plot type for second movie */ -#define DRAW_BONDS 0 /* set to 1 to draw bonds between neighbours */ +#define DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ #define FILL_TRIANGLES 1 /* set to 1 to fill triangles between neighbours */ #define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */ -#define COLOR_SEG_GROUPS 1 /* set to 1 to collor segment groups differently */ +#define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */ /* Color schemes */ @@ -178,36 +181,40 @@ #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 3.0e2 /* energy of particle with hottest color */ +#define PARTICLE_EMAX 1.2e3 /* 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 HUE_TYPE1 150.0 /* hue of particles of type 1 */ +#define HUE_TYPE2 190.0 /* hue of particles of type 2 */ +#define HUE_TYPE3 220.0 /* hue of particles of type 3 */ +// #define HUE_TYPE1 310.0 /* hue of particles of type 1 */ +// #define HUE_TYPE2 150.0 /* hue of particles of type 2 */ +// #define HUE_TYPE3 180.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 4.0 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST 3.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 INITIAL_DAMPING 35.0 /* damping coefficient of particles during initial phase */ +#define REPEL_RADIUS 15.0 /* radius in which repelling force acts (in units of particle radius) */ +#define DAMPING 0.0 /* damping coefficient of particles */ +#define INITIAL_DAMPING 0.0 /* damping coefficient of particles during initial phase */ #define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ -#define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU */ -#define PARTICLE_INERTIA_MOMENT 0.2 /* moment of inertia of particle */ +#define PARTICLE_MASS_B 2.0 /* mass of particle of radius MU */ +#define PARTICLE_INERTIA_MOMENT 0.02 /* 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.01 /* initial inverse temperature */ #define MU_XI 0.01 /* friction constant in thermostat */ #define KSPRING_BOUNDARY 1.0e7 /* confining harmonic potential outside simulation region */ #define KSPRING_OBSTACLE 1.0e11 /* harmonic potential of obstacles */ -#define NBH_DIST_FACTOR 7.5 /* radius in which to count neighbours */ -#define GRAVITY 1500.0 /* gravity acting on all particles */ +#define NBH_DIST_FACTOR 10.0 /* radius in which to count neighbours */ +// #define NBH_DIST_FACTOR 7.5 /* radius in which to count neighbours */ +#define GRAVITY 0.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 */ @@ -215,7 +222,7 @@ #define GRAVITY_INITIAL_TIME 200 /* time at start of simulation with constant gravity */ #define GRAVITY_RESTORE_TIME 700 /* time at end of simulation with gravity restored to initial value */ -#define ROTATION 0 /* set to 1 to include rotation of particles */ +#define ROTATION 1 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ #define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */ #define KTORQUE 100.0 /* force constant in angular dynamics */ @@ -232,8 +239,8 @@ #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 1600 /* final phase in which temperature is constant */ -#define FINAL_DECREASE_PHASE 1500 /* final phase in which temperature decreases */ +#define MIDDLE_CONSTANT_PHASE 2000 /* final phase in which temperature is constant */ +#define FINAL_DECREASE_PHASE 1300 /* 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 */ @@ -271,6 +278,7 @@ #define ADD_PARTICLES 0 /* set to 1 to add particles */ #define ADD_TIME 0 /* time at which to add first particle */ #define ADD_PERIOD 10000 /* time interval between adding further particles */ +#define N_ADD_PARTICLES 20 /* number of particles to add */ #define FINAL_NOADD_PERIOD 200 /* final period where no particles are added */ #define SAFETY_FACTOR 2.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ @@ -289,40 +297,48 @@ #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 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.0 /* initial position of segments */ -#define SEGMENTS_Y0 0.8 /* initial position of segments */ +#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 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 750.0 /* mass of segment group */ -#define SEGMENT_GROUP_I 500.0 /* moment of inertia of segment group */ +#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 0 /* set to 1 to add repulsion between centers of mass 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 0 /* set to 1 for view to track group of segments */ +#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 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) */ -#define RD_TYPES 3 /* number of types in reaction-diffusion equation */ -#define REACTION_PROB 0.0045 /* probability controlling reaction term */ +#define REACTION_DIFFUSION 1 /* set to 1 to simulate a chemical reaction (particles may change type) */ +#define RD_REACTION 8 /* type of reaction, see list in global_ljones.c */ +#define RD_TYPES 9 /* number of types in reaction-diffusion equation */ +#define RD_INITIAL_COND 2 /* initial condition of particles */ +#define REACION_DIST 3.5 /* maximal distance for reaction to occur */ +#define REACTION_PROB 0.5 /* probability controlling reaction term */ +#define DISSOCIATION_PROB 0.005 /* probability controlling dissociation reaction */ +#define CENTER_COLLIDED_PARTICLES 0 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */ +#define COLLISION_TIME 25 /* time during which collisions are shown */ #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ +#define PLOT_PARTICLE_NUMBER 1 /* set to 1 to make of plot of particle number over time */ +#define PARTICLE_NB_PLOT_FACTOR 1.0 /* expected final number of particles over initial number */ #define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ #define PLOT_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */ #define PLOT_TRAJECTORIES 0 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */ @@ -345,12 +361,12 @@ #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 FMAX 1.0e10 /* maximal force */ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 120 /* size of hashgrid in x direction */ -#define HASHY 60 /* size of hashgrid in y direction */ +#define HASHX 90 /* size of hashgrid in x direction */ +#define HASHY 45 /* 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 */ @@ -385,6 +401,7 @@ int thermostat_on = 1; /* thermostat switch used when VARY_THERMOSTAT is on * #include "sub_lj.c" #include "sub_hashgrid.c" +FILE *lj_time_series, *lj_final_position; /*********************/ /* animation part */ @@ -1032,12 +1049,14 @@ 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, speed_ratio; + fboundary = 0.0, pleft = 0.0, pright = 0.0, entropy[2], mean_energy, gravity = GRAVITY, speed_ratio, + ymin, ymax; 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, - group, gshift; + group, gshift, n_total_active = 0, ncollisions = 0; + int *particle_numbers; static int imin, imax; static short int first = 1; t_particle *particle; @@ -1046,6 +1065,7 @@ void animation() t_group_segments *segment_group; t_tracer *trajectory; t_group_data *group_speeds; + t_collision *collisions; t_hashgrid *hashgrid; char message[100]; @@ -1058,7 +1078,8 @@ void animation() 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)); + 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 */ @@ -1070,6 +1091,17 @@ void animation() pangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); pressure = (double *)malloc(N_PRESSURES*sizeof(double)); + if (REACTION_DIFFUSION) + { + collisions = (t_collision *)malloc(2*NMAXCIRCLES*sizeof(t_collision)); + for (i=0; i<2*NMAXCIRCLES; i++) collisions[i].time = 0; + } + + if (SAVE_TIME_SERIES) + { + lj_time_series = fopen("lj_time_series.dat", "w"); + lj_final_position = fopen("lj_final_position.dat", "w"); + } /* initialise positions and radii of circles */ init_particle_config(particle); @@ -1081,7 +1113,7 @@ void animation() if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle); if (ADD_FIXED_SEGMENTS) init_segment_config(segment); - if (MOVE_SEGMENT_GROUPS) + if ((MOVE_SEGMENT_GROUPS)&&(ADD_FIXED_SEGMENTS)) { for (i=0; i ymax) ymax = particle[j].yc; + } + for (j=0; j OBSTACLE_INITIAL_TIME)) evolve_segment_groups(segment, i, segment_group); } /* end of for (n=0; nINITIAL_TIME)&&(SAVE_TIME_SERIES)) + { + n_total_active = 0; + for (j=0; j INITIAL_TIME)&&(REACTION_DIFFUSION)) + { + ncollisions = update_types(particle, collisions, ncollisions); + nactive = 0; + for (j=0; j ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) + { + for (k=0; k INITIAL_TIME + WALL_TIME)&&(PRINT_ENTROPY)) { @@ -1335,8 +1418,13 @@ void animation() if (MOVE_BOUNDARY) print_segments_speeds(vxsegments, vysegments); else print_segment_group_speeds(segment_group); } + if ((i > INITIAL_TIME)&&(PLOT_PARTICLE_NUMBER)) + { + count_particle_number(particle, particle_numbers, i - INITIAL_TIME); + draw_particle_nb_plot(particle_numbers, i - INITIAL_TIME); + } - glutSwapBuffers(); + if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers(); if (MOVIE) { @@ -1365,14 +1453,19 @@ void animation() if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) { if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); - draw_particles(particle, PLOT_B, beta); + draw_particles(particle, PLOT_B, beta, collisions, ncollisions); 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); + else if (PRINT_PARTICLE_NUMBER) + { + if (REACTION_DIFFUSION) + print_particle_types_number(particle, RD_TYPES); + else 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_segment_group_speeds(segment_group); @@ -1381,6 +1474,7 @@ void animation() save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); counter++; } + else if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); /* it seems that saving too many files too fast can cause trouble with the file system */ /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ @@ -1393,20 +1487,38 @@ void animation() } } + + if (SAVE_TIME_SERIES) + { + n_total_active = 0; + for (j=0; j +#include +#include +#include +#include +#include +#include /* Sam Leffler's libtiff library. */ +#include +#include + +#define MOVIE 1 /* set to 1 to generate movie */ +#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define SAVE_MEMORY 1 /* set to 1 to save memory while saving frames */ +#define NO_EXTRA_BUFFER_SWAP 0 /* some OS require one less buffer swap when recording images */ + +#define TIME_LAPSE 1 /* set to 1 to add a time-lapse movie at the end */ + /* so far incompatible with double movie */ +#define TIME_LAPSE_FACTOR 3 /* factor of time-lapse movie */ +#define TIME_LAPSE_FIRST 1 /* set to 1 to show time-lapse version first */ + +#define SAVE_TIME_SERIES 1 /* set to 1 to save time series of particle positions */ + +/* General geometrical parameters */ + +#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 INITXMIN -2.0 +#define INITXMAX -0.07 /* x interval for initial condition */ +#define INITYMIN -1.125 +#define INITYMAX 0.75 /* y interval for initial condition */ + +#define BCXMIN -2.05 +#define BCXMAX 2.2 /* x interval for boundary condition */ +#define BCYMIN -1.125 +#define BCYMAX 1.25 /* y interval for boundary condition */ + +#define OBSXMIN -2.0 +#define OBSXMAX 2.0 /* x interval for motion of obstacle */ + +#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 181 /* 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 12 /* 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 2 /* shape of nozzle, see list in global_ljones.c */ +#define NOZZLE_SHAPE_B 4 /* 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 TYPE_PROPORTION 0.55 /* 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 */ +#define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ + +#define INTERACTION 1 /* particle interaction, see list in global_ljones.c */ +#define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */ +#define SPIN_INTER_FREQUENCY 5.0 /* angular frequency of spin-spin interaction */ +#define SPIN_INTER_FREQUENCY_B 2.0 /* angular frequency of spin-spin interaction for second particle type */ + +#define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ +#define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ +#define PDISC_DISTANCE 2.7 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_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.008 /* parameter controlling radius of particles */ +// #define MU 0.012 /* parameter controlling radius of particles */ +#define MU_B 0.012 /* parameter controlling radius of particles of second type */ +#define NPOLY 25 /* number of sides of polygon */ +#define APOLY 0.666666666 /* angle by which to turn polygon, in units of Pi/2 */ +#define MDEPTH 4 /* depth of computation of Menger gasket */ +#define MRATIO 3 /* ratio defining Menger gasket */ +#define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ +#define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ +#define FOCI 1 /* set to 1 to draw focal points of ellipse */ +#define NGRIDX 90 /* number of grid point for grid of disks */ +#define NGRIDY 100 /* number of grid point for grid of disks */ +#define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */ +#define EHRENFEST_WIDTH 0.035 /* width of tube for Ehrenfest urn configuration */ +#define TWO_CIRCLES_RADIUS_RATIO 0.8 /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */ +#define DAM_WIDTH 0.05 /* width of dam for S_DAM segment configuration */ + +#define X_SHOOTER -0.2 +#define Y_SHOOTER -0.6 +#define X_TARGET 0.4 +#define Y_TARGET 0.7 /* shooter and target positions in laser fight */ + +/* Parameters for length and speed of simulation */ + +#define NSTEPS 4000 /* number of frames of movie */ +// #define NSTEPS 6000 /* number of frames of movie */ +#define NVID 60 /* number of iterations between images displayed on screen */ +#define NSEG 250 /* number of segments of boundary */ +#define INITIAL_TIME 200 /* time after which to start saving frames */ +#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 */ + +#define PAUSE 1000 /* number of frames after which to pause */ +#define PSLEEP 1 /* sleep time during pause */ +#define SLEEP1 1 /* initial sleeping time */ +#define SLEEP2 1 /* final sleeping time */ +#define MID_FRAMES 20 /* number of still frames between parts of two-part movie */ +#define END_FRAMES 100 /* number of still frames at end of movie */ + +/* Boundary conditions, see list in global_ljones.c */ + +#define BOUNDARY_COND 0 + +/* Plot type, see list in global_ljones.c */ + +#define PLOT 8 +#define PLOT_B 0 /* plot type for second movie */ + +#define DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */ +#define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ +#define FILL_TRIANGLES 1 /* set to 1 to fill triangles between neighbours */ +#define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */ +#define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */ + +/* Color schemes */ + +#define COLOR_PALETTE 10 /* Color palette, see list in global_ljones.c */ + +#define BLACK 1 /* background */ + +#define COLOR_SCHEME 3 /* choice of color scheme, see list in global_ljones.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 ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ + +#define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ +#define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ +#define LUMMEAN 0.5 /* amplitude of luminosity variation for scheme C_LUM */ +#define LUMAMP 0.3 /* amplitude of luminosity variation for scheme C_LUM */ +#define HUEMEAN 220.0 /* mean value of hue for color scheme C_HUE */ +#define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */ + +/* particle properties */ + +#define ENERGY_HUE_MIN 330.0 /* color of original particle */ +#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.2e3 /* 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 180.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 4.5 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST_B 3.5 /* Lennard-Jones equilibrium distance for second type of particle */ +#define REPEL_RADIUS 15.0 /* radius in which repelling force acts (in units of particle radius) */ +#define DAMPING 15.0 /* damping coefficient of particles */ +#define INITIAL_DAMPING 100.0 /* damping coefficient of particles during initial phase */ +#define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ +#define PARTICLE_MASS_B 1.5 /* mass of particle of radius MU */ +#define PARTICLE_INERTIA_MOMENT 0.02 /* 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 VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ +#define SIGMA 5.0 /* noise intensity in thermostat */ +#define BETA 0.01 /* initial inverse temperature */ +#define MU_XI 0.01 /* friction constant in thermostat */ +#define KSPRING_BOUNDARY 1.0e7 /* confining harmonic potential outside simulation region */ +#define KSPRING_OBSTACLE 1.0e11 /* harmonic potential of obstacles */ +// #define NBH_DIST_FACTOR 10.0 /* radius in which to count neighbours */ +#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 /* 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 */ +#define GRAVITY_INITIAL_TIME 200 /* time at start of simulation with constant gravity */ +#define GRAVITY_RESTORE_TIME 700 /* time at end of simulation with gravity restored to initial value */ + +#define ROTATION 0 /* set to 1 to include rotation of particles */ +#define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ +#define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */ +#define KTORQUE 100.0 /* force constant in angular dynamics */ +#define KTORQUE_B 10.0 /* force constant in angular dynamics */ +#define KTORQUE_DIFF 150.0 /* force constant in angular dynamics for different particles */ +#define DRAW_SPIN 0 /* set to 1 to draw spin vectors of particles */ +#define DRAW_SPIN_B 0 /* set to 1 to draw spin vectors of particles */ +#define DRAW_CROSS 1 /* set to 1 to draw cross on particles of second type */ +#define SPIN_RANGE 7.0 /* range of spin-spin interaction */ +#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.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 2000 /* final phase in which temperature is constant */ +#define FINAL_DECREASE_PHASE 1300 /* 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 */ +#define COMPRESSION_RATIO 0.3 /* final size of container */ +#define RESTORE_CONTAINER_SIZE 1 /* set to 1 to restore container to initial size at end of simulation */ +#define RESTORE_TIME 700 /* time before end of sim at which to restore size */ + +#define MOVE_OBSTACLE 0 /* set to 1 to have a moving obstacle */ +#define CENTER_VIEW_ON_OBSTACLE 0 /* set to 1 to center display on moving obstacle */ +#define RESAMPLE_Y 0 /* set to 1 to resample y coordinate of moved particles (for shock waves) */ +#define NTRIALS 2000 /* number of trials when resampling */ +#define OBSTACLE_RADIUS 0.12 /* radius of obstacle for circle boundary conditions */ +#define FUNNEL_WIDTH 0.25 /* funnel width for funnel boundary conditions */ +#define OBSTACLE_XMIN 0.0 /* initial position of obstacle */ +#define OBSTACLE_XMAX 3.0 /* final position of obstacle */ +#define RECORD_PRESSURES 0 /* set to 1 to record pressures on obstacle */ +#define N_PRESSURES 100 /* number of intervals to record pressure */ +#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 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 */ +#define PARTIAL_THERMO_HEIGHT 0.2 /* vertical size of partial thermostat coupling */ + +#define INCREASE_KREPEL 0 /* set to 1 to increase KREPEL during simulation */ +#define KREPEL_FACTOR 1000.0 /* factor by which to change KREPEL during simulation */ + +#define PART_AT_BOTTOM 0 /* set to 1 to include "seed" particles at bottom */ +#define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */ +#define NPART_BOTTOM 100 /* number of particles at the bottom */ + +#define ADD_PARTICLES 1 /* set to 1 to add particles */ +#define ADD_TIME 1000 /* time at which to add first particle */ +#define ADD_PERIOD 3 /* time interval between adding further particles */ +#define N_ADD_PARTICLES 8 /* number of particles to add */ +#define FINAL_NOADD_PERIOD 1250 /* final period where no particles are added */ +#define SAFETY_FACTOR 2.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ + +#define TRACER_PARTICLE 0 /* set to 1 to have a tracer particle */ +#define N_TRACER_PARTICLES 3 /* number of tracer particles */ +#define TRAJECTORY_LENGTH 8000 /* length of recorded trajectory */ +#define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */ +#define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */ + +#define ROTATE_BOUNDARY 0 /* set to 1 to rotate the repelling segments */ +#define SMOOTH_ROTATION 1 /* set to 1 to update segments at each time step (rather than at each movie frame) */ +#define PERIOD_ROTATE_BOUNDARY 1000 /* period of rotating boundary */ +#define ROTATE_INITIAL_TIME 0 /* initial time without rotation */ +#define ROTATE_FINAL_TIME 100 /* final time without rotation */ +#define ROTATE_CHANGE_TIME 0.33 /* relative duration of acceleration/deceleration phases */ +#define OMEGAMAX 100.0 /* maximal rotation speed */ +#define PRINT_OMEGA 0 /* set to 1 to print angular speed */ +#define PRINT_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */ +#define PRINT_SEGMENTS_SPEEDS 1 /* set to 1 to print velocity of moving segments */ + +#define 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 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 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 0 /* 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 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) */ +#define RD_REACTION 2 /* type of reaction, see list in global_ljones.c */ +#define RD_TYPES 3 /* number of types in reaction-diffusion equation */ +#define RD_INITIAL_COND 3 /* initial condition of particles */ +#define REACION_DIST 4.0 /* maximal distance for reaction to occur */ +#define REACTION_PROB 1.0 /* probability controlling reaction term */ +// #define REACTION_PROB 0.0045 /* probability controlling reaction term */ +#define DISSOCIATION_PROB 0.005 /* probability controlling dissociation reaction */ +#define CENTER_COLLIDED_PARTICLES 0 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */ +#define COLLISION_TIME 25 /* time during which collisions are shown */ + +#define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ +#define PLOT_PARTICLE_NUMBER 0 /* set to 1 to make of plot of particle number over time */ +#define PARTICLE_NB_PLOT_FACTOR 1.0 /* expected final number of particles over initial number */ +#define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ +#define PLOT_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */ +#define PLOT_TRAJECTORIES 0 /* set to 1 to add a plot of obstacle trajectories (e.g. for rockets) */ +#define VMAX_PLOT_SPEEDS 0.6 /* vertical scale of plot of obstacle speeds */ + +#define EHRENFEST_COPY 0 /* set to 1 to add equal number of larger particles (for Ehrenfest model) */ + +#define LID_MASS 1000.0 /* mass of lid for BC_RECTANGLE_LID b.c. */ +#define LID_WIDTH 0.1 /* width of lid for BC_RECTANGLE_LID b.c. */ +#define WALL_MASS 2000.0 /* mass of wall for BC_RECTANGLE_WALL b.c. */ +#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 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.0e10 /* maximal force */ +#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 150 /* size of hashgrid in x direction */ +#define HASHY 75 /* 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 */ + +#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 8.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 12.0 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ + +#define PARTICLE_COLOR_SCHEME 1 /* choice of particle colors */ + +#define PCS_HEIGHT 0 /* color depends on final height */ +#define PCS_TREE 1 /* final position is a pine tree */ +#define PCS_CIRCLE 2 /* final position is a circle */ +#define PCS_NUMBER 3 /* final position is a number */ +#define PCS_IMAGE 4 /* final position is taken from an image file (.ppm format) */ + +#define ADD_MESSAGE 0 /* set to 1 to add message at the end */ +#define USE_RGB_COLORS 1 /* set to 1 to use RGB colors instead of hue */ +#define BACKGROUND_HUE 300.0 /* hue of particles not in image */ + +#define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO)) +#define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE)) +#define TWO_OBSTACLES ((SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)||(SEGMENT_PATTERN == S_TWO_ROCKETS)) + +double xshift = 0.0; /* x shift of shown window */ +double xspeed = 0.0; /* x speed of obstacle */ +double ylid = 0.9; /* y coordinate of lid (for BC_RECTANGLE_LID b.c.) */ +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) */ +double vysegments[2] = {SEGMENTS_VY0, SEGMENTS_VY0}; /* vy coordinate of segments (for option MOVE_BOUNDARY) */ +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" + +FILE *lj_time_series, *lj_final_position, *image_file; + +/*********************/ +/* animation part */ +/*********************/ + + +double repel_schedule(int i) +{ + static double kexponent; + static int first = 1; + double krepel; + + if (first) + { + kexponent = log(KREPEL_FACTOR)/(double)(INITIAL_TIME + NSTEPS); + first = 0; + } + krepel = KREPEL*exp(kexponent*(double)i); + printf("krepel = %.3lg\n", krepel); + return(krepel); +} + + +double temperature_schedule(int i) +{ + static double bexponent, omega, bexp2; + static int first = 1, t1, t2, t3; + double beta; + + if (first) + { + t1 = NSTEPS - MIDDLE_CONSTANT_PHASE - FINAL_DECREASE_PHASE - FINAL_CONSTANT_PHASE; + t2 = NSTEPS - FINAL_DECREASE_PHASE - FINAL_CONSTANT_PHASE; + t3 = NSTEPS - FINAL_CONSTANT_PHASE; + bexponent = log(BETA_FACTOR)/(double)(t1); + omega = N_TOSCILLATIONS*DPI/(double)(t1); + bexp2 = -log(BETA_FACTOR)/(double)(FINAL_DECREASE_PHASE); + first = 0; + } + if (i < INITIAL_TIME) beta = BETA; + else if (i < INITIAL_TIME + t1) + { + beta = BETA*exp(bexponent*(double)(i - INITIAL_TIME)); + if (!NO_OSCILLATION) beta = beta*2.0/(1.0 + cos(omega*(double)(i - INITIAL_TIME))); + } + else if (i < INITIAL_TIME + t2) beta = BETA*BETA_FACTOR; + else if (i < INITIAL_TIME + t3) + { + beta = BETA*exp(bexp2*(double)(i - INITIAL_TIME - t3)); + } + else beta = BETA; + printf("beta = %.3lg\n", beta); + return(beta); +} + +double container_size_schedule(int i) +{ + if ((i < INITIAL_TIME)||(i > INITIAL_TIME + NSTEPS - RESTORE_TIME)) return(INITXMIN); + else + return(INITXMIN + (1.0-COMPRESSION_RATIO)*(INITXMAX-INITXMIN)*(double)(i-INITIAL_TIME)/(double)(NSTEPS-RESTORE_TIME)); +} + +double obstacle_schedule_old(int i) +{ + double time; + static double t1 = 0.5, t2 = 0.75, t3 = 0.875; + + if (i < INITIAL_TIME) return(OBSTACLE_XMIN); + else + { + time = (double)(i-INITIAL_TIME)/(double)(NSTEPS); + + if (time < t1) return(OBSTACLE_XMIN + (OBSTACLE_XMAX - OBSTACLE_XMIN)*time/t1); + else if (time < t2) return(OBSTACLE_XMIN + (OBSTACLE_XMAX - OBSTACLE_XMIN)*(time - t1)/(t2 - t1)); + else if (time < t3) return(OBSTACLE_XMIN + (OBSTACLE_XMAX - OBSTACLE_XMIN)*(time - t2)/(t3 - t2)); + else return(OBSTACLE_XMAX); + } +} + +double obstacle_schedule(int i) +{ + double time, acceleration = 40.0; + double x; +// static double t1 = 0.5, t2 = 0.75, t3 = 0.875; + + if (i < INITIAL_TIME) return(OBSTACLE_XMIN); + else + { + time = (double)(i-INITIAL_TIME)/(double)(NSTEPS); + + x = OBSTACLE_XMIN + 0.5*acceleration*time*time; + xspeed = acceleration*time; +// while (x > OBSXMAX) x += OBSXMIN - OBSXMAX; + return(x); + } +} + +double obstacle_schedule_smooth(int i, int j) +{ + double time, acceleration = 50.0; + double x; +// static double t1 = 0.5, t2 = 0.75, t3 = 0.875; + + if (i < INITIAL_TIME) return(OBSTACLE_XMIN); + else + { + time = ((double)(i-INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS); + + x = OBSTACLE_XMIN + 0.5*acceleration*time*time; + xspeed = acceleration*time; +// while (x > OBSXMAX) x += OBSXMIN - OBSXMAX; + return(x); + } +} + +double gravity_schedule(int i, int j) +{ + double time, gravity, x, y; + + switch (GRAVITY_SCHEDULE){ + + case (G_INCREASE_RELEASE): + { + if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY); + else + { + time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) + + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME); + gravity = GRAVITY*(1.0 + time*(GRAVITY_FACTOR - 1.0)); + return(gravity); + } + break; + } + + case (G_INCREASE_DECREASE): + { + if ((i < INITIAL_TIME + GRAVITY_INITIAL_TIME)||(i > NSTEPS + INITIAL_TIME - GRAVITY_RESTORE_TIME)) return(GRAVITY); + else + { + time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) + + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_RESTORE_TIME - GRAVITY_INITIAL_TIME); + x = 2.0 - cos(DPI*time); + y = 0.5*((GRAVITY_FACTOR - 1.0)*x + 3.0 - GRAVITY_FACTOR); + gravity = GRAVITY*y; + return(gravity); + } + break; + } + } +} + +double rotation_angle(double phase) +{ + /* case of rotating hourglass */ +// while (phase > DPI) phase -= DPI; +// return(phase - 0.5*sin(2.0*phase)); + + /* case of centrifuge */ +// while (phase > 1.0) phase -= 1.0; +// phase *= DPI; +// angular_speed = 0.5*OMEGAMAX*(1.0 - cos(phase)); +// return(0.5*OMEGAMAX*(phase - sin(phase))); + + /* case of centrifuge remaining at constant speed for a while */ + if (phase < ROTATE_CHANGE_TIME) + { +// angular_speed = 0.5*OMEGAMAX*(1.0 - cos(phase*PI/ROTATE_CHANGE_TIME)); + return(0.5*OMEGAMAX*(phase - (ROTATE_CHANGE_TIME/PI)*sin(phase*PI/ROTATE_CHANGE_TIME))); + } + else if (phase < 1.0 - ROTATE_CHANGE_TIME) + { +// angular_speed = OMEGAMAX; + return(0.5*OMEGAMAX*(2.0*phase - ROTATE_CHANGE_TIME)); + } + else + { +// angular_speed = 0.5*OMEGAMAX*(1.0 + cos((phase - 1.0 + ROTATE_CHANGE_TIME)*PI/ROTATE_CHANGE_TIME)); + return(0.5*OMEGAMAX*(2.0 - 2.0*ROTATE_CHANGE_TIME + phase - 1.0 + (ROTATE_CHANGE_TIME/PI)*sin((1.0-phase)*PI/ROTATE_CHANGE_TIME))); + } +} + +double rotation_schedule(int i) +{ + double phase; + static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME; + + if (i < imin) + { + angular_speed = 0.0; + return(0.0); + } + else + { + if (i > imax) i = imax; + phase = (DPI/(double)PERIOD_ROTATE_BOUNDARY)*(double)(i - imin); + return(rotation_angle(phase)); + } +} + +double rotation_schedule_smooth(int i, int j) +{ + double phase, angle, phase1, angle1; + static int imin = INITIAL_TIME + ROTATE_INITIAL_TIME, imax = INITIAL_TIME + NSTEPS - ROTATE_FINAL_TIME; + + if (i < imin) + { + angular_speed = 0.0; + return(0.0); + } + else + { + if (i > imax) + { + angle = rotation_angle(1.0); + angular_speed = 0.0; + } + else + { + phase = (1.0/(double)(imax - imin))*((double)(i - imin) + (double)j/(double)NVID); + angle = rotation_angle(phase); + + phase1 = (1.0/(double)(imax - imin))*((double)(i + 1 - imin) + (double)j/(double)NVID); + angle1 = rotation_angle(phase1); + angular_speed = 25.0*(angle1 - angle); + } + + return(angle); + } +} + +int thermostat_schedule(int i) +{ + if (i < INITIAL_TIME) return(1); + else return(0); +} + +double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], + double qx[NMAXCIRCLES], double qy[NMAXCIRCLES], double qangle[NMAXCIRCLES], + double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], + double beta, int *nactive, int *nsuccess, int *nmove, int initial_phase) +{ + double a, totalenergy = 0.0, damping; + static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; + int j, move; + + if (initial_phase) damping = INITIAL_DAMPING; + else damping = DAMPING; + + #pragma omp parallel for private(j,xi,totalenergy,a,move) + for (j=0; j BCXMAX) + else if (particle[j].xc > xshift + BCXMAX) + { + particle[j].xc += BCXMIN - BCXMAX; + } + if (particle[j].yc > BCYMAX) particle[j].yc += BCYMIN - BCYMAX; + else if (particle[j].yc < BCYMIN) particle[j].yc += BCYMAX - BCYMIN; + } + else if (!NO_WRAP_BC) + { + move += wrap_particle(&particle[j], &px[j], &py[j]); + } +// if (move > 0) +// { +// compute_relative_positions(particle, hashgrid); +// update_hashgrid(particle, hashgrid, 0); /* REDUNDANT ? */ +// } + } + + return(totalenergy); +} + + +void evolve_lid(double fboundary) +{ + double force; + + force = fboundary - GRAVITY*LID_MASS; + if (ylid > BCYMAX + LID_WIDTH) force -= KSPRING_BOUNDARY*(ylid - BCYMAX - LID_WIDTH); + vylid += force*DT_PARTICLE/LID_MASS; + ylid += vylid*DT_PARTICLE; +} + +void evolve_wall(double fboundary) +{ + double force; + + force = fboundary; + if (xwall > BCYMAX - WALL_WIDTH) force -= KSPRING_BOUNDARY*(xwall - BCYMAX + WALL_WIDTH); + else if (xwall < BCYMIN + WALL_WIDTH) force += KSPRING_BOUNDARY*(BCYMIN + WALL_WIDTH - xwall); + + force -= vxwall*WALL_FRICTION; + + vxwall += fboundary*DT_PARTICLE/WALL_MASS; + + if (vxwall > WALL_VMAX) vxwall = WALL_VMAX; + else if (vxwall < -WALL_VMAX) vxwall = -WALL_VMAX; + + xwall += vxwall*DT_PARTICLE; +// printf("fboundary = %.3lg, xwall = %.3lg, vxwall = %.3lg\n", fboundary, xwall, vxwall); +} + + +void evolve_segments(t_segment segment[NMAXSEGMENTS], int time) +{ + int i, nactive = 0, group; + double fx[2] = {0.0, 0.0}, fy[2] = {0.0, 0.0}, x, y, padding = 3.0*MU, mass2 = SEGMENTS_MASS; + + if (SEGMENT_PATTERN == S_TWO_CIRCLES_EXT) mass2 = SEGMENTS_MASS*TWO_CIRCLES_RADIUS_RATIO; + + for (group=0; group<2; group++) + { + fx[group] = 0.0; + fy[group] = 0.0; + } + for (i=0; i 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); + } + if (group == 0) fy[group] -= GRAVITY*SEGMENTS_MASS; + else fy[group] -= GRAVITY*mass2; + } + if (nactive > 0) for (group=0; group<2; group++) + { + fx[group] = fx[group]/(double)nactive; + fy[group] = fy[group]/(double)nactive; + } + if (FLOOR_FORCE) + { + if (fx[0] > FMAX) fx[0] = FMAX; + else if (fx[0] < -FMAX) fx[0] = -FMAX; + if (fy[0] > FMAX) fy[0] = FMAX; + else if (fy[0] < -FMAX) fy[0] = -FMAX; + } + vxsegments[0] += fx[0]*DT_PARTICLE/SEGMENTS_MASS; + vysegments[0] += fy[0]*DT_PARTICLE/SEGMENTS_MASS; + xsegments[0] += vxsegments[0]*DT_PARTICLE; + ysegments[0] += vysegments[0]*DT_PARTICLE; + if (TWO_OBSTACLES) + { + if (FLOOR_FORCE) + { + if (fx[1] > FMAX) fx[1] = FMAX; + else if (fx[1] < -FMAX) fx[1] = -FMAX; + if (fy[1] > FMAX) fy[1] = FMAX; + else if (fy[1] < -FMAX) fy[1] = -FMAX; + } + vxsegments[1] += fx[1]*DT_PARTICLE/mass2; + vysegments[1] += fy[1]*DT_PARTICLE/mass2; + xsegments[1] += vxsegments[1]*DT_PARTICLE; + ysegments[1] += vysegments[1]*DT_PARTICLE; + } + + /* add some damping if y coordinate is small (for lunar landing) */ + if (DAMP_SEGS_AT_NEGATIVE_Y) + for (group=0; group<2; group++) + if (ysegments[group] < 0.1) + { + vysegments[group] *= exp(-DAMPING*DT_PARTICLE); + vxsegments[group] *= exp(-DAMPING*DT_PARTICLE); + } + + /* to avoid numerical instabilities */ + for (group=0; group<2; group++) + { + if (xsegments[group] + 1.0 > BCXMAX) + { + xsegments[group] = BCXMAX - 1.0; + vxsegments[group] = 0.0; + } + if ((RELEASE_ROCKET_AT_DEACTIVATION)&&((BOUNDARY_COND == BC_REFLECT_ABS)||(BOUNDARY_COND == BC_ABSORBING))) + { +// ysegments[group] = SEGMENTS_Y0; + if (time < SEGMENT_DEACTIVATION_TIME) vysegments[group] = 0.0; + else if ((ysegments[group] < SEGMENTS_Y0)&&(vysegments[group] < 0.0)) + vysegments[group] = -0.5*vysegments[group]; + } + } + + translate_segments(segment, xsegments, ysegments); +} + + +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; + padding = 0.1; + + 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 draw_particles_movie(t_particle particle[NMAXCIRCLES], int plot, double beta) +{ + int i, j, k, m, width, nnbg, nsides; + double ej, hue, huex, huey, rgb[3], rgbx[3], rgby[3], radius, x1, y1, x2, y2, angle, ca, sa, length, linkcolor, sign = 1.0, angle1, signy = 1.0, periodx, periody, x, y, lum, logratio; + char message[100]; + + if (!TRACER_PARTICLE) blank(); + glColor3f(1.0, 1.0, 1.0); + + /* show region of partial thermostat */ +// if ((PARTIAL_THERMO_COUPLING)&&(PARTIAL_THERMO_REGION == TH_INBOX)) +// { +// if (INCREASE_BETA) +// { +// logratio = log(beta/BETA)/log(0.5*BETA_FACTOR); +// if (logratio > 1.0) logratio = 1.0; +// else if (logratio < 0.0) logratio = 0.0; +// if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio; +// else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio; +// } +// else hue = 0.25*PARTICLE_HUE_MIN + 0.75*PARTICLE_HUE_MAX; +// erase_area_hsl_turbo(0.0, YMIN, 2.0*PARTIAL_THERMO_WIDTH, PARTIAL_THERMO_HEIGHT*(YMAX - YMIN), hue, 0.9, 0.15); +// +// } + + /* draw "altitude lines" */ + if (ALTITUDE_LINES) draw_altitude_lines(); + + /* draw the bonds first */ +// if ((DRAW_BONDS)||(plot == P_BONDS)) +// { +// glLineWidth(LINK_WIDTH); +// for (j=0; j= 0.0)) + for (i=-1; i<2; i++) + for (k=-1; k<2; k++) + { + x = x1 + (double)i*periodx; + y = y1 + (double)k*periody; + if (x < 1.2*particle[j].radius) + draw_one_particle(particle[j], x, y, radius, angle, nsides, width, rgb); + } + else if ((x1 >= 0.0)&&(y1 < 0.0)) + for (i=-1; i<2; i++) + for (k=-1; k<2; k++) + { + x = x1 + (double)i*periodx; + y = y1 + (double)k*periody; + if (y < 1.2*particle[j].radius) + draw_one_particle(particle[j], x, y, radius, angle, nsides, width, rgb); + } + } + } + +// /* draw spin vectors */ + if ((DRAW_SPIN)||(DRAW_SPIN_B)) + { + glLineWidth(width); + for (j=0; j 1.0) return(0.0); + if (y < 1.0/3.0) return(slope - 2.0*h - slope*y); + if (y < 2.0/3.0) return(slope - h - slope*y); + return(slope - slope*y); +} + +int tree_test(double x, double y, double ymin, double ymax) +{ + double xmax; + + if (y < ymin) return(0); + if (y > ymax) return(0); + + xmax = tree_side((y-ymin)/(ymax-ymin)); + if (vabs(x) < xmax) return(1); + else return(0); +} + +int number_test(double x, double y, double xmin, double xmax, double ymin, double ymax) +{ + int i, n; + double x1, y1, xx[4], dx, r; + + if (x < xmin) return(0); + if (x > xmax) return(0); + if (y < ymin) return(0); + if (y > ymax) return(0); + + dx = 0.25*(xmax-xmin); + for (i=0; i<4; i++) xx[i] = xmin + (double)i*dx; + + n = (int)(4.0*(x-xmin)/(xmax-xmin)); + x1 = (x - xx[n])/dx; + y1 = (y-ymin)/(ymax-ymin); + + switch(n) { + case (0): /* number 2 */ + { + if (x1 < 0.1) return(0); + if (x1 > 0.9) return(0); + if (y1 < 0.05) return(0); + if (y1 < 0.2) return(1); + if (y1 < 0.425) return(x1 < 0.25); + if (y1 < 0.575) return(1); + if (y1 < 0.8) return(x1 > 0.75); + if (y1 < 0.95) return(1); + return(0); + } + case (1): /* number 0 */ + { + r = module2(x1 - 0.5, y1 - 0.5); + if ((r < 0.45)&&(r > 0.3)) return(2); + else return(0); + } + case (2): /* number 2 */ + { + if (x1 < 0.1) return(0); + if (x1 > 0.9) return(0); + if (y1 < 0.05) return(0); + if (y1 < 0.2) return(3); + if (y1 < 0.425) return(3*(x1 < 0.25)); + if (y1 < 0.575) return(3); + if (y1 < 0.8) return(3*(x1 > 0.75)); + if (y1 < 0.95) return(3); + return(0); + } + case (3): /* number 3 */ + { + if (x1 < 0.1) return(0); + if (x1 > 0.9) return(0); + if (y1 < 0.05) return(0); + if (y1 < 0.2) return(4); + if (y1 < 0.425) return(4*(x1 > 0.75)); + if (y1 < 0.575) return(4*(x1 > 0.25)); + if (y1 < 0.8) return(4*(x1 > 0.75)); + if (y1 < 0.95) return(4); + return(0); + } + } + + return(n+1); +} + +int choose_colors(int* active_particles, t_particle *particle) +{ + int i, j, k, i1, j1, n, scan, nactive, nx, ny, maxrgb, rgbval, nmaxpixels = 100000, rgbint[3]; + double x, y, ymin = YMAX, ymax = YMIN, xmin, xmax, x1, y1, rgb[3]; + double *x_values, *y_values; + int *rgb_values; + + x_values = (double *)malloc(NMAXCIRCLES*sizeof(double)); + y_values = (double *)malloc(NMAXCIRCLES*sizeof(double)); + + scan = fscanf(lj_final_position,"%i\n", &nactive); + for (i=0; i ymax) ymax = y; + } + + if (PARTICLE_COLOR_SCHEME == PCS_IMAGE) + { + image_file = fopen("Lennard-Jones_image.ppm", "r"); + rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int)); + + scan = fscanf(image_file,"%i %i\n", &nx, &ny); + scan = fscanf(image_file,"%i\n", &maxrgb); + + printf("%i columns, %i rows, %i colors\n", nx, ny, maxrgb); + + if (nx*ny > nmaxpixels) + { + printf("Image too large, increase nmaxpixels in choose_colors()\n"); + exit(0); + } + + for (j=0; j xmax)||(y < ymin)||(y > ymax)) + { + particle[j].color_hue = BACKGROUND_HUE; + for (k=0; k<3; k++) particle[j].color_rgb[k] = rgbint[k]; + } + + else + { + x1 = (x-xmin)/(xmax-xmin); + y1 = 1.0 - (y-ymin)/(ymax-ymin); +// y1 = (y-ymin)/(ymax-ymin); + + i1 = (int)(x1*(double)nx); + j1 = (int)(y1*(double)ny); + + n = 3*(j1*nx+i1); + + if (USE_RGB_COLORS) for (k=0; k<3; k++) + { + particle[j].color_rgb[k] = rgb_values[n+k]; +// particle[j].color_rgb[k] = 256.0*(double)rgb_values[n+k]/(double)maxrgb; + } + + else particle[j].color_hue = 359.0*(1.0 - (double)rgb_values[n]/(double)maxrgb); + } + break; + } + default: particle[j].color_hue = 180.0; + } + } + + free(x_values); + free(y_values); + + if (PARTICLE_COLOR_SCHEME == PCS_IMAGE) + { + free(rgb_values); + fclose(image_file); + } + + return(nactive); +} + +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, speed_ratio, + ymin, ymax; + 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, + group, gshift, n_total_active = 0, scan, ncollisions = 0; + int *active_particles; + 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_collision *collisions; + 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 */ + 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)); + qy = (double *)malloc(NMAXCIRCLES*sizeof(double)); + px = (double *)malloc(NMAXCIRCLES*sizeof(double)); + py = (double *)malloc(NMAXCIRCLES*sizeof(double)); + qangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); + pangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); + pressure = (double *)malloc(N_PRESSURES*sizeof(double)); + + if (REACTION_DIFFUSION) + { + collisions = (t_collision *)malloc(2*NMAXCIRCLES*sizeof(t_collision)); + for (i=0; i<2*NMAXCIRCLES; i++) collisions[i].time = 0; + } + + lj_time_series = fopen("lj_time_series.dat", "r"); + lj_final_position = fopen("lj_final_position.dat", "r"); + active_particles = (int *)malloc(NMAXCIRCLES*sizeof(int)); + + + /* 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)&&(ADD_FIXED_SEGMENTS)) + { + for (i=0; i ymax) ymax = particle[j].yc; +// } +// for (j=0; j FMAX) particle[j].fx = FMAX; +// if (particle[j].fx < -FMAX) particle[j].fx = -FMAX; +// if (particle[j].fy > FMAX) particle[j].fy = FMAX; +// if (particle[j].fy < -FMAX) particle[j].fy = -FMAX; +// if (particle[j].torque > FMAX) particle[j].torque = FMAX; +// if (particle[j].torque < -FMAX) particle[j].torque = -FMAX; +// } +// } + + /* timestep of thermostat algorithm */ +// totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, beta, &nactive, &nsuccess, &nmove, i < INITIAL_TIME); + + /* evolution of lid coordinate */ +// if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(fboundary); +// if (BOUNDARY_COND == BC_RECTANGLE_WALL) +// { +// if (i < INITIAL_TIME + WALL_TIME) evolve_wall(fboundary); +// else xwall = 0.0; +// } +// if ((MOVE_BOUNDARY)&&(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)) + { + nthermo = partial_thermostat_coupling(particle, xshift + PARTIAL_THERMO_SHIFT); + printf("%i particles coupled to thermostat out of %i active\n", nthermo, nactive); + mean_energy = compute_mean_energy(particle); + } + else mean_energy = totalenergy/(double)ncircles; + +// if (CENTER_PX) center_momentum(px); +// if (CENTER_PY) center_momentum(py); +// if (CENTER_PANGLE) center_momentum(pangle); +// if (FLOOR_OMEGA) floor = floor_momentum(pangle); + +// printf("pressure left %.5lg, pressure right %.5lg\n", pleft, pright); +// for (j=0; j DPI) particle[j].angle -= DPI; + while (particle[j].angle < 0.0) particle[j].angle += DPI; + } + + + /* update tracer particle trajectory */ + if ((TRACER_PARTICLE)) + { + for (j=0; j= TRAJECTORY_LENGTH) traj_position = 0; + traj_length++; + if (traj_length >= TRAJECTORY_LENGTH) traj_length = TRAJECTORY_LENGTH - 1; +// for (j=0; j max_nb) max_nb = particle[j].neighb; +// if (particle[j].neighb < min_nb) min_nb = particle[j].neighb; +// } +// printf("Mean number of neighbours: %.3f\n", (double)total_neighbours/(double)ncircles); +// printf("Min number of neighbours: %i\n", min_nb); +// printf("Max number of neighbours: %i\n", max_nb); + + if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); + printf("nactive = %i\n", nactive); + draw_particles_movie(particle, PLOT, beta); + draw_container(xmincontainer, xmaxcontainer, obstacle, segment, wall); + + /* add a particle */ +// if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) +// { +// for (k=0; k INITIAL_TIME + WALL_TIME)&&(PRINT_ENTROPY)) + { + compute_entropy(particle, entropy); + printf("Entropy 1 = %.5lg, Entropy 2 = %.5lg\n", entropy[0], entropy[1]); + 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) + { + if (MOVE_BOUNDARY) print_segments_speeds(vxsegments, vysegments); + else print_segment_group_speeds(segment_group); + } + + glutSwapBuffers(); + + if (MOVIE) + { + if (i >= 0) + { + if (TIME_LAPSE_FIRST) + { + if ((TIME_LAPSE)&&((i - 0)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) + { + save_frame_lj(); + } + save_frame_lj_counter(NSTEPS/TIME_LAPSE_FACTOR + MID_FRAMES + i - 0); + } + else + { + save_frame_lj(); + if ((TIME_LAPSE)&&((i - 0)%TIME_LAPSE_FACTOR == 0)&&(!DOUBLE_MOVIE)) + { + save_frame_lj_counter(NSTEPS + END_FRAMES + (i - 0)/TIME_LAPSE_FACTOR); + } + } + } + else printf("Initial phase time %i of %i\n", i, 0); + + + if ((i >= 0)&&(DOUBLE_MOVIE)) + { + if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); + draw_particles(particle, PLOT_B, beta, collisions, ncollisions); + 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_segment_group_speeds(segment_group); +// print_segments_speeds(vxsegments, vysegments); + glutSwapBuffers(); + save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); + counter++; + } + + /* it seems that saving too many files too fast can cause trouble with the file system */ + /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ + if (i % PAUSE == PAUSE - 1) + { + printf("Making a short pause\n"); + sleep(PSLEEP); + s = system("mv lj*.tif tif_ljones/"); + } + } + + } + +// if (SAVE_TIME_SERIES) +// { +// n_total_active = 0; +// for (j=0; j #include #include /* Sam Leffler's libtiff library. */ +#include #include -#define MOVIE 1 /* set to 1 to generate movie */ +#define MOVIE 0 /* set to 1 to generate movie */ +#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ @@ -54,10 +56,10 @@ /* Choice of the billiard table, see global_particles.c */ -#define B_DOMAIN 30 /* choice of domain shape */ +#define B_DOMAIN 31 /* choice of domain shape */ #define CIRCLE_PATTERN 1 /* pattern of circles */ -#define POLYLINE_PATTERN 11 /* pattern of polyline */ +#define POLYLINE_PATTERN 15 /* pattern of polyline */ #define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ @@ -72,7 +74,7 @@ #define SDEPTH 1 /* Sierpinski gastket depth */ #define LAMBDA 1.5 /* parameter controlling shape of domain */ -#define MU 0.5 /* second parameter controlling shape of billiard */ +#define MU 0.005 /* second parameter controlling shape of billiard */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ @@ -87,27 +89,37 @@ /* Simulation parameters */ -#define NPART 10000 /* number of particles */ -// #define NPART 2000 /* number of particles */ +// #define NPART 1 /* number of particles */ +#define NPART 20000 /* number of particles */ +// #define NPART 10000 /* number of particles */ #define NPARTMAX 100000 /* maximal number of particles after resampling */ #define LMAX 0.01 /* minimal segment length triggering resampling */ #define DMIN 0.02 /* minimal distance to boundary for triggering resampling */ #define CYCLE 1 /* set to 1 for closed curve (start in all directions) */ -#define SHOWTRAILS 1 /* set to 1 to keep trails of the particles */ +#define SHOWTRAILS 0 /* set to 1 to keep trails of the particles */ +#define HEATMAP 1 /* set to 1 to show heat map of particles */ +#define DRAW_HEATMAP_PARTICLES 1 /* set to 1 to draw particles in heat map */ +#define HEATMAP_MAX_PART_BY_CELL 0 /* to draw only limited number of particles in cell */ +#define PLOT_HEATMAP_AVERAGE 0 /* set to 1 to plot average number of particles in heat map */ #define SHOWZOOM 0 /* set to 1 to show zoom on specific area */ #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print number of particles */ #define PRINT_LEFT_RIGHT_PARTICLE_NUMBER 1 /* set to 1 to print number of particles on left and right side */ +#define PRINT_CIRCLE_PARTICLE_NUMBER 0 /* set to 1 to print number of particles outside circular maze */ #define PRINT_COLLISION_NUMBER 0 /* set to 1 to print number of collisions */ #define TEST_ACTIVE 1 /* set to 1 to test whether particle is in billiard */ #define TEST_INITIAL_COND 0 /* set to 1 to allow only initial conditions that pass a test */ -#define NSTEPS 4500 /* number of frames of movie */ +#define NSTEPS 12300 /* number of frames of movie */ +// #define NSTEPS 6500 /* number of frames of movie */ #define TIME 1500 /* time between movie frames, for fluidity of real-time simulation */ +// #define TIME 750 /* time between movie frames, for fluidity of real-time simulation */ #define DPHI 0.00002 /* integration step */ +// #define DPHI 0.00005 /* integration step */ +// #define DPHI 0.00001 /* integration step */ #define NVID 25 /* number of iterations between images displayed on screen */ // #define NVID 100 /* number of iterations between images displayed on screen */ -#define END_FRAMES 25 /* number of still frames at the end of the movie */ +#define END_FRAMES 100 /* number of still frames at the end of the movie */ /* Decreasing TIME accelerates the animation and the movie */ /* For constant speed of movie, TIME*DPHI should be kept constant */ @@ -117,15 +129,17 @@ /* Colors and other graphical parameters */ -#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ -#define NCOLORS 128 /* number of colors */ +#define NCOLORS 1000 /* number of colors */ #define COLORSHIFT 0 /* hue of initial color */ +#define COLOR_HUEMIN 0 /* minimal color hue */ +#define COLOR_HUEMAX 150 /* maximal color hue */ #define RAINBOW_COLOR 0 /* set to 1 to use different colors for all particles */ #define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */ #define NSEG 100 /* number of segments of boundary */ -// #define LENGTH 0.01 /* length of velocity vectors */ -#define LENGTH 0.04 /* length of velocity vectors */ +#define LENGTH 0.025 /* length of velocity vectors */ +// #define LENGTH 0.04 /* length of velocity vectors */ #define BILLIARD_WIDTH 2 /* width of billiard */ #define PARTICLE_WIDTH 2 /* width of particles */ #define FRONT_WIDTH 3 /* width of wave front */ @@ -141,14 +155,14 @@ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ -#define NXMAZE 16 /* width of maze */ -#define NYMAZE 16 /* height of maze */ -// #define NXMAZE 10 /* width of maze */ -// #define NYMAZE 10 /* height of maze */ -#define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ -#define RAND_SHIFT 1 /* seed of random number generator */ +#define NXMAZE 24 /* width of maze */ +#define NYMAZE 20 /* height of maze */ +#define MAZE_MAX_NGBH 6 /* max number of neighbours of maze cell */ +#define RAND_SHIFT 11 /* seed of random number generator */ +// #define RAND_SHIFT 0 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ - +#define MAZE_RANDOM_FACTOR 0.1 /* randomization factor for S_MAZE_RANDOM */ +#define MAZE_CORNER_RADIUS 0.5 /* radius of tounded corners in maze */ #include "global_particles.c" #include "sub_maze.c" @@ -389,7 +403,7 @@ void draw_zoom(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTM // if ((active[i])&&(vabs(x1) < 1.0)&&(vabs(y1) < 1.0)&&(vabs(x2) < 1.0)&&(vabs(y2) < 1.0)) if (((active[i])&&(vabs(x1) < 1.0)&&(vabs(y1) < 1.0))||((vabs(x2) < 1.0)&&(vabs(y2) < 1.0))) { - rgb_color_scheme(color[i], rgb); + rgb_color_scheme_minmax(color[i], rgb); glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_LINE_STRIP); @@ -447,7 +461,7 @@ void draw_config_showtrails(int color[NPARTMAX], double *configs[NPARTMAX], int if (active[i]) { - rgb_color_scheme(color[i], rgb); + rgb_color_scheme_minmax(color[i], rgb); glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_LINE_STRIP); @@ -483,6 +497,105 @@ void draw_config_showtrails(int color[NPARTMAX], double *configs[NPARTMAX], int } +void draw_config_heatmap(double *configs[NPARTMAX], int active[NPARTMAX], int heatmap_number[NXMAZE*NYMAZE+1], int heatmap_total[NXMAZE*NYMAZE+1], short int heatmap_visited[NXMAZE*NYMAZE+1], int draw_particles) +/* draw a heat map of particle distribution (for mazes) */ +{ + int i, j, n, part_number; + double x, y, cosphi, sinphi, rgb[3], len, padding = 0.02; + double *xtable, *ytable; + short int *drawtable; + static int time, first = 1; + static double minprop; + + if (first) + { + time = 0; + first = 0; + minprop = 0.01; +// if (PLOT_HEATMAP_AVERAGE) minprop = 0.005; +// else minprop = 0.01; + } + time++; + + drawtable = malloc(sizeof(short int)*(NPARTMAX)); + xtable = malloc(sizeof(double)*(NPARTMAX)); + ytable = malloc(sizeof(double)*(NPARTMAX)); + + glutSwapBuffers(); + blank(); + if (PAINT_INT) paint_billiard_interior(); + + for (i=0; i configs[i][3] - padding) len = configs[i][3] - padding; + if (len < 1.0e-10) len = 1.0e-10; +// if (len < 0.0) len = 1.0e-10; + + x = configs[i][4] + len*cosphi; + y = configs[i][5] + len*sinphi; + + xtable[i] = x; + ytable[i] = y; + + /* test whether particle does not escape billiard */ + if ((TEST_ACTIVE)&&(active[i])) active[i] = xy_in_billiard(x, y); + + if (active[i]) + { + n = find_maze_cell(x, y); + + if ((n > -1)&&(n < NXMAZE*NYMAZE+1)) + { + heatmap_number[n]++; + heatmap_total[n]++; + heatmap_visited[n] = 1; + } + + if (HEATMAP_MAX_PART_BY_CELL > 0) + { + drawtable[i] = ((n == -2)||((n >= -1)&&(heatmap_number[n] <= HEATMAP_MAX_PART_BY_CELL))); + } + else drawtable[i] = (n >= -1); + } + +// printf("Particle %i is in maze cell %i\n", i, n); + } + + for (n=0; n= NCOLORS) color[i] -= NCOLORS; - } - if ((ABSORBING_CIRCLES)&&(c < 0)) active[i] = 0; - } +// if (configs[i][2]<0.0) +// { +// c = vbilliard(configs[i]); +// if (!RAINBOW_COLOR) +// { +// color[i]++; +// if (color[i] >= NCOLORS) color[i] -= NCOLORS; +// } +// if ((ABSORBING_CIRCLES)&&(c < 0)) active[i] = 0; +// } - configs[i][2] += DPHI; +// configs[i][2] += DPHI; cosphi = (configs[i][6] - configs[i][4])/configs[i][3]; sinphi = (configs[i][7] - configs[i][5])/configs[i][3]; @@ -520,11 +635,11 @@ void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPAR y2 = configs[i][5] + (configs[i][2] + LENGTH)*sinphi; /* test whether particle does not escape billiard */ - if (active[i]) active[i] = xy_in_billiard(x1, y1); + if ((TEST_ACTIVE)&&(active[i])) active[i] = xy_in_billiard(x1, y1); if (active[i]) { - rgb_color_scheme(color[i], rgb); + rgb_color_scheme_minmax(color[i], rgb); glColor3f(rgb[0], rgb[1], rgb[2]); glBegin(GL_LINE_STRIP); @@ -580,7 +695,7 @@ void draw_config(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPAR glLineWidth(3.0); } - if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3]; +// if (configs[i][2] > configs[i][3] - DPHI) configs[i][2] -= configs[i][3]; } if (DRAW_BILLIARD) draw_billiard(); @@ -608,6 +723,7 @@ void graph_movie(int time, int color[NPARTMAX], double *configs[NPARTMAX], int a for (j=0; j=0) color[i]++; if ((!RAINBOW_COLOR)&&(c>=0)) color[i]++; @@ -658,41 +775,71 @@ void print_part_number(double *configs[NPARTMAX], int active[NPARTMAX], double x } -void print_left_right_part_number(double *configs[NPARTMAX], int active[NPARTMAX], double xl, double yl, double xr, double yr, double xmin, double xmax) +void print_left_right_part_number(double *configs[NPARTMAX], int active[NPARTMAX], double xl, double yl, double xr, double yr, t_exit exits[NPARTMAX]) { char message[50]; int i, nleft = 0, nright = 0; double rgb[3], x1, y1, cosphi, sinphi; - static int maxnleft = 0, maxnright = 0; + static int first = 1; + static double xmin, xmax; + + if (first) + { + compute_maze_boundaries(POLYLINE_PATTERN, &xmin, &xmax); + first = 0; + } /* count active particles, using the fact that absorbed particles have been given dummy coordinates */ - for (i=0; i= DUMMY_ABSORBING) + for (i=0; i xmax) + if ((x1 < xmin)&&(x1 > XMIN)&&(y1 < YMAX)&&(y1 > YMIN)) exits[i].left = 1; + else if ((x1 > xmax)&&(x1 < XMAX)&&(y1 < YMAX)&&(y1 > YMIN)) { - nright++; - printf("Detected leaving particle %i at (%.2f, %2f)\n", i, x1, y1); + exits[i].right = 1; + printf("Detected leaving particle %i at (%.2f, %2f)\n\n\n", i, x1, y1); } } - if (nleft > maxnleft) maxnleft = nleft; - if (nright > maxnright) maxnright = nright; - + + for (i=0; i 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_circle_part_number(double *configs[NPARTMAX], int active[NPARTMAX], double xr, double yr) +{ + char message[50]; + int i, npart = 0; + double rgb[3], x1, y1, cosphi, sinphi; + + /* count active particles, using the fact that absorbed particles have been given dummy coordinates */ + for (i=0; i= DUMMY_ABSORBING) npart++; + + hsl_to_rgb(0.0, 0.0, 0.0, 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", maxnleft); - else sprintf(message, "%i particle", maxnleft); - write_text(xl, yl, message); - if (nright > 1) sprintf(message, "%i particles", maxnright); - else sprintf(message, "%i particle", maxnright); + if (npart > 1) sprintf(message, "%i particles", npart); + else sprintf(message, "%i particle", npart); write_text(xr, yr, message); } @@ -715,7 +862,9 @@ void animation() double time, dt, alpha, r, rgb[3], alphamax; double *configs[NPARTMAX]; int i, j, resamp = 1, s, i1, i2, c, lengthmax; - int *color, *newcolor, *active; + int *color, *newcolor, *active, *heatmap_number, *heatmap_total; + short int *heatmap_visited; + t_exit *exits; // t_circle *circles; /* experimental */ /* Since NPARTMAX can be big, it seemed wiser to use some memory allocation here */ @@ -724,6 +873,26 @@ void animation() active = malloc(sizeof(int)*(NPARTMAX)); // circles = malloc(sizeof(t_circle)*(NMAXCIRCLES)); /* experimental */ + if (HEATMAP) + { + heatmap_number = malloc(sizeof(int)*(NXMAZE*NYMAZE+1)); + heatmap_total = malloc(sizeof(int)*(NXMAZE*NYMAZE+1)); + heatmap_visited = malloc(sizeof(short int)*(NXMAZE*NYMAZE+1)); + for (i=0; i /* Sam Leffler's libtiff library. */ #define MOVIE 0 /* set to 1 to generate movie */ +#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define WINWIDTH 1280 /* window width */ #define WINHEIGHT 720 /* window height */ @@ -117,6 +118,8 @@ // #define NCOLORS 256 /* number of colors */ #define NCOLORS 48 /* number of colors */ #define COLORSHIFT 2 /* hue of initial color */ +#define COLOR_HUEMIN 0 /* minimal color hue */ +#define COLOR_HUEMAX 360 /* maximal color hue */ #define RAINBOW_COLOR 1 /* set to 1 to use different colors for all particles */ #define SINGLE_COLOR 1 /* set to 1 to make all particles a single color */ #define FLOWER_COLOR 0 /* set to 1 to adapt initial colors to flower billiard (tracks vs core) */ @@ -147,6 +150,8 @@ #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 MAZE_RANDOM_FACTOR 0.1 /* randomization factor for S_MAZE_RANDOM */ +#define MAZE_CORNER_RADIUS 0.5 /* radius of tounded corners in maze */ #define NPATHBINS 200 /* number of bins for path length histogramm */ #define PATHLMAX 1.8 /* max free path on graph */ diff --git a/rde.c b/rde.c index 552b12a..75da199 100644 --- a/rde.c +++ b/rde.c @@ -41,57 +41,64 @@ #define MOVIE 0 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ -#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ +#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ +#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ /* General geometrical parameters */ -#define WINWIDTH 1920 /* window width */ -#define WINHEIGHT 1000 /* window height */ -#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 1920 /* window width */ +// #define WINHEIGHT 1150 /* window height */ +// #define NX 960 /* number of grid points on x axis */ +// #define NY 575 /* 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 -1.0 +// #define XMAX 3.0 /* x interval */ +// #define YMIN -1.197916667 +// #define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ -// #define WINWIDTH 1280 /* window width */ -// #define WINHEIGHT 720 /* window height */ -// -// // #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 960 /* number of grid points on x axis */ -// // #define NY 540 /* number of grid points on y axis */ -// -// // #define NX 1280 /* number of grid points on x axis */ -// // #define NY 720 /* number of grid points on y axis */ -// -// #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 WINWIDTH 1280 /* window width */ +#define WINHEIGHT 720 /* window height */ + +// #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 960 /* number of grid points on x axis */ +// #define NY 540 /* number of grid points on y axis */ + +// #define NX 1280 /* number of grid points on x axis */ +// #define NY 720 /* number of grid points on y axis */ + +#define XMIN -1.0 +#define XMAX 3.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 6 /* choice of reaction term, see list in global_3d.c */ -#define NFIELDS 2 /* number of fields in reaction-diffusion equation */ -#define NLAPLACIANS 1 /* number of fields for which to compute Laplacian */ +#define RDE_EQUATION 7 /* choice of reaction term, see list in global_3d.c */ +#define NFIELDS 3 /* number of fields in reaction-diffusion equation */ +#define NLAPLACIANS 0 /* number of fields for which to compute Laplacian */ #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 ADD_FORCE_FIELD 1 /* set to 1 to add a foce field (for compressible Euler equation) */ +#define POTENTIAL 7 /* type of potential or vector potential, see list in global_3d.c */ +#define FORCE_FIELD 1 /* type of force field, see list in global_3d.c */ #define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */ +#define ADAPT_STATE_TO_BC 1 /* set to 1 to smoothly adapt initial state to obstacles */ +#define OBSTACLE_GEOMETRY 3 /* geometry of obstacles, as in B_DOMAIN */ +#define BC_STIFFNESS 50.0 /* controls region of boundary condition control */ #define JULIA_SCALE 0.5 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 197 /* choice of domain shape, see list in global_pdes.c */ +// #define B_DOMAIN 3 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 99 /* pattern of circles, see list in global_pdes.c */ @@ -99,8 +106,8 @@ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.7 /* parameter controlling the dimensions of domain */ -#define MU 0.15 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ +#define MU 0.3 /* parameter controlling the dimensions of domain */ #define NPOLY 5 /* number of sides of polygon */ #define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 7 /* depth of computation of Menger gasket */ @@ -127,10 +134,11 @@ /* Physical patameters of wave equation */ +#define DT 0.00000025 // #define DT 0.00000002 // #define DT 0.00000003 // #define DT 0.000000011 -#define DT 0.0000012 +// #define DT 0.0000012 // #define DT 0.000001 #define VISCOSITY 2.0 @@ -148,17 +156,21 @@ #define BZQ 0.0008 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */ #define B_FIELD 10.0 /* magnetic field */ +#define G_FIELD 0.01 /* gravity */ #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 K_EULER_INC 0.5 /* constant in incompressible Euler equation */ -#define SMOOTHEN_VORTICITY 1 /* set to 1 to smoothen vorticity field in Euler equation */ +#define SMOOTHEN_VORTICITY 0 /* set to 1 to smoothen vorticity field in Euler equation */ +#define SMOOTHEN_VELOCITY 1 /* set to 1 to smoothen velocity field in Euler equation */ +// #define SMOOTHEN_PERIOD 5 /* period between smoothenings */ #define SMOOTHEN_PERIOD 10 /* period between smoothenings */ -// #define SMOOTH_FACTOR 0.05 /* factor by which to smoothen */ -#define SMOOTH_FACTOR 0.03 /* factor by which to smoothen */ +#define SMOOTH_FACTOR 0.1 /* factor by which to smoothen */ +// #define SMOOTH_FACTOR 0.035 /* factor by which to smoothen */ // #define SMOOTH_FACTOR 0.015 /* factor by which to smoothen */ // #define SMOOTH_FACTOR 0.01 /* factor by which to smoothen */ -#define ADD_TRACERS 1 /* set to 1 to add tracer particles (for Euler equations) */ +#define ADD_TRACERS 0 /* set to 1 to add tracer particles (for Euler equations) */ #define N_TRACERS 1000 /* number of tracer particles */ #define T_OUT 2.0 /* outside temperature */ @@ -183,25 +195,33 @@ #define RPSLZB_INITIAL_TIME 0 /* initial time during which rpslzb remains constant */ #define RPSLZB_FINAL_TIME 500 /* final time during which rpslzb remains constant */ +#define CHANGE_FLOW_SPEED 0 /* set to 1 to change speed of laminar flow */ +#define IN_OUT_FLOW_BC 4 /* type of in-flow/out-flow boundary conditions for Euler equation */ + /* see list in global_pdes.c */ +#define IN_OUT_FLOW_MIN_AMP 0.05 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) */ +#define IN_OUT_FLOW_AMP 0.3 /* amplitude of in-flow/out-flow boundary conditions (for Euler equation) */ +#define LAMINAR_FLOW_MODULATION 0.05 /* asymmetry of laminar flow */ +#define LAMINAR_FLOW_YPERIOD 1.0 /* period of laminar flow in y direction */ + +#define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */ /* Boundary conditions, see list in global_pdes.c */ #define B_COND 1 -#define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */ - /* Parameters for length and speed of simulation */ -#define NSTEPS 2250 /* number of frames of movie */ -// #define NSTEPS 500 /* number of frames of movie */ +// #define NSTEPS 1000 /* number of frames of movie */ +#define NSTEPS 2200 /* number of frames of movie */ // #define NVID 90 /* number of iterations between images displayed on screen */ -#define NVID 120 /* number of iterations between images displayed on screen */ +#define NVID 100 /* number of iterations between images displayed on screen */ +// #define NVID 200 /* number of iterations between images displayed on screen */ // #define NVID 1100 /* number of iterations between images displayed on screen */ #define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */ #define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */ #define MAX_DT 0.024 /* maximal value of integration step */ #define NSEG 100 /* number of segments of boundary */ -#define BOUNDARY_WIDTH 5 /* width of billiard boundary */ +#define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PAUSE 100 /* number of frames after which to pause */ #define PSLEEP 2 /* sleep time during pause */ @@ -214,23 +234,23 @@ /* Visualisation */ -#define PLOT_3D 0 /* controls whether plot is 2D or 3D */ +#define PLOT_3D 1 /* 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 */ -#define DRAW_PERIODICISED 1 /* set to 1 to repeat wave periodically in x and y directions */ +#define DRAW_PERIODICISED 0 /* set to 1 to repeat wave periodically in x and y directions */ /* Plot type - color scheme */ -#define CPLOT 52 -#define CPLOT_B 51 +#define CPLOT 62 +#define CPLOT_B 61 /* Plot type - height of 3D plot */ -#define ZPLOT 52 /* z coordinate in 3D plot */ +#define ZPLOT 62 /* z coordinate in 3D plot */ // #define ZPLOT 32 /* z coordinate in 3D plot */ -#define ZPLOT_B 51 /* z coordinate in second 3D plot */ +#define ZPLOT_B 61 /* 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 */ @@ -248,6 +268,7 @@ #define PRINT_RPSLZB 0 /* set to 1 to print rpslzb parameter */ #define PRINT_PROBABILITIES 0 /* set to 1 to print probabilities (for Ehrenfest urn configuration) */ #define PRINT_NOISE 0 /* set to 1 to print noise intensity */ +#define PRINT_FLOW_SPEED 0 /* set to 1 to print speed of flow */ #define DRAW_FIELD_LINES 0 /* set to 1 to draw field lines */ #define FIELD_LINE_WIDTH 1 /* width of field lines */ @@ -266,8 +287,8 @@ /* Color schemes, see list in global_pdes.c */ -#define COLOR_PALETTE 14 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ @@ -277,12 +298,15 @@ #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 1.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_AMPLITUDE 15.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define CURL_SCALE 0.000015 /* scaling factor for curl representation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define SLOPE_SCHROD_LUM 50.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ #define MIN_SCHROD_LUM 0.2 /* minimal luminosity in color scheme Z_ARGUMENT*/ +#define VSCALE_PRESSURE 0.5 /* additional scaling factor for color scheme Z_EULER_PRESSURE */ +#define PRESSURE_SHIFT 25.0 /* shift for color scheme Z_EULER_PRESSURE */ +#define PRESSURE_LOG_SHIFT -2.5 /* shift for color scheme Z_EULER_PRESSURE */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ #define COLORDRIFT 0.0 /* how much the color hue drifts during the whole simulation */ @@ -295,6 +319,11 @@ #define LOG_SCALE 0.5 /* scaling factor for energy log representation */ #define LOG_SHIFT 1.0 #define LOG_MIN 1.0e-3 /* floor value for log vorticity plot */ +#define VSCALE_SPEED 1.5 /* additional scaling factor for color scheme Z_EULER_SPEED */ +#define VMEAN_SPEED 0.0 /* mean value around which to scale for color scheme Z_EULER_SPEED */ +#define VSCALE_DENSITY 10.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ +#define VSCALE_VORTICITY 50.0 /* additional scaling factor for color scheme Z_EULERC_VORTICITY */ +#define VORTICITY_SHIFT 0.3 /* vertical shift of vorticity */ #define NXMAZE 7 /* width of maze */ #define NYMAZE 7 /* height of maze */ @@ -303,8 +332,8 @@ #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 3.0 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 3.0 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ /* only for compatibility with wave_common.c */ @@ -316,6 +345,8 @@ #define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ #define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ #define VSCALE_ENERGY 200.0 /* additional scaling factor for color scheme P_3D_ENERGY */ +// #define VSCALE_SPEED 5.0 /* additional scaling factor for color scheme Z_EULER_SPEED */ +// #define VMEAN_SPEED 0.0 /* mean value around which to scale for color scheme Z_EULER_SPEED */ #define PHASE_FACTOR 20.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */ @@ -325,6 +356,7 @@ #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 */ +#define FLUX_WINDOW 20 /* averaging window for energy flux */ /* end of constants added only for compatibility with wave_common.c */ @@ -335,21 +367,24 @@ double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "lig double observer[3] = {8.0, 8.0, 8.0}; /* location of observer for REP_PROJ_3D representation */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ -#define Z_SCALING_FACTOR 0.08 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define Z_SCALING_FACTOR 2.4 /* overall scaling factor of z axis for REP_PROJ_3D representation */ #define XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ -#define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D 0.1 /* overall y shift for REP_PROJ_3D representation */ +#define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */ #define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */ /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ -#define VMAX 10.0 /* max value of wave amplitude */ +#define VMAX 1000.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)) +#define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)) +#define COMPUTE_PRESSURE ((ZPLOT == Z_EULER_PRESSURE)||(CPLOT == Z_EULER_PRESSURE)||(ZPLOT_B == Z_EULER_PRESSURE)||(CPLOT_B == Z_EULER_PRESSURE)) + +#define ASYM_SPEED_COLOR (VMEAN_SPEED == 0.0) #include "global_pdes.c" #include "global_3d.c" /* constants and global variables */ @@ -473,7 +508,37 @@ void compute_vector_potential(int i, int j, double *ax, double *ay) } } - +void compute_gfield(int i, int j, double *gx, double *gy) +/* initialize the exterior field, for the compressible Euler equation */ +{ + double x, y, xy[2], r, f; + + ij_to_xy(i, j, xy); + x = xy[0]; + y = xy[1]; + + switch (FORCE_FIELD) { + case (GF_VERTICAL): + { + *gx = 0.0; + *gy = -G_FIELD; + break; + } + case (GF_CIRCLE): + { + r = module2(x,y) + 1.0e-2; + f = 0.5*(1.0 - tanh(BC_STIFFNESS*(r - LAMBDA))); + *gx = G_FIELD*f*x/r; + *gy = G_FIELD*f*y/r; + break; + } + default: + { + *gx = 0.0; + *gy = 0.0; + } + } +} void initialize_potential(double potential_field[NX*NY]) /* initialize the potential field, e.g. for the Schrödinger equation */ { @@ -500,21 +565,69 @@ void initialize_vector_potential(double vpotential_field[2*NX*NY]) } } -void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY], double vector_potential_field[2*NX*NY]) +void initialize_gfield(double gfield[2*NX*NY]) +/* initialize the exterior field, e.g. for the compressible Euler equation */ +{ + int i, j; + + #pragma omp parallel for private(i,j) + for (i=0; i= SMOOTHEN_PERIOD) smooth = 0; + } + } + if (TEST_GRADIENT) { + test = 0.0; for (i=0; i<2*NX*NY; i++){ - test += nabla_omega[i]*nabla_omega[i]; - test += nabla_psi[i]*nabla_psi[i]; + test += nabla_v[i]*nabla_v[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; i 1)&&(j < NY - 1)) - { - phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*stiffness*(delta_phi[0][i*NY+j] + phi_in[1][i*NY+j]*dx*dx); + phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*stiffness*(delta_phi[0][i*NY+j] + phi_in[1][i*NY+j]*dx*dx); // phi_out[0][i*NY+j] += intstep*EULER_GRADIENT_YSHIFT; - phi_out[1][i*NY+j] = phi_in[1][i*NY+j] - intstep*K_EULER*(nabla_omega[i*NY+j]*nabla_psi[NX*NY+i*NY+j]); - phi_out[1][i*NY+j] += intstep*K_EULER*(nabla_omega[NX*NY+i*NY+j]*nabla_psi[i*NY+j]); + phi_out[1][i*NY+j] = phi_in[1][i*NY+j] - intstep*K_EULER*(nabla_omega[i*NY+j]*nabla_psi[NX*NY+i*NY+j]); + phi_out[1][i*NY+j] += intstep*K_EULER*(nabla_omega[NX*NY+i*NY+j]*nabla_psi[i*NY+j]); -// if ((i == 0)&&(j%10 == 0)) printf("j = %i, psi = %.5lg\n", j, phi_out[0][i*NY+j]); + if (COMPUTE_PRESSURE) + { + phi_out[2][i*NY+j] = phi_in[2][i*NY+j] + intstep*stiffness*(delta_pressure[i*NY+j] - delta_p[i*NY+j]); + phi_out[2][i*NY+j] *= exp(-2.0e-3); + } + break; + } + case (E_EULER_COMP): + { + rho = phi_in[0][i*NY+j]; + if (rho == 0.0) rho = 1.0e-1; + u = phi_in[1][i*NY+j]; + v = phi_in[2][i*NY+j]; + rhox = nabla_rho[i*NY+j]; + rhoy = nabla_rho[NX*NY+i*NY+j]; +// ux = nabla_u[i*NY+j]; +// uy = nabla_u[NX*NY+i*NY+j]; +// vx = nabla_v[i*NY+j]; +// vy = nabla_v[NX*NY+i*NY+j]; + + ux = rde[i*NY+j].dxu; + uy = rde[i*NY+j].dyu; + vx = rde[i*NY+j].dxv; + vy = rde[i*NY+j].dyv; + + phi_out[0][i*NY+j] = rho - intstep*(u*rhox + v*rhoy + rho*(ux + vy)); + phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + K_EULER_INC*rhox/rho); + phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + K_EULER_INC*rhoy/rho); + + if (ADD_FORCE_FIELD) + { + phi_out[1][i*NY+j] += intstep*gfield[i*NY+j]; + phi_out[2][i*NY+j] += intstep*gfield[NX*NY+i*NY+j]; } break; } @@ -652,12 +833,45 @@ void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short i } } - if (TEST_GRADIENT) { - test = 0.0; - for (i=0; i 0)) + { + switch (IN_OUT_FLOW_BC) { + case (BCE_TOPBOTTOM): + { + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, NX, 0, 10); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, NX, NY-10, NY); + break; + } + case (BCE_TOPBOTTOMLEFT): + { + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, NX, 0, 10); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, NX, NY-10, NY); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, -0.1, phi_out, xy_in, 0, 2, 0, NY); + break; + } + case (BCE_CHANNELS): + { + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, imin-5, imin+10, NY - y_maze_entry, y_maze_entry); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, imax-10, imax+5, NY- y_maze_entry, y_maze_entry); + break; + } + case (BCE_MIDDLE_STRIP): + { + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, 0, NX, NY/2 - 10, NY/2 + 10); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, 0, 2, 0, NY); + set_boundary_laminar_flow(flow_speed, LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi_out, xy_in, NX-2, NX, 0, NY); + break; + } } - printf("Delta psi + omega = %.5lg\n", test/((double)NX*NY)); + } + + if (TEST_GRADIENT) { +// test = 0.0; +// for (i=0; i 0.0)&&(v < 0.1)) @@ -783,8 +1024,10 @@ void print_parameters(t_rde rde[NX*NY], short int xy_in[NX*NY], double time, sho } else { - xbox = XMAX - 0.39; - xtext = XMAX - 0.55; + xbox = XMAX - 0.49; + xtext = XMAX - 0.65; +// xbox = XMAX - 0.39; +// xtext = XMAX - 0.55; } } else @@ -798,8 +1041,10 @@ void print_parameters(t_rde rde[NX*NY], short int xy_in[NX*NY], double time, sho } else { - xbox = XMAX - 0.39; - xtext = XMAX - 0.61; + xbox = XMAX - 0.49; + xtext = XMAX - 0.71; +// xbox = XMAX - 0.39; +// xtext = XMAX - 0.61; } } @@ -834,6 +1079,7 @@ void print_parameters(t_rde rde[NX*NY], short int xy_in[NX*NY], double time, sho else if (PRINT_VISCOSITY) sprintf(message, "Viscosity %.3f", viscosity); else if (PRINT_RPSLZB) sprintf(message, "b = %.3f", rpslzb); else if (PRINT_NOISE) sprintf(message, "noise %.3f", noise); + else if (PRINT_FLOW_SPEED) sprintf(message, "Speed %.3f", flow_speed); if (PLOT_3D) write_text(xtext, y, message); else { @@ -851,14 +1097,16 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl if (PLOT_3D) { if (ROTATE_COLOR_SCHEME) - draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); + draw_color_scheme_palette_3d(XMIN + 0.3, YMIN + 0.1, XMAX - 0.3, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value); +// draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); else draw_color_scheme_palette_3d(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value); } else { if (ROTATE_COLOR_SCHEME) - draw_color_scheme_palette_fade(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); + draw_color_scheme_palette_fade(XMIN + 0.8, YMIN + 0.1, XMAX - 0.8, YMIN + 0.1 + width, plot, -range, range, palette, fade, fade_value); +// draw_color_scheme_palette_fade(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); else draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value); } @@ -902,6 +1150,14 @@ double rpslzb_schedule(int i) } } +double flow_speed_schedule(int i) +{ + double ratio; + + ratio = (double)i/(double)NSTEPS; + return (IN_OUT_FLOW_MIN_AMP + (IN_OUT_FLOW_AMP - IN_OUT_FLOW_MIN_AMP)*ratio); +} + void viewpoint_schedule(int i) /* change position of observer */ @@ -930,7 +1186,7 @@ void animation() { double time = 0.0, scale, dx, var, jangle, cosj, sinj, sqrintstep, intstep0, viscosity_printed, fade_value, noise = NOISE_INTENSITY; - double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers; + double *phi[NFIELDS], *phi_tmp[NFIELDS], *potential_field, *vector_potential_field, *tracers, *gfield, *bc_field; short int *xy_in; int i, j, k, s, nvid, field; static int counter = 0; @@ -965,7 +1221,20 @@ void animation() initialize_vector_potential(vector_potential_field); } - if (ADD_TRACERS) tracers = (double *)malloc(2*NSTEPS*N_TRACERS*sizeof(double)); + if (ADD_FORCE_FIELD) + { + gfield = (double *)malloc(2*NX*NY*sizeof(double)); + initialize_gfield(gfield); + } + if (ADAPT_STATE_TO_BC) + { + bc_field = (double *)malloc(NX*NY*sizeof(double)); + initialize_bcfield(bc_field); + } + + +// if (ADD_TRACERS) tracers = (double *)malloc(2*NSTEPS*N_TRACERS*sizeof(double)); + if (ADD_TRACERS) tracers = (double *)malloc(4*NSTEPS*N_TRACERS*sizeof(double)); dx = (XMAX-XMIN)/((double)NX); intstep = DT/(dx*dx); @@ -987,11 +1256,18 @@ void animation() // add_coherent_state(-0.75, -0.75, 0.0, 5.0, 0.1, phi, xy_in); // init_fermion_state(-0.5, 0.5, 2.0, 0.0, 0.1, phi, xy_in); // init_boson_state(-0.5, 0.5, 2.0, 0.0, 0.1, phi, xy_in); + +// init_vortex_state(0.1, 0.4, 0.0, 0.3, -0.1, phi, xy_in); +// add_vortex_state(0.1, -0.4, 0.0, 0.3, 0.1, phi, xy_in); // init_shear_flow(1.0, 0.02, 0.15, 1, 1, phi, xy_in); -// init_laminar_flow(1.0, 0.1, 0.5, 0.0, phi, xy_in); +// init_laminar_flow(flow_speed_schedule(0), LAMINAR_FLOW_MODULATION, LAMINAR_FLOW_YPERIOD, 0.0, phi, xy_in); + init_laminar_flow(IN_OUT_FLOW_AMP, LAMINAR_FLOW_MODULATION, 0.02, 0.1, 1.0, 0.0, 0.1, phi, xy_in); - init_shear_flow(-1.0, 0.0, 0.1, 1, 1, 0.0, phi, xy_in); +// init_shear_flow(-1.0, 0.1, 0.2, 1, 1, 0.2, phi, xy_in); +// init_shear_flow(1.0, 0.02, 0.15, 1, 1, 0.0, phi, xy_in); + + if (ADAPT_STATE_TO_BC) adapt_state_to_bc(phi, bc_field, xy_in); init_cfield_rde(phi, xy_in, CPLOT, rde, 0); if (PLOT_3D) init_zfield_rde(phi, xy_in, ZPLOT, rde, 0); @@ -1043,6 +1319,8 @@ void animation() } } if (CHANGE_RPSLZB) rpslzb = rpslzb_schedule(i); + if (CHANGE_FLOW_SPEED) flow_speed = flow_speed_schedule(i); + else flow_speed = IN_OUT_FLOW_AMP; if (ROTATE_VIEW) { @@ -1065,7 +1343,9 @@ void animation() // printf("Integration step %.5lg\n", intstep); printf("Evolving wave\n"); - for (j=0; j= INITIAL_TIME)&&(DOUBLE_MOVIE)) @@ -1136,9 +1422,14 @@ void animation() if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); glutSwapBuffers(); +// if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); counter++; } + else if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); + + /* TEST */ +// if (ADAPT_STATE_TO_BC) adapt_state_to_bc(phi, bc_field, xy_in); /* it seems that saving too many files too fast can cause trouble with the file system */ /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ @@ -1161,7 +1452,8 @@ void animation() if (ADD_TRACERS) draw_tracers(phi, tracers, NSTEPS, 0, 1.0); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(rde, xy_in, time, 0, viscosity_printed, noise); - if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); + if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0); +// if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); glutSwapBuffers(); if (!FADE) for (i=0; i 360.0) hue = 360.0; + return(hue); + } + } +} + + +void set_type_color(int type, double lum, double rgb[3]) +{ + int hue; + + hue = type_hue(type); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE); + glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); +} + void init_particle_config(t_particle particles[NMAXCIRCLES]) @@ -948,6 +992,14 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) break; } + case (O_HLINE_HOLE_SPOKES): + { + radius = 2.0*MU; + + for (i=-3; i<4; i+=2) + add_obstacle(0.5*(double)i, YMIN + 0.3, radius, obstacle); + break; + } default: { printf("Function init_obstacle_config not defined for this pattern \n"); @@ -1838,6 +1890,51 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS]) break; } + case (S_HLINE_HOLE_SPOKES): + { + x = 0.15; + x1 = XMAX + 1.0; + y1 = 0.0; + width = 0.05; + + add_rectangle_to_segments(x1, y1 - width, x, y1, segment, 0); + add_rectangle_to_segments(-x, y1 - width, -x1, y1, segment, 0); + + /* closing segment */ + segment[nsegments].x1 = -x; + segment[nsegments].y1 = y1 - 0.5*width; + segment[nsegments].x2 = x; + segment[nsegments].y2 = y1 - 0.5*width; + nsegments++; + + /* spokes */ + for (i=-3; i<4; i+=2) + { + x = 0.5*(double)i; + segment[nsegments].x1 = x - 0.1; + segment[nsegments].y1 = BCYMIN - 0.1; + segment[nsegments].x2 = x; + segment[nsegments].y2 = YMIN + 0.3; + nsegments++; + segment[nsegments].x1 = x; + segment[nsegments].y1 = YMIN + 0.3; + segment[nsegments].x2 = x + 0.1; + segment[nsegments].y2 = BCYMIN - 0.1; + nsegments++; + } + + + cycle = 0; + concave = 0; /* add_rectangle_to_segments already deals with concave corners */ + + for (i=0; i 1.0 - LAMBDA + padding)&&(vabs(x) < LAMBDA - padding)) return(1); + return(0); + } default: return(1); } } @@ -2848,7 +2968,7 @@ int compute_repelling_force(int i, int j, double force[2], double *torque, t_par double distance, ca, sa, cj, sj, ca_rel, sa_rel, f[2], ff[2], torque1, ck, sk, ck_rel, sk_rel; static double distmin = 10.0*((XMAX - XMIN)/HASHX + (YMAX - YMIN)/HASHY); int interact, k; - + if (BOUNDARY_COND == BC_GENUS_TWO) distmin *= 2.0; k = particle[i].hashneighbour[j]; @@ -3147,10 +3267,7 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl } case (P_TYPE): { - if (particle.type <= 1) hue = HUE_TYPE0; - else if (particle.type == 2) hue = HUE_TYPE1; - else if (particle.type == 3) hue = HUE_TYPE2; - else hue = HUE_TYPE3; + hue = type_hue(particle.type); *radius = particle.radius; *width = BOUNDARY_WIDTH; break; @@ -3287,7 +3404,6 @@ void set_segment_group_color(int group, double lum, double rgb[3]) glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); } - void draw_altitude_lines() { int i, i1, i2; @@ -3371,10 +3487,7 @@ void draw_triangles(t_particle particle[NMAXCIRCLES], int plot) { if (plot == P_TYPE) { - if (t0 <= 1) hue = HUE_TYPE0; - else if (t0 == 2) hue = HUE_TYPE1; - else if (t0 == 3) hue = HUE_TYPE2; - else hue = HUE_TYPE3; + hue = type_hue(t0); hsl_to_rgb(hue, 0.9, 0.3, rgb); } else @@ -3423,11 +3536,143 @@ void draw_one_particle_links(t_particle particle) } } + +int draw_special_particle(t_particle particle, double xc1, double yc1, double radius, double angle, int nsides, double rgb[3], short int fill) +/* draw special particles shapes, for chemical reactions */ +{ + double x1, y1, omega; + int wsign, i; + + switch(RD_REACTION) + { + case (CHEM_AAB): + { + if (particle.type == 2) for (wsign = -1; wsign <= 1; wsign+=2) + { + x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle); + y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle); + if (fill) draw_colored_polygon(x1, y1, 0.7*radius, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 0.7*radius, nsides, angle + APOLY*PID); + return(0); + } + break; + } + case (CHEM_ABC): + { + if (particle.type == 3) for (wsign = -1; wsign <= 1; wsign+=2) + { + x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle); + y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle); + if (wsign == 1) + { + if (fill) draw_colored_polygon(x1, y1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 1.2*MU_B, nsides, angle + APOLY*PID); + } + else + { + if (fill) draw_colored_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID); + } + return(0); + } + break; + } + case (CHEM_A2BC): + { + if (particle.type == 3) + { + draw_colored_polygon(xc1, yc1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); + for (wsign = -1; wsign <= 1; wsign+=2) + { + x1 = xc1 + 1.5*radius*cos(particle.angle + 0.6*(double)wsign*PI); + y1 = yc1 + 1.5*radius*sin(particle.angle + 0.6*(double)wsign*PI); + if (fill) draw_colored_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID); + return(0); + } + } + break; + } + case (CHEM_CATALYSIS): + { + if (particle.type == 3) for (wsign = -1; wsign <= 1; wsign+=2) + { + x1 = xc1 + (double)wsign*0.7*radius*cos(particle.angle); + y1 = yc1 + (double)wsign*0.7*radius*sin(particle.angle); + if (fill) draw_colored_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 1.2*MU, nsides, angle + APOLY*PID); + return(0); + } + break; + } + case (CHEM_BAA): + { + if (particle.type == 2) for (wsign = -1; wsign <= 1; wsign+=2) + { + x1 = xc1 + (double)wsign*1.2*radius*cos(particle.angle); + y1 = yc1 + (double)wsign*1.2*radius*sin(particle.angle); + if (fill) draw_colored_polygon(x1, y1, 0.9*radius, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 0.9*radius, nsides, angle + APOLY*PID); + return(0); + } + break; + } + case (CHEM_AABAA): + { + if (particle.type == 2) for (wsign = -1; wsign <= 1; wsign+=2) + { + x1 = xc1 + (double)wsign*1.2*radius*cos(particle.angle); + y1 = yc1 + (double)wsign*1.2*radius*sin(particle.angle); + if (fill) draw_colored_polygon(x1, y1, 0.9*radius, nsides, angle + APOLY*PID, rgb); + else draw_polygon(x1, y1, 0.9*radius, nsides, angle + APOLY*PID); + return(0); + } + break; + } + case (CHEM_POLYMER): + { + if (particle.type >= 3) + { + omega = DPI/(double)(particle.type - 2); + draw_colored_polygon(xc1, yc1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); + for (i=0; i= 3) + { + omega = DPI/(double)(particle.type - 2); + draw_colored_polygon(xc1, yc1, 1.2*MU_B, nsides, angle + APOLY*PID, rgb); + for (i=0; i 0) + { + lum = (double)collisions[i].time/(double)COLLISION_TIME; + if (collisions[i].color == 0.0) for (j=0; j<3; j++) rgb[j] = lum; + else hsl_to_rgb_palette(collisions[i].color, 0.9, lum, rgb, COLOR_PALETTE); + draw_colored_polygon(collisions[i].x, collisions[i].y, 5.0*MU, NSEG, 0.0, rgb); + collisions[i].time--; + } + +} void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length) /* draw tracer particle trajectory */ @@ -3488,15 +3755,16 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int i, j, time; double x1, x2, y1, y2, rgb[3], lum; - blank(); +// blank(); glLineWidth(TRAJECTORY_WIDTH); printf("drawing trajectory\n"); for (j=0; j 0)) draw_collisions(collisions, ncollisions); /* determine particle color and size */ for (j=0; j TPYE_PROPORTION)) + if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION)) { particle[i].type = 2; particle[i].radius = MU_B; @@ -5280,10 +5597,66 @@ int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgr if (!in_segment_region(particle[i].xc, particle[i].yc)) particle[i].active = 0; - /* case of reaction-diffusion equation */ + /* case of reaction-diffusion equation/chemical reactions */ if (REACTION_DIFFUSION) for (i=0; i< ncircles; i++) { - particle[i].type = 1 + (int)(RD_TYPES*(double)rand()/(double)RAND_MAX); + switch (RD_INITIAL_COND) { + case (IC_UNIFORM): + { + particle[i].type = 1; + break; + } + case (IC_UNIFORM2): + { + particle[i].type = 2; + particle[i].radius = MU_B; + particle[i].omega = OMEGA_INITIAL*gaussian(); + break; + } + case (IC_RANDOM_UNIF): + { + particle[i].type = 1 + (int)(RD_TYPES*(double)rand()/(double)RAND_MAX); + break; + } + case (IC_RANDOM_TWO): + { + if ((double)rand()/(double)RAND_MAX < TYPE_PROPORTION) particle[i].type = 1; + else + { + particle[i].type = 2; + particle[i].radius = MU_B; + particle[i].mass_inv = 1.0/PARTICLE_MASS_B; + } + break; + } + case (IC_CIRCLE): + { + if (module2(particle[i].xc,particle[i].yc) < LAMBDA) particle[i].type = 1; + else + { + particle[i].type = 2; + particle[i].radius = MU_B; + particle[i].mass_inv = 1.0/PARTICLE_MASS_B; + } + break; + } + case (IC_CATALYSIS): + { + if ((particle[i].xc > 0.0)||((double)rand()/(double)RAND_MAX < TYPE_PROPORTION)) + { + particle[i].type = 1; + particle[i].radius = MU; + particle[i].mass_inv = 1.0/PARTICLE_MASS; + } + else + { + particle[i].type = 2; + particle[i].radius = MU_B; + particle[i].mass_inv = 1.0/PARTICLE_MASS_B; + } + break; + } + } } /* keep only active particles */ @@ -5303,6 +5676,8 @@ int add_particles(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], doub /* add several particles to the system */ { static int i = 0; + double x, y; + // add_particle(XMIN + 0.1, 0.0, 50.0, 0.0, 3.0, 0, particle); // px[ncircles - 1] = particle[ncircles - 1].vx; // py[ncircles - 1] = particle[ncircles - 1].vy; @@ -5330,9 +5705,15 @@ int add_particles(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], doub // add_particle(XMIN - 0.5*MU, 0.0, 50.0 + 5.0*(double)i, 0.0, 2.0*PARTICLE_MASS, 0, particle); +// x = INITXMIN + (INITXMAX - INITXMIN)*(double)rand()/(double)RAND_MAX; +// y = INITYMIN + (INITYMAX - INITYMIN)*(double)rand()/(double)RAND_MAX; + printf("Adding a particle\n\n"); - add_particle(BCXMIN + 0.1, 0.5*(BCYMIN + BCYMAX), 200.0, 0.0, PARTICLE_MASS, 0, particle); + x = BCXMIN + (BCXMAX - BCXMIN)*(double)rand()/(double)RAND_MAX; + y = YMAX + 0.5*(BCYMAX - YMAX)*(double)rand()/(double)RAND_MAX; + add_particle(x, y, 0.0, 0.0, PARTICLE_MASS, 0, particle); +// add_particle(BCXMIN + 0.1, 0.5*(BCYMIN + BCYMAX), 200.0, 0.0, PARTICLE_MASS, 0, particle); // i++; particle[ncircles - 1].radius = MU; @@ -5432,23 +5813,499 @@ double compute_mean_energy(t_particle particle[NMAXCIRCLES]) return(total_energy/(double)nactive); } -void update_types(t_particle particle[NMAXCIRCLES]) + +void compute_inverse_masses(double inv_masses[RD_TYPES+1]) +/* compute inverse masses of molecules in chemical reactions */ +{ + int type; + double mass; + + /* default values that apply in most cases */ + inv_masses[1] = 1.0/PARTICLE_MASS; + inv_masses[2] = 1.0/PARTICLE_MASS_B; + + switch (RD_REACTION) { + case (CHEM_AAB): + { + inv_masses[2] = 0.5/PARTICLE_MASS; + break; + } + case (CHEM_ABC): + { + inv_masses[3] = 1.0/(PARTICLE_MASS + PARTICLE_MASS_B); + break; + } + case (CHEM_A2BC): + { + inv_masses[3] = 1.0/(2.0*PARTICLE_MASS + PARTICLE_MASS_B); + break; + } + case (CHEM_CATALYSIS): + { + inv_masses[3] = 0.5/PARTICLE_MASS; + break; + } + case (CHEM_BAA): + { + inv_masses[1] = 2.0/PARTICLE_MASS_B; + break; + } + case (CHEM_AABAA): + { + inv_masses[2] = 0.5/PARTICLE_MASS; + break; + } + case (CHEM_POLYMER): + { + for (type = 3; type < RD_TYPES+1; type++) + { + mass = PARTICLE_MASS_B + (double)(type-2)*PARTICLE_MASS; + inv_masses[type] = 1.0/(1.0/PARTICLE_MASS + 1.0/mass); + } + break; + } + case (CHEM_POLYMER_DISS): + { + for (type = 3; type < RD_TYPES+1; type++) + { + mass = PARTICLE_MASS_B + (double)(type-2)*PARTICLE_MASS; + inv_masses[type] = 1.0/(1.0/PARTICLE_MASS + 1.0/mass); + } + break; + } + } +} + +int chem_merge_AAB(int i, int newtype, t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, double inv_masses[RD_TYPES+1]) +/* merging particle i with a particle of same type into a particle of type newtype */ +/* particular case of chem_merge */ +{ + int j, k, type1; + short int search = 1; + double distance, rx, ry, mr, newmass_inv; + + type1 = particle[i].type; + newmass_inv = inv_masses[newtype]; + mr = inv_masses[newtype]/inv_masses[type1]; + + for (j=0; j 10) + { + printf("Error: need to increase size of k1 table in chem_multi_merge()\n"); + exit(1); + } + + type1 = particle[i].type; + newmass_inv = inv_masses[newtype]; + mr1 = newmass_inv/inv_masses[type1]; + mr2 = newmass_inv/inv_masses[type2]; + + n1 = 0; + n2 = 0; + for (j=0; j 2)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) + ncollisions = chem_split(i, 1, particle[i].type - 1, particle, collisions, ncollisions, inv_masses); + } + + printf("%i collisions\n", ncollisions); + return(ncollisions); + } + } } @@ -5657,6 +6514,101 @@ void draw_trajectory_plot(t_group_data *group_speeds, int i) write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.1, message); } +void draw_particle_nb_plot(int *particle_numbers, int i) +/* draw plot of number of particles as a function of time */ +{ + int j, type, power; + 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, nmax, ygrad; + + if (first) + { + xmin = XMAX - 1.05; + xmax = XMAX - 0.05; + ymin = YMAX - 1.0; + ymax = YMAX - 0.05; + + xmid = 0.5*(xmin + xmax); + ymid = 0.5*(ymin + ymax); + + dx = 0.5*(xmax - xmin); + dy = 0.5*(ymax - ymin); + + plotxmin = xmin + 0.18; + plotxmax = xmax - 0.1; + plotymin = ymin + 0.07; + plotymax = ymax - 0.15; + + nmax = 0; + for (j=1; j nmax) nmax = particle_numbers[RD_TYPES+1+j]; + + nmax *= PARTICLE_NB_PLOT_FACTOR; + + power = (int)(log((double)nmax)/log(10.0)) - 1.0; + ygrad = (int)ipow(10.0, power); + + first = 0; + } + + erase_area_hsl(xmid, ymid, dx, dy, 0.0, 0.9, 0.0); + + glLineWidth(2); + + /* plot particle number */ + for (type=1; type= 1000) write_text_fixedwidth(plotxmin - 0.17, y1 - 0.015, message); + else write_text_fixedwidth(plotxmin - 0.14, y1 - 0.015, message); + } + j++; + } + + sprintf(message, "time"); + write_text_fixedwidth(plotxmax - 0.1, plotymin - 0.05, 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 */ diff --git a/sub_maze.c b/sub_maze.c index d2b441e..a5a3301 100644 --- a/sub_maze.c +++ b/sub_maze.c @@ -1,4 +1,4 @@ -/* Warning: the function init_maze does not always return a maze with a solution */ +/* The function init_maze has been improved and should return a maze with a solution */ /* The current algorithm uses a self-avoiding random walk. A better option may be */ /* to give random weights to the dual graph, and finite a maximal spanning tree */ @@ -10,6 +10,7 @@ typedef struct int neighb[MAZE_MAX_NGBH]; /* neighbour cells */ short int directions[MAZE_MAX_NGBH]; /* direction of neighbours */ short int north, east, south, west; /* closed walls */ + short int northeast, northwest, southeast, southwest; /* closed walls */ short int active; /* takes value 1 if currently active in RW path */ short int tested; /* takes value 1 if tested */ short int connected; /* takes value 1 if connected to exit */ @@ -133,7 +134,361 @@ void init_maze_graph(t_maze maze[NXMAZE*NYMAZE]) } } -int find_maze_path(t_maze maze[NXMAZE*NYMAZE], int n0, int *path, int *pathlength) + +void init_circular_maze_graph(t_maze maze[NXMAZE*NYMAZE]) +/* initialise graph of circular maze */ +/* NXMAZE should be a power of 2, and NYMAZE a multiple of NXMAZE */ +/* row number i represents the radial coordinate, and column number j the angular one */ +/* the maze is split into square blocks of size NXMAZE times NXMAZE */ +/* in each block, the first column has 1 cell, the second has 2 cells */ +/* then there are 2 columns of 4 cells, 4 columns of 8 cells, etc */ +{ + int i, j, k, n, p, q, block, nblocks, width, nextblock, prevblock; + + printf("Initializing maze\n"); + if (MAZE_MAX_NGBH < 5) + { + printf("Error: MAZE_MAX_NGBH should be at least 5 for circular maze\n"); + exit(0); + } + if (NYMAZE%NXMAZE != 0) printf("Warning: NYMAZE should be a multiple of NXMAZE\n"); + + /* set dummy variables for potentially unused cells */ + for (i=0; i 0) printf("Cell (%i, %i)\n", i, j); +// for (k = 0; k pathlength) j = 0; printf("j = %i\n", j); // while (deadend) - if (!maze[path[j]].connected) deadend = find_maze_path(maze, path[j], newpath, &newpathlength); + if (!maze[path[j]].connected) deadend = find_maze_path(maze, path[j], newpath, &newpathlength, type); if (!deadend) for (n=0; n 0) +// { +// printf("Cell (%i, %i)\t", i, j); +// if (maze[n].north) printf("North "); +// if (maze[n].northeast) printf("N-E "); +// if (maze[n].east) printf("East "); +// if (maze[n].south) printf("South "); +// if (maze[n].west) printf("West "); +// printf("\n"); +// } +// } +// sleep(5); +} + +void init_hex_maze(t_maze maze[NXMAZE*NYMAZE]) +/* init a maze with hexagonal cells */ +{ + init_maze_oftype(maze, 2); +} + void init_maze_exit(int nx, int ny, t_maze maze[NXMAZE*NYMAZE]) /* init a maze with exit at (nx, ny) */ { @@ -307,12 +772,12 @@ void init_maze_exit(int nx, int ny, t_maze maze[NXMAZE*NYMAZE]) for (i=0; i pathlength) j = 0; printf("j = %i\n", j); // while (deadend) - if (!maze[path[j]].connected) deadend = find_maze_path(maze, path[j], newpath, &newpathlength); + if (!maze[path[j]].connected) deadend = find_maze_path(maze, path[j], newpath, &newpathlength, 0); if (!deadend) for (n=0; n= 12) - { - free(image); - mem_counter = 0; - } - + if (SAVE_MEMORY) free(image); +// mem_counter++; +// if (mem_counter >= 12) +// { +// free(image); +// mem_counter = 0; +// } return 0; } @@ -162,6 +173,52 @@ void rgb_color_scheme(int i, double rgb[3]) /* color scheme */ } +void rgb_color_scheme_minmax(int i, double rgb[3]) +/* color scheme with specified color interval */ +{ + double hue, y, r; + int delta_hue; + + delta_hue = COLOR_HUEMAX - COLOR_HUEMIN; + + hue = (double)(COLOR_HUEMIN + i*delta_hue/NCOLORS); + r = 0.9; + + while (hue < COLOR_HUEMIN) hue += delta_hue; + while (hue >= COLOR_HUEMAX) hue -= delta_hue; + + hsl_to_rgb(hue, r, 0.5, rgb); + /* saturation = r, luminosity = 0.5 */ +} + +void rgb_color_scheme_density(int part_number, double rgb[3], double minprop) +/* color scheme with specified color interval */ +{ + int k; + double hue, y, r, logprop, logmin; + + if (part_number < 0) for (k=0; k<3; k++) rgb[k] = 0.25; /* visited cells */ + else + { + if (part_number<=1) hue = COLOR_HUEMAX; + else + { + logmin = log(minprop*(double)NPART); + if (logmin > 0.0) logprop = log((double)part_number)/logmin + 0.01; + else (logprop = 0.0); + + if (logprop > 1.0) logprop = 1.0; + if (logprop < 0.0) logprop = 0.0; + + hue = COLOR_HUEMAX + (COLOR_HUEMIN - COLOR_HUEMAX)*logprop; + } + r = 0.9; + hsl_to_rgb(hue, r, 0.5, rgb); + /* saturation = r, luminosity = 0.5 */ + } +} + + void rgb_color_scheme_lum(int i, double lum, double rgb[3]) /* color scheme */ { double hue, y, r; @@ -347,6 +404,28 @@ void draw_colored_circle(double x, double y, double r, int nseg, double rgb[3]) } +void draw_filled_sector(double x, double y, double rmin, double rmax, double phi1, double dphi, int nseg, double rgb[3]) +{ + int i; + double alpha, dalpha; + + dalpha = dphi/(double)nseg; + + glColor3f(rgb[0], rgb[1], rgb[2]); + for (i=0; i<=nseg; i++) + { + alpha = phi1 + (double)i*dalpha; + glBegin(GL_TRIANGLE_FAN); + glVertex2d(x + rmin*cos(alpha), y + rmin*sin(alpha)); + glVertex2d(x + rmax*cos(alpha), y + rmax*sin(alpha)); + glVertex2d(x + rmax*cos(alpha+dalpha), y + rmax*sin(alpha+dalpha)); + glVertex2d(x + rmin*cos(alpha+dalpha), y + rmin*sin(alpha+dalpha)); + glEnd(); + } + +} + + void draw_initial_condition_circle(double x, double y, double r, int color) /* draws a colored circle to mark initial condition */ { @@ -579,7 +658,7 @@ void draw_billiard() /* draws the billiard boundary */ } break; } - case D_STADIUM: + case (D_STADIUM): { glBegin(GL_LINE_LOOP); for (i=0; i<=NSEG; i++) @@ -613,7 +692,7 @@ void draw_billiard() /* draws the billiard boundary */ break; } - case D_SINAI: + case (D_SINAI): { glColor3f(0.5, 0.5, 0.5); glBegin(GL_TRIANGLE_FAN); @@ -643,7 +722,7 @@ void draw_billiard() /* draws the billiard boundary */ break; } - case D_DIAMOND: + case (D_DIAMOND): { alpha = atan(1.0 - 1.0/LAMBDA); dphi = (PID - 2.0*alpha)/(double)NSEG; @@ -1351,6 +1430,67 @@ void draw_billiard() /* draws the billiard boundary */ glVertex2d(polyline[k].x2, polyline[k].y2); glEnd(); } + + if (FOCI) + { + switch (POLYLINE_PATTERN) { + case (P_TOKARSKY): + { + glLineWidth(2); + rgb[0] = 1.0; rgb[1] = 0.0; rgb[2] = 0.0; + draw_colored_circle(-0.95, 0.0, r, NSEG, rgb); + rgb[0] = 0.0; rgb[1] = 0.8; rgb[2] = 0.2; + draw_colored_circle(0.95, 0.0, r, NSEG, rgb); + break; + } + case (P_TOKA_PRIME): + { + glLineWidth(2); + rgb[0] = 1.0; rgb[1] = 0.0; rgb[2] = 0.0; + draw_colored_circle(-polyline[84].x1, polyline[84].y1, r, NSEG, rgb); + rgb[0] = 0.0; rgb[1] = 0.8; rgb[2] = 0.2; + draw_colored_circle(polyline[84].x1, polyline[84].y1, r, NSEG, rgb); + break; + } + case (P_TOKA_NONSELF): + { + glLineWidth(2); + rgb[0] = 0.0; rgb[1] = 0.8; rgb[2] = 0.2; + draw_colored_circle(0.0, 0.0, r, NSEG, rgb); + break; + } + } + } + if (ABSORBING_CIRCLES) + { + rgb[0] = 0.7; rgb[1] = 0.7; rgb[2] = 0.7; + for (k=0; k narcs) return(0); + + len = conf[0] - (double)nside; + angle = arcs[narc].angle1 + len*(arcs[narc].dangle); + + pos[0] = arcs[narc].xc + arcs[narc].radius*cos(angle); + pos[1] = arcs[narc].yc + arcs[narc].radius*sin(angle); + + *alpha = angle + PID + conf[1]; + + return(narc + BOUNDARY_SHIFT); + } + } + + + int vpolyline_arc_xy(double config[8], double alpha, double pos[2]) + /* determine initial configuration for start at point pos = (x,y) */ + { + double c0, s0, x1, y1, a, b, c, t, dx, delta, s, xi, yi, angle, margin = 1.0e-12, tmin, rlarge = 1.0e8; + double tval[nsides + ncircles], xint[nsides + ncircles], yint[nsides + ncircles], sint[nsides + ncircles], + aint[nsides + ncircles]; + int i, nt = 0, nsegment[nsides + ncircles], narc[nsides + ncircles], ntmin, sign, type[nsides + ncircles]; + + c0 = cos(alpha); + s0 = sin(alpha); + + /* look for intersections with line segments */ + for (i=0; i margin) /* there is an intersection with the line containing segment i */ + { + t = -(a*pos[0] + b*pos[1] + c)/delta; + if (t > margin) + { + xi = pos[0] + t*c0; + yi = pos[1] + t*s0; + dx = polyline[i].x2 - polyline[i].x1; + + if (vabs(dx) > margin) s = (xi - polyline[i].x1)/dx; + else s = (yi - polyline[i].y1)/(polyline[i].y2 - polyline[i].y1); +// printf("s = %.2f\n", s); + + if ((s >= 0.0)&&(s <= 1.0)) + /* the intersection is on the segment */ + { + nsegment[nt] = i; + tval[nt] = t; + sint[nt] = s; +// xint[nt] = pos[0] + t*c0; +// yint[nt] = pos[1] + t*s0; + xint[nt] = xi; + yint[nt] = yi; + type[nt] = 0; +// printf("s = %.2f, x = %.2f, y = %.2f\n", s, xint[nt], yint[nt]); + nt++; + } + } + } + } + + /* look for intersections with arcs */ + for (i=0; i margin) for (sign=-1; sign<=1; sign+=2) + { + t = -b + (double)sign*sqrt(b*b - c); + if (t > margin) + { + xi = pos[0] + t*c0; + yi = pos[1] + t*s0; + + angle = argument(xi - arcs[i].xc, yi - arcs[i].yc); + if (angle < 0.0) angle += DPI; + + if ((angle > arcs[i].angle1)&&(angle < arcs[i].angle1 + arcs[i].dangle)) + { + narc[nt] = i; + tval[nt] = t; + sint[nt] = (angle - arcs[i].angle1)/(arcs[i].dangle); + xint[nt] = xi; + yint[nt] = yi; + aint[nt] = angle; + type[nt] = 1; +// printf("s = %.2f, x = %.2f, y = %.2f\n", sint[nt], xint[nt], yint[nt]); + nt++; + } + } + } + } + + if (ABSORBING_CIRCLES) for (i=0; i margin) /* there is an intersection with circle i */ + { + t = -b - sqrt(delta); + if (t > margin) + { + nsegment[nt] = -1-i; + + tval[nt] = t; + xint[nt] = pos[0] + t*c0; + yint[nt] = pos[1] + t*s0; + sint[nt] = argument(xint[nt] - circles[i].xc, yint[nt] - circles[i].yc); + + nt++; + } + } + } + + if (nt > 0) /* there is at least one intersection */ + { + /* find earliest intersection */ + tmin = tval[0]; + ntmin = 0; + for (i=1; i= 0)&&(type[ntmin] == 0)) /* first intersection is with a segment */ + { + config[0] = (double)nsegment[ntmin] + sint[ntmin]; + config[1] = polyline[nsegment[ntmin]].angle - alpha; + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= PI) config[1] -= DPI; + } +// else if ((narc[ntmin] >= 0)&&(type[ntmin] == 1)) /* first intersection is with a segment */ + else if ((type[ntmin] == 1)) /* first intersection is with an arc */ + { + config[0] = BOUNDARY_SHIFT + (double)narc[ntmin] + sint[ntmin]; + config[1] = PID + aint[ntmin] - alpha; + if (config[1] < 0.0) config[1] += DPI; + if (config[1] >= DPI) config[1] -= DPI; + +// printf("s = %.2f, alpha = %.2f, normal = %.2f, theta = %.2f\n", config[0], alpha*180.0/PI, aint[ntmin]*180.0/PI, config[1]*180.0/PI); + } + /* set dummy coordinates if circles are absorbing */ + else if ((ABSORBING_CIRCLES)&&(nsegment[ntmin] < 0)) + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + config[2] = 0.0; /* running time */ + config[3] = module2(xint[ntmin]-pos[0], yint[ntmin]-pos[1]); /* distance to collision */ + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = xint[ntmin]; /* position of collision */ + config[7] = yint[ntmin]; + + +// print_config(config); + + /* returned value should be positive for end of trajectory detection */ + if (nsegment[ntmin] > 0) return(nsegment[ntmin]); + else return(BOUNDARY_SHIFT-nsegment[ntmin]); + } + else /* there is no intersection - set dummy values */ + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + config[2] = 0.0; + config[3] = rlarge; + config[4] = pos[0]; /* start position */ + config[5] = pos[1]; + config[6] = rlarge*cos(alpha); + config[7] = rlarge*sin(alpha); + + return(-1); + } + } + + int vpolyline_arc(double config[8]) + /* determine initial configuration when starting from boundary */ + { + double pos[2], alpha; + int c; + + c = pos_polyline_arc(config, pos, &alpha); + + vpolyline_arc_xy(config, alpha, pos); +// c = vpolyline_xy(config, alpha, pos); + + return(c); + } + /****************************************************************************************/ /* general billiard */ /****************************************************************************************/ @@ -4891,6 +5259,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(pos_polyline(conf, pos, alpha)); break; } + case (D_POLYLINE_ARCS): + { + return(pos_polyline_arc(conf, pos, alpha)); + break; + } default: { printf("Function pos_billiard not defined for this billiard \n"); @@ -5004,6 +5377,11 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(vpolyline_xy(config, alpha, pos)); break; } + case (D_POLYLINE_ARCS): + { + return(vpolyline_arc_xy(config, alpha, pos)); + break; + } default: { printf("Function vbilliard_xy not defined for this billiard \n"); @@ -5160,6 +5538,13 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ return(vpolyline(config)); break; } + case (D_POLYLINE_ARCS): + { + c = pos_polyline_arc(config, pos, &alpha); + + return(vpolyline_arc(config)); + break; + } default: { printf("Function vbilliard not defined for this billiard \n"); @@ -5378,7 +5763,21 @@ void print_colors(int color[NPARTMAX]) /* for debugging purposes */ { /* not easy to implement for non-convex polygons */ if (POLYLINE_PATTERN == P_MAZE) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); - else if (POLYLINE_PATTERN == P_MAZE_DIAG) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); + else if ((POLYLINE_PATTERN == P_MAZE_DIAG)||(POLYLINE_PATTERN == P_MAZE_RANDOM)) + return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); + else if ((POLYLINE_PATTERN == P_MAZE_CIRCULAR)||(POLYLINE_PATTERN == P_MAZE_HEX)) + return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); + else return(1); + break; + } + case D_POLYLINE_ARCS: + { + /* not easy to implement for non-convex polygons */ + if (POLYLINE_PATTERN == P_MAZE) return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); + else if ((POLYLINE_PATTERN == P_MAZE_DIAG)||(POLYLINE_PATTERN == P_MAZE_RANDOM)) + return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); + else if (POLYLINE_PATTERN == P_MAZE_CIRCULAR) + return ((vabs(x) < 1.1*XMAX)&&(vabs(y) < 1.1*YMAX)); else return(1); break; } @@ -5744,14 +6143,403 @@ int axial_symmetry_tsegment(t_segment z1, t_segment z2, t_segment z, t_segment * return(1); } -void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) + +void init_polyline_circular_maze(t_segment* polyline, t_circle* circles, t_arc* arcs, t_maze* maze) { - int i, j, k, l, n, z, ii, jj, terni[SDEPTH], ternj[SDEPTH], quater[SDEPTH], cond; + int nblocks, block, i, j, n, p, q; + double rmin, rmax, angle, r, dr, phi, dphi, ww; + + /* build walls of maze */ + nblocks = NYMAZE/NXMAZE; + rmin = 0.15; + rmax = 1.0; + angle = DPI/(double)nblocks; + + dr = (rmax - rmin)/(double)(NXMAZE); + + nsides = 0; + ncircles = 0; + + /* add straight walls */ + for (block = 0; block < nblocks; block++) + { + dphi = angle; + + /* first circle */ + n = nmaze(0, block*NXMAZE); + r = rmin; + phi = (double)block*angle; + + if (maze[n].south) + { + polyline[nsides].x1 = r*cos(phi) + MAZE_XSHIFT; + polyline[nsides].y1 = r*sin(phi); + polyline[nsides].x2 = (r+dr)*cos(phi) + MAZE_XSHIFT; + polyline[nsides].y2 = (r+dr)*sin(phi); + polyline[nsides].angle = phi; + nsides++; + } + + /* second circle */ + r = rmin + dr; + dphi *= 0.5; + for (q=0; q<2; q++) + { + n = nmaze(1, block*NXMAZE + q); + phi = (double)(block)*angle + (double)q*dphi; + + if (maze[n].south) + { + polyline[nsides].x1 = r*cos(phi) + MAZE_XSHIFT; + polyline[nsides].y1 = r*sin(phi); + polyline[nsides].x2 = (r+dr)*cos(phi) + MAZE_XSHIFT; + polyline[nsides].y2 = (r+dr)*sin(phi); + polyline[nsides].angle = phi; + nsides++; + } + } + + /* other circles */ + ww = 2; + i = 2; + while (ww < NXMAZE) + { + dphi *= 0.5; + for (p = 0; p < ww; p++) + { + r = rmin + (double)i*dr; +// printf("Segment, i = %i, dphi = %.2lg, r = %.2lg\n", i, dphi, r); + for (q = 0; q < 2*ww; q++) + { + j = block*NXMAZE + q; + n = nmaze(i,j); + phi = (double)(block)*angle + (double)q*dphi; + + if (maze[n].south) + { + polyline[nsides].x1 = r*cos(phi) + MAZE_XSHIFT; + polyline[nsides].y1 = r*sin(phi); + polyline[nsides].x2 = (r+dr)*cos(phi) + MAZE_XSHIFT; + polyline[nsides].y2 = (r+dr)*sin(phi); + polyline[nsides].angle = phi; + nsides++; + } + } + i++; + } + ww *= 2; + } + + } + + /* add circular arcs */ + narcs = 0; + + for (block = 0; block < nblocks; block++) + { + dphi = angle; + + /* first circle */ + n = nmaze(0, block*NXMAZE); + r = rmin; + phi = (double)block*angle; + + if ((block > 0)&&(maze[n].west)) + { + arcs[narcs].xc = MAZE_XSHIFT; + arcs[narcs].yc = 0.0; + arcs[narcs].radius = r; + arcs[narcs].angle1 = phi; + arcs[narcs].dangle = dphi; + narcs++; + } + + /* second circle */ + r = rmin + dr; + dphi *= 0.5; + for (q=0; q<2; q++) + { + n = nmaze(1, block*NXMAZE + q); + phi = (double)(block)*angle + (double)q*dphi; + + if (maze[n].west) + { + arcs[narcs].xc = MAZE_XSHIFT; + arcs[narcs].yc = 0.0; + arcs[narcs].radius = r; + arcs[narcs].angle1 = phi; + arcs[narcs].dangle = dphi; + narcs++; + } + } + + /* other circles */ + ww = 2; + i = 2; + while (ww < NXMAZE) + { + dphi *= 0.5; + for (p = 0; p < ww; p++) + { + r = rmin + (double)i*dr; + printf("Circle, i = %i, dphi = %.2lg, r = %.2lg\n", i, dphi, r); + for (q = 0; q < 2*ww; q++) + { + j = block*NXMAZE + q; + n = nmaze(i,j); + phi = (double)(block)*angle + (double)q*dphi; + + if (maze[n].west) + { + arcs[narcs].xc = MAZE_XSHIFT; + arcs[narcs].yc = 0.0; + arcs[narcs].radius = r; + arcs[narcs].angle1 = phi; + arcs[narcs].dangle = dphi; + narcs++; + } + } + i++; + } + ww *= 2; + } + } + + /* outer boundary of maze */ + arcs[narcs].xc = MAZE_XSHIFT; + arcs[narcs].yc = 0.0; + arcs[narcs].radius = rmax; + arcs[narcs].angle1 = dphi; + arcs[narcs].dangle = DPI - dphi; + narcs++; + +} + +void init_polyline_hex_maze(t_segment* polyline, t_circle* circles, t_maze* maze) +{ + int i, j, n; + double r, r1, x1, y1, padding = 0.02, h; + + r = 2.0*(YMAX - YMIN)/(3.0*((double)(NXMAZE)-0.5)); + r1 = (YMAX - YMIN)/(sqrt(3.0)*(double)(NYMAZE+1)); + if (r1 < r) r = r1; + h = 0.5*sqrt(3.0)*r; + + for (i = 0; i < NXMAZE; i++) + for (j = 0; j < NYMAZE; j++) + { + n = nmaze(i, j); + x1 = YMIN + padding + 1.5*(double)i*r + MAZE_XSHIFT; + y1 = YMIN + padding + h + 2.0*(double)j*h; + if (i%2 == 0) y1 += h; + + if (((i>0)||(j!=NYMAZE/2))&&(maze[n].southwest)) + { + polyline[nsides].x1 = x1 - 0.5*r; + polyline[nsides].y1 = y1 - h; + polyline[nsides].x2 = x1 - r; + polyline[nsides].y2 = y1; + polyline[nsides].angle = DPI/3.0; + nsides++; + } + + if (maze[n].south) + { + polyline[nsides].x1 = x1 - 0.5*r; + polyline[nsides].y1 = y1 - h; + polyline[nsides].x2 = x1 + 0.5*r; + polyline[nsides].y2 = y1 - h; + polyline[nsides].angle = 0.0; + nsides++; + } + + if (maze[n].northwest) + { + polyline[nsides].x1 = x1 - r; + polyline[nsides].y1 = y1; + polyline[nsides].x2 = x1 - 0.5*r; + polyline[nsides].y2 = y1 + h; + polyline[nsides].angle = PI/3.0; + nsides++; + } + + if (((j==0)||(i==NXMAZE-1))&&(maze[n].southeast)) + { + polyline[nsides].x1 = x1 + 0.5*r; + polyline[nsides].y1 = y1 - h; + polyline[nsides].x2 = x1 + r; + polyline[nsides].y2 = y1; + polyline[nsides].angle = PI/3.0; + nsides++; + } + + if ((j==NYMAZE-1)&&(maze[n].north)) + { + polyline[nsides].x1 = x1 - 0.5*r; + polyline[nsides].y1 = y1 + h; + polyline[nsides].x2 = x1 + 0.5*r; + polyline[nsides].y2 = y1 + h; + polyline[nsides].angle = 0.0; + nsides++; + } + + if (((j==NYMAZE-1)||((i==NXMAZE-1)&&(j!=NYMAZE/2-1)))&&(maze[n].northeast)) + { + polyline[nsides].x1 = x1 + r; + polyline[nsides].y1 = y1; + polyline[nsides].x2 = x1 + 0.5*r; + polyline[nsides].y2 = y1 + h; + polyline[nsides].angle = DPI/3.0; + nsides++; + } + } + /* left side of maze */ + x1 = YMIN + padding + MAZE_XSHIFT - 0.5*r; + y1 = YMIN + padding + h; + polyline[nsides].x1 = x1 + 1.5*r; + polyline[nsides].y1 = YMIN - MAZE_VERTICAL_EXTENSION; + polyline[nsides].x2 = x1 + 1.5*r; + polyline[nsides].y2 = y1 - h; + polyline[nsides].angle = PID; + nsides++; + + y1 = YMIN + padding + 3.0*h + 2.0*(double)(NYMAZE-1)*h; + polyline[nsides].x1 = x1; + polyline[nsides].y1 = y1; + polyline[nsides].x2 = x1; + polyline[nsides].y2 = YMAX + MAZE_VERTICAL_EXTENSION; + polyline[nsides].angle = PID; + nsides++; + + /* right side of maze */ + x1 = YMIN + padding + 1.5*(double)(NXMAZE-1)*r + MAZE_XSHIFT + 0.5*r; + y1 = YMIN + padding; + polyline[nsides].x1 = x1; + polyline[nsides].y1 = YMIN - MAZE_VERTICAL_EXTENSION; + polyline[nsides].x2 = x1; + polyline[nsides].y2 = y1; + polyline[nsides].angle = PID; + nsides++; + + y1 = YMIN + padding + 2.0*h + 2.0*(double)(NYMAZE-1)*h; + polyline[nsides].x1 = x1 - 1.5*r; + polyline[nsides].y1 = y1 + h; + polyline[nsides].x2 = x1 - 1.5*r; + polyline[nsides].y2 = YMAX + MAZE_VERTICAL_EXTENSION; + polyline[nsides].angle = PID; + nsides++; + + /* add circular arcs in corners */ + if (B_DOMAIN == D_POLYLINE_ARCS) + { + narcs = 0; + + for (i=0; i 0)||(j!=NYMAZE/2))&&(maze[n].north)&&(maze[n].northwest)) + { + arcs[narcs].xc = x1; + arcs[narcs].yc = y1; + arcs[narcs].radius = h; + arcs[narcs].dangle = PI/3.0; + arcs[narcs].angle1 = PID; + narcs++; + } + if (((i > 0)||(j!=NYMAZE/2))&&(maze[n].northwest)&&(maze[n].southwest)) + { + arcs[narcs].xc = x1; + arcs[narcs].yc = y1; + arcs[narcs].radius = h; + arcs[narcs].dangle = PI/3.0; + arcs[narcs].angle1 = 5.0*PI/6.0; + narcs++; + } + if (((i > 0)||(j!=NYMAZE/2))&&(maze[n].southwest)&&(maze[n].south)) + { + arcs[narcs].xc = x1; + arcs[narcs].yc = y1; + arcs[narcs].radius = h; + arcs[narcs].dangle = PI/3.0; + arcs[narcs].angle1 = 7.0*PI/6.0; + narcs++; + } + if (((i > 0)||(j!=NYMAZE/2))&&(maze[n].south)&&(maze[n].southeast)) + { + arcs[narcs].xc = x1; + arcs[narcs].yc = y1; + arcs[narcs].radius = h; + arcs[narcs].dangle = PI/3.0; + arcs[narcs].angle1 = 3.0*PID; + narcs++; + } + if (((i < NXMAZE-1)||(j!=NYMAZE/2-1))&&(maze[n].southeast)&&(maze[n].northeast)) + { + arcs[narcs].xc = x1; + arcs[narcs].yc = y1; + arcs[narcs].radius = h; + arcs[narcs].dangle = PI/3.0; + arcs[narcs].angle1 = 11.0*PI/6.0; + narcs++; + } + } + } +} + +// (()||())&& + +void compute_maze_boundaries(int type, double *xmin, double *xmax) +{ + double r, r1, h, padding = 0.02; + + switch (type) { + case (P_MAZE_HEX): + { + r = 2.0*(YMAX - YMIN)/(3.0*((double)(NXMAZE)-0.5)); + r1 = 2.0*(YMAX - YMIN)/(sqrt(3.0)*(double)(NYMAZE)); + if (r1 < r) r = r1; + + *xmin = YMIN + padding + MAZE_XSHIFT - 2.0*r; + *xmax = YMIN + padding + 1.5*(double)(NXMAZE-1)*r + MAZE_XSHIFT + 2.0*r; + + printf("xmin = %.3lg, xmax = %.3lg\n", *xmin, *xmax); + break; + } + default: + { + *xmin = YMIN + MAZE_XSHIFT; + *xmax = YMAX + MAZE_XSHIFT; + break; + } + } + +} + +void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES], t_arc arcs[NMAXCIRCLES]) +{ + int i, j, k, l, n, z, ii, jj, terni[SDEPTH], ternj[SDEPTH], quater[SDEPTH], cond, p, q, block, nblocks, ww; short int vkoch[NMAXCIRCLES], turnright; - 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; + double ratio, omega, angle, sw, length, dist, x, y, ta, tb, a, b, ra, rb, r, + x1, y1, x2, y2, dx, dy, padding = 0.02, width = 0.01, xright, yright, xtop, ytop, rand_factor, rmin, rmax, dr, phi, dphi, rscat; + double *maze_coords; t_maze* maze; + narcs = 0; /* default value */ + switch (POLYLINE_PATTERN) { case (P_RECTANGLE): { @@ -6203,10 +6991,8 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) 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++; } - nsides++; if (maze[n].south) { @@ -6215,11 +7001,8 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) 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++; } - - nsides++; } /* top side of maze */ @@ -6262,6 +7045,63 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) polyline[nsides].y2 = YMAX + 1000.0; polyline[nsides].angle = PID; nsides++; + + /* add circular arcs in corners */ + if (B_DOMAIN == D_POLYLINE_ARCS) + { + narcs = 0; + + for (i=0; i0)||(j!=NYMAZE/2))&&(maze[n].south)&&(maze[n].west)) + { + arcs[narcs].xc = x1 + ra*dx; + arcs[narcs].yc = y1 + ra*dy; + arcs[narcs].radius = ra*dx; + arcs[narcs].angle1 = PI; + arcs[narcs].dangle = PID; + narcs++; + } + + if (((i0)||(j!=NYMAZE/2))&&(maze[n].north)&&(maze[n].west)) + { + arcs[narcs].xc = x1 + ra*dx; + arcs[narcs].yc = y1 + rb*dy; + arcs[narcs].radius = ra*dx; + arcs[narcs].angle1 = PID; + arcs[narcs].dangle = PID; + narcs++; + } + } + } free(maze); break; @@ -6278,6 +7118,9 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) nsides = 0; ncircles = 0; + + /* deactivate sides for openings */ + maze[nmaze(0,NYMAZE/2)].west = 0; for (i=0; i0)||(j!=NYMAZE/2))&&(maze[n].south)&&(maze[n].west)) { polyline[nsides].x1 = x1 + 0.25*dx; polyline[nsides].y1 = y1; @@ -6327,8 +7170,8 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) nsides++; } - if ((maze[n].south)&&(maze[n].east)) - { +// if ((maze[n].south)&&(maze[n].east)) + if (((i0)||(j!=NYMAZE/2))&&(maze[n].north)&&(maze[n].west)) { polyline[nsides].x1 = x1; polyline[nsides].y1 = y1 + 0.75*dy; @@ -6402,6 +7246,249 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES]) free(maze); break; } + case (P_MAZE_RANDOM): + { + maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze)); + maze_coords = (double *)malloc(2*NXMAZE*NYMAZE*sizeof(double)); + + init_maze_exit(0, NYMAZE/2, maze); + + /* build walls of maze */ + dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE); + dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE); + + /* compute maze coordinates with random shift */ + for (i=0; i0)||(j!=NYMAZE/2))&&(maze[n].west)) + { + polyline[nsides].x1 = x1; + polyline[nsides].y1 = y1; + polyline[nsides].x2 = xtop; + polyline[nsides].y2 = ytop; + polyline[nsides].angle = PID; + nsides++; + } + + if (maze[n].south) + { + polyline[nsides].x1 = x1; + polyline[nsides].y1 = y1; + polyline[nsides].x2 = xright; + polyline[nsides].y2 = yright; + polyline[nsides].angle = 0.0; + nsides++; + } + } +// sleep(3); + + /* 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 - 1000.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 + 1000.0; + polyline[nsides].angle = PID; + nsides++; + + /* left side of maze */ + x1 = YMIN + padding + MAZE_XSHIFT; + polyline[nsides].x1 = x1; + polyline[nsides].y1 = YMIN - 1000.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 + 1000.0; + polyline[nsides].angle = PID; + nsides++; + +// ncircles = nsides; +// for (i=0; i NXMAZE-1) return(-1); + if (j < 0) return(-1); + if (j > NYMAZE-1) return(-1); + + x1 = x - YMIN - padding - MAZE_XSHIFT - (double)i*dx; + y1 = y - YMIN - padding - (double)j*dy; + + /* avoid finding wrong cell for particles close to wall */ + if (x1 < tolerance*dx) return(-10); + if (x1 > dx - tolerance*dx) return(-10); + if (y1 < tolerance*dy) return(-10); + if (y1 > dy - tolerance*dy) return(-10); + + n = nmaze(i, j); + + if (n < NXMAZE*NYMAZE) return(n); + else return(-1); + } + case (1): /* circular maze */ + { + r = module2(x - MAZE_XSHIFT,y); + phi = argument(x - MAZE_XSHIFT,y); + if (phi < 0.0) phi += DPI; + if (phi >= DPI) phi -= DPI; + + if (r < rmin*(1.0 - tolerance)) return(NXMAZE*NYMAZE); + if (r < rmin*(1.0 + tolerance)) return(-10); + if (r > rmax) return(-2); + + i = (int)((r - rmin)/dr); + if (i > NXMAZE-1) return(-1); + + /* avoid finding wrong cell for particles close to wall */ + if (r - rmin - (double)i*dr < tolerance*dr) return(-10); + if (r - rmin - (double)(i+1)*dr > -tolerance*dr) return(-10); + + block = (int)(phi/angle); + phi1 = phi - (double)block*angle; + + /* avoid finding wrong cell for particles close to wall */ + if (phi1 < tolerance*angle) + return(-10); + if (phi1 - angle > -tolerance*angle) + return(-10); + + if (i == 0) j = block*NXMAZE; + else + { + w = 1 + (int)(log((double)i)/log(2.0)); +// dphi = angle/ipow(2.0, w); + dphi = dphi_table[w]; + q = (int)(phi1/dphi); + j = block*NXMAZE + q; + + /* avoid finding wrong cell for particles close to wall */ + if (phi1 - (double)q*dphi < tolerance*dphi) + return(-1); + if (phi1 - (double)(q+1)*dphi > - tolerance*dphi) + return(-1); + } + + n = nmaze(i, j); + + if (n < NXMAZE*NYMAZE) return(n); + else return(-1); + } + case (2): /* honeycomb maze */ + { + /* transformation to coordinate system along hexagon centers */ + u = (x-x0)/(1.5*rr); + v = 0.5*(u + (y-y0)/h); + + i = (int)u; + j = (int)v; + + j -= i/2; + + if (i < 0) return(-1); + if (i > NXMAZE-1) return(-1); + if (j < 0) return(-1); + if (j > NYMAZE-1) return(-1); + + for (k=-1; k<2; k++) if ((i+k >= 0)&&(i+k < NXMAZE)) + for (l=-1; l<2; l++) if ((j+l >= 0)&&(j+l < NYMAZE)) + { + x1 = x0 + 1.5*(double)(i+k)*rr; + y1 = y0 + 2.0*(double)(j+l)*h; + if ((i+k)%2 == 0) y1 += h; + if (in_polygon(x - x1, y - y1, rrtol, 6, 0.0)) + { + n = nmaze(i+k, j+l); + + if (i+k < 0) return(-1); + if (i+k > NXMAZE-1) return(-1); + if (j+l < 0) return(-1); + if (j+l > NYMAZE-1) return(-1); + + if (n < NXMAZE*NYMAZE) return(n); + else return(-1); + } + } + return(-10); + } + } +} + +void draw_maze_cell(int n, int part_number, double minprop) +/* give specified color to maze cell */ +{ + int i, j, block, q; + double x, y, padding = 0.02, rgb[3], w; + static double dx, dy, rmin, rmax, angle, r, r1, dr, phi1, dphi, h, x0, y0; + static int first = 1, nblocks; + + if (first) switch (maze_type(POLYLINE_PATTERN)) { + case (0): /* maze with square or rectangular cells */ + { + dx = (YMAX - YMIN - 2.0*padding)/(double)(NXMAZE); + dy = (YMAX - YMIN - 2.0*padding)/(double)(NYMAZE); + break; + } + case (1): /* circular maze */ + { + nblocks = NYMAZE/NXMAZE; + rmin = 0.15; + rmax = 1.0; + angle = DPI/(double)nblocks; + + dr = (rmax - rmin)/(double)(NXMAZE); + + break; + } + case (2): /* honeycomb maze */ + { + r = 2.0*(YMAX - YMIN)/(3.0*((double)(NXMAZE)-0.5)); + r1 = (YMAX - YMIN)/(sqrt(3.0)*(double)(NYMAZE+1)); + if (r1 < r) r = r1; + h = 0.5*sqrt(3.0)*r; + x0 = YMIN + padding + MAZE_XSHIFT; + y0 = YMIN + padding + h; + break; + } + first = 0; + } + + if (part_number != 0) switch (maze_type(POLYLINE_PATTERN)) { + case (0): /* maze with square or rectangular cells */ + { + i = n%NXMAZE; + j = n/NXMAZE; + + x = YMIN + padding + (double)i*dx + MAZE_XSHIFT; + y = YMIN + padding + (double)j*dy; + + rgb_color_scheme_density(part_number, rgb, minprop); + erase_rectangle(x, y, x + dx, y + dy, rgb); + break; + } + case (1): /* circular maze */ + { + if (n == NXMAZE*NYMAZE) /* inner circle */ + { + rgb_color_scheme_density(part_number, rgb, minprop); + draw_filled_sector(MAZE_XSHIFT, 0.0, 0.0, rmin, 0.0, DPI, NSEG, rgb); + break; + } + + i = n%NXMAZE; + j = n/NXMAZE; + block = j/NXMAZE; + + r = rmin + (double)i*dr; + r1 = r + dr; + if (i == 0) /* first circle */ + { + dphi = angle; + phi1 = (double)block*dphi; + } + else /* other circles */ + { + w = 1 + (int)(log((double)i)/log(2.0)); + dphi = angle*ipow(0.5, w); + q = j - block*NXMAZE; + phi1 = (double)block*angle + (double)q*dphi; + } + rgb_color_scheme_density(part_number, rgb, minprop); + draw_filled_sector(MAZE_XSHIFT, 0.0, r, r1, phi1, dphi, NSEG, rgb); + break; + } + case (2): /* honeycomb maze */ + { + i = n%NXMAZE; + j = n/NXMAZE; + + x = x0 + 1.5*(double)i*r; + y = y0 + 2.0*(double)j*h; + if (i%2 == 0) y += h; + + rgb_color_scheme_density(part_number, rgb, minprop); + draw_colored_circle(x, y, r, 6, rgb); + break; + } + } + +} + + +/* can probably be removed */ +void print_left_right_part_number_old(double *configs[NPARTMAX], int active[NPARTMAX], double xl, double yl, double xr, double yr, double xmin, double xmax) +{ + char message[50]; + int i, nleft = 0, nright = 0; + double rgb[3], x1, y1, cosphi, sinphi; + static int maxnleft = 0, maxnright = 0; + + /* count active particles, using the fact that absorbed particles have been given dummy coordinates */ + for (i=0; i= DUMMY_ABSORBING)) + { + cosphi = (configs[i][6] - configs[i][4])/configs[i][3]; + sinphi = (configs[i][7] - configs[i][5])/configs[i][3]; + x1 = configs[i][4] + configs[i][2]*cosphi; + y1 = configs[i][5] + configs[i][2]*cosphi; + if (x1 < xmin) nleft++; +// else if ((x1 > xmax)&&(y1 <= YMAX)&&(y1 >= YMIN)) + else if ((x1 > xmax)) + { + nright++; + printf("Detected leaving particle %i at (%.2f, %2f)\n", i, x1, y1); + } + } + if (nleft > maxnleft) maxnleft = nleft; + if (nright > maxnright) maxnright = nright; + + hsl_to_rgb(0.0, 0.0, 0.0, rgb); + + erase_area(xl, yl - 0.03, 0.45, 0.12, rgb); + erase_area(xr, yr - 0.03, 0.35, 0.12, rgb); + + glColor3f(1.0, 1.0, 1.0); + if (maxnleft > 1) sprintf(message, "%i particles", maxnleft); + else sprintf(message, "%i particle", maxnleft); + write_text(xl, yl, message); + if (maxnright > 1) sprintf(message, "%i particles", maxnright); + else sprintf(message, "%i particle", maxnright); + write_text(xr, yr, message); +} diff --git a/sub_rde.c b/sub_rde.c index bdb6b86..16dc790 100644 --- a/sub_rde.c +++ b/sub_rde.c @@ -4,6 +4,7 @@ double intstep; /* integration step */ double intstep1; /* integration step used in absorbing boundary conditions */ double viscosity; /* viscosity (parameter in front of Laplacian) */ double rpslzb; /* second parameter in Rock-Paper-Scissors-Lizard-Spock equation */ +double flow_speed; /* flow speed for laminar boundary conditions in Euler equation */ double gaussian() /* returns standard normal random variable, using Box-Mueller algorithm */ @@ -313,70 +314,146 @@ void antisymmetrize_wave_function(double *phi[NFIELDS], short int xy_in[NX*NY]) } } -void init_vortex_state(double x, double y, double scale, double *phi[NFIELDS], short int xy_in[NX*NY]) -/* initialise field with vortex at position (x,y) with variance scale */ -/* phi[0] is stream function, phi[1] is vorticity */ +void init_vortex_state(double amp, double x, double y, double scale, double density_mod, double *phi[NFIELDS], short int xy_in[NX*NY]) +/* initialise field with vortex at position (x,y) with amplitude and variance scale */ +/* for incompressible Euler, phi[0] is stream function, phi[1] is vorticity */ +/* for compressible Euler, phi[0] is the density, phi[1] and phi[2] are velocity components */ { int i, j; double xy[2], dist2, module, phase, scale2, sign; // scale2 = scale*scale; - for (i=0; i 0.0) sign = 1.0; - else sign = -1.0; + if (scale > 0.0) sign = 1.0; + else sign = -1.0; - phi[1][i*NY+j] += sign*module; - phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */ - } - else - { - phi[0][i*NY+j] = 0.0; - phi[1][i*NY+j] = 0.0; - } + phi[1][i*NY+j] += sign*module; + phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */ + } + else + { + phi[0][i*NY+j] = 0.0; + phi[1][i*NY+j] = 0.0; + } + } + break; } + case (E_EULER_COMP): + { + for (i=0; i 0.0) sign = 1.0; + else sign = -1.0; + +// phi[0][i*NY+j] = 1.0; + /* nonconstant density to make things more interesting */ + phi[0][i*NY+j] = 0.5 + density_mod*module/amp; + phi[1][i*NY+j] = -sign*module*(xy[1]-y)/vabs(scale); + phi[2][i*NY+j] = sign*module*(xy[0]-x)/vabs(scale); + } + else + { + phi[0][i*NY+j] = 1.0; + phi[1][i*NY+j] = 0.0; + phi[2][i*NY+j] = 0.0; + } + } + break; + } + } } -void add_vortex_state(double x, double y, double scale, double *phi[NFIELDS], short int xy_in[NX*NY]) +void add_vortex_state(double amp, double x, double y, double scale, double density_mod, double *phi[NFIELDS], short int xy_in[NX*NY]) /* add vortex at position (x,y) with variance scale to field */ -/* phi[0] is stream function, phi[1] is vorticity */ +/* for incompressible Euler, phi[0] is stream function, phi[1] is vorticity */ +/* for compressible Euler, phi[0] is the density, phi[1] and phi[2] are velocity components */ { int i, j; double xy[2], dist2, module, phase, scale2, sign; // scale2 = scale*scale; - for (i=0; i 0.0) sign = 1.0; - else sign = -1.0; + if (scale > 0.0) sign = 1.0; + else sign = -1.0; - phi[1][i*NY+j] += sign*module; - phi[0][i*NY+j] -= sign*module; /* approximate, stream function should solve Poisson equation */ - } - else - { - phi[0][i*NY+j] = 0.0; - phi[1][i*NY+j] = 0.0; - } + phi[1][i*NY+j] -= sign*module; + phi[0][i*NY+j] = sign*module; /* approximate, stream function should solve Poisson equation */ + } + else + { + phi[0][i*NY+j] = 0.0; + phi[1][i*NY+j] = 0.0; + } + } + break; } + case (E_EULER_COMP): + { + for (i=0; i 0.0) sign = 1.0; + else sign = -1.0; + +// phi[0][i*NY+j] = 1.0; + /* nonconstant density to make things more interesting */ + phi[0][i*NY+j] += 0.5 + density_mod*sign*module/amp; + phi[1][i*NY+j] += sign*module*(xy[1]-y)/vabs(scale); + phi[2][i*NY+j] -= sign*module*(xy[0]-x)/vabs(scale); + } + else + { + phi[0][i*NY+j] = 1.0; + phi[1][i*NY+j] = 0.0; + phi[2][i*NY+j] = 0.0; + } + } + break; + } + } } @@ -421,17 +498,20 @@ void init_shear_flow(double amp, double delta, double rho, int nx, int ny, doubl phi[1][i*NY+j] = amp*(f - 1.0/(rho*cplus*cplus)); phi[0][i*NY+j] = amp*(f/(a*a) - rho/(cplus*cplus)); } + + phi[2][i*NY+j] = 0.0; } else { phi[0][i*NY+j] = 0.0; phi[1][i*NY+j] = 0.0; + phi[2][i*NY+j] = 0.0; } } } -void init_laminar_flow(double amp, double modulation, double period, double yshift, double *phi[NFIELDS], short int xy_in[NX*NY]) -/* initialise field with a laminar flow in x direction */ +void set_boundary_laminar_flow(double amp, double modulation, double period, double yshift, double *phi[NFIELDS], short int xy_in[NX*NY], int imin, int imax, int jmin, int jmax) +/* enfoce laminar flow in x direction on top and bottom boundaries */ /* phi[0] is stream function, phi[1] is vorticity */ /* amp is global amplitude */ { @@ -439,31 +519,121 @@ void init_laminar_flow(double amp, double modulation, double period, double yshi double xy[2], y1, a, b, f, cplus, cminus; a = period*PI/YMAX; -// a = PID/YMAX; - for (i=0; i 0)&&(maze[n].west)) + { + polyarc[na].xc = MAZE_XSHIFT; + polyarc[na].yc = 0.0; + polyarc[na].r = r; + polyarc[na].angle1 = phi; + polyarc[na].dangle = dphi; + polyarc[na].width = width; + na++; + } + + /* second circle */ + r = rmin + dr; + dphi *= 0.5; + for (q=0; q<2; q++) + { + n = nmaze(1, block*NXMAZE + q); + phi = (double)(block)*angle + (double)q*dphi; + + if (maze[n].west) + { + polyarc[na].xc = MAZE_XSHIFT; + polyarc[na].yc = 0.0; + polyarc[na].r = r; + polyarc[na].angle1 = phi; + polyarc[na].dangle = dphi; + polyarc[na].width = width; + na++; + } + } + + /* other circles */ + ww = 2; + i = 2; + while (ww < NXMAZE) + { + dphi *= 0.5; + for (p = 0; p < ww; p++) + { + r = rmin + (double)i*dr; + printf("Circle, i = %i, dphi = %.2lg, r = %.2lg\n", i, dphi, r); + for (q = 0; q < 2*ww; q++) + { + j = block*NXMAZE + q; + n = nmaze(i,j); + phi = (double)(block)*angle + (double)q*dphi; + + if (maze[n].west) + { + polyarc[na].xc = MAZE_XSHIFT; + polyarc[na].yc = 0.0; + polyarc[na].r = r; + polyarc[na].angle1 = phi; + polyarc[na].dangle = dphi; + polyarc[na].width = width; + na++; + } + } + i++; + } + ww *= 2; + } + } + + /* outer boundary of maze */ + polyarc[na].xc = MAZE_XSHIFT; + polyarc[na].yc = 0.0; + polyarc[na].r = rmax; + polyarc[na].angle1 = dphi; + polyarc[na].dangle = DPI - dphi; + polyarc[na].width = width; + na++; + + *npolyrect_rot = np; + *npolyarc = na; + + free(maze); +} + int init_polyline(int depth, t_vertex polyline[NMAXPOLY]) /* initialise variable polyline, for certain polygonal domain shapes */ { @@ -1820,6 +2045,10 @@ int init_polyrect(t_rectangle polyrect[NMAXPOLY]) { return(compute_maze_coordinates(polyrect, 1)); } + case (D_MAZE_CHANNELS): + { + return(compute_maze_coordinates(polyrect, 2)); + } default: { if ((ADD_POTENTIAL)&&(POTENTIAL == POT_MAZE)) return(compute_maze_coordinates(polyrect, 1)); @@ -1828,6 +2057,17 @@ int init_polyrect(t_rectangle polyrect[NMAXPOLY]) } } +void init_polyrect_arc(t_rect_rotated polyrectrot[NMAXPOLY], t_arc polyarc[NMAXPOLY], int *npolyrect, int *npolyarc) +/* initialise variables polyrectrot and polyarc, for certain domain shapes */ +{ + switch (B_DOMAIN) { + case (D_MAZE_CIRCULAR): + { + compute_circular_maze_coordinates(polyrectrot, polyarc, npolyrect, npolyarc); + break; + } + } +} void isospectral_initial_point(double x, double y, double left[2], double right[2]) /* compute initial coordinates in isospectral billiards */ @@ -1941,14 +2181,56 @@ int ij_in_polyrect(double i, double j, t_rectangle rectangle) return(1); } +int xy_in_rectrotated(double x, double y, t_rect_rotated rectrot) +/* returns 1 if (x,y) is in rectangle */ +{ + double l, u1, u2, v1, v2, pscal, h2; + + l = module2(rectrot.x2 - rectrot.x1, rectrot.y2 - rectrot.y1); + if (l == 0.0) return(0); + + /* unit vector along axis */ + u1 = (rectrot.x2 - rectrot.x1)/l; + u2 = (rectrot.y2 - rectrot.y1)/l; + + /* vector from one extremity to (x,y) */ + v1 = x - rectrot.x1; + v2 = y - rectrot.y1; + + /* inner product */ + pscal = u1*v1 + u2*v2; + if (pscal < 0.0) return(0); + if (pscal > l) return(0); + + h2 = v1*v1 + v2*v2 - pscal*pscal; + return(4.0*h2 <= rectrot.width*rectrot.width); +} + +int xy_in_arc(double x, double y, t_arc arc) +/* returns 1 if (x,y) is in arc */ +{ + double rho, phi, alpha; + + rho = module2(x - arc.xc, y - arc.yc); + + if (vabs(rho - arc.r) > 0.5*arc.width) return(0); + + phi = argument(x - arc.xc, y - arc.yc); + + alpha = phi - arc.angle1; + while (alpha < 0.0) alpha += DPI; + while (alpha > DPI) alpha -= DPI; + + return(alpha <= arc.dangle); +} 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, r, height; + 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, ca, sa, l, ht; int i, j, k, k1, k2, condition = 0, m; static int first = 1, nsides; - static double h, hh, ra, rb; + static double h, hh, ra, rb, ll, salpha; switch (b_domain) { case (D_NOTHING): @@ -2123,7 +2405,7 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ omega = DPI/((double)NPOLY); for (k=0; k polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(0); return(1); } + case (D_MAZE_CHANNELS): + { + 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_CIRCULAR): + { + for (i=0; i 5.0) return(1); return(0); } + case (D_FUNNELS): + { + y1 = y; + if (y > 0.5*YMAX) y1 -= YMAX; + if (y < -0.5*YMAX) y1 += YMAX; + y1 = vabs(y1 - MU*x); + y1 = vabs(y1 - 0.5*YMAX)*2.0/YMAX; + if (y1 > 0.25*(1.0 + LAMBDA + x*x)) return(0); + return(1); + } + case (D_ONE_FUNNEL): + { + y1 = vabs(y); + if (y1 > MU + LAMBDA*x*x) return(0); + return(1); + } + case (D_LENSES_RING): + { + if (first) + { + salpha = DPI/(double)NPOLY; + h = LAMBDA*tan(PI/(double)NPOLY); + if (h < MU) ll = sqrt(MU*MU - h*h); + else ll = 0.0; + first = 0; + } + for (i=0; i 0) draw_line(x1, y1, x, y); + x1 = x; + y1 = y; + } + } + } + break; + } + case (D_ONE_FUNNEL): + { + for (k=-1; k<2; k+=2) + { + x1 = XMIN; + y1 = (double)k*(MU + LAMBDA*x1*x1); + for (i=0; i<=NSEG; i++) + { + x = XMIN + (XMAX - XMIN)*(double)i/(double)NSEG; + y = (double)k*(MU + LAMBDA*x*x); + draw_line(x1, y1, x, y); + x1 = x; + y1 = y; + } + } + break; + } + case (D_LENSES_RING): + { + if (first) + { + salpha = DPI/(double)NPOLY; + h = LAMBDA*tan(PI/(double)NPOLY); + if (h < MU) ll = sqrt(MU*MU - h*h); + else ll = 0.0; + arcangle = atan(h/ll); + first = 0; + } + for (i=0; i 1.0) amp -= 2.0; - while (amp < -1.0) amp += 2.0; - amp_to_rgb(0.5*(1.0 + amp), rgb); +// value = min + 1.0*dy*(double)(j - jmin); +// amp = 0.7*color_amplitude_linear(value, 1.0, 1); +// while (amp > 1.0) amp -= 2.0; +// while (amp < -1.0) amp += 2.0; +// amp_to_rgb(0.5*(1.0 + amp), rgb); + value = dy_phase*(double)(j - jmin); + color_scheme_palette(C_ONEDIM_LINEAR, COLOR_PALETTE, value, 1.0, 1, rgb); break; } case (P_PHASE): @@ -4277,7 +4690,7 @@ void draw_color_scheme(double x1, double y1, double x2, double y2, int plot, dou void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int plot, double min, double max, int palette) { int j, k, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax; - double y, dy, dy_e, rgb[3], value, lum, amp; + double y, dy, dy_e, rgb[3], value, lum, amp, dy_phase; xy_to_ij(x1, y1, ij_botleft); xy_to_ij(x2, y2, ij_topright); @@ -4304,6 +4717,7 @@ void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int p glBegin(GL_QUADS); dy = (max - min)/((double)(jmax - jmin)); dy_e = max/((double)(jmax - jmin)); + dy_phase = 1.0/((double)(jmax - jmin)); for (j = jmin; j < jmax; j++) { @@ -4356,11 +4770,13 @@ void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int p } case (P_TOTAL_ENERGY_FLUX): { - value = min + 1.0*dy*(double)(j - jmin); - amp = 0.7*color_amplitude_linear(value, 1.0, 1); - while (amp > 1.0) amp -= 2.0; - while (amp < -1.0) amp += 2.0; - amp_to_rgb(0.5*(1.0 + amp), rgb); +// value = min + 1.0*dy*(double)(j - jmin); +// amp = 0.7*color_amplitude_linear(value, 1.0, 1); +// while (amp > 1.0) amp -= 2.0; +// while (amp < -1.0) amp += 2.0; +// amp_to_rgb(0.5*(1.0 + amp), rgb); + value = dy_phase*(double)(j - jmin); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; } case (P_PHASE): @@ -4404,7 +4820,7 @@ void draw_color_scheme_palette(double x1, double y1, double x2, double y2, int p void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, int plot, double min, double max, int palette, int fade, double fade_value) { int j, k, ij_botleft[2], ij_topright[2], imin, imax, jmin, jmax; - double y, dy, dy_e, rgb[3], value, lum, amp; + double y, dy, dy_e, rgb[3], value, lum, amp, dy_phase; xy_to_ij(x1, y1, ij_botleft); xy_to_ij(x2, y2, ij_topright); @@ -4431,6 +4847,7 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, glBegin(GL_QUADS); dy = (max - min)/((double)(jmax - jmin)); dy_e = max/((double)(jmax - jmin)); + dy_phase = 1.0/((double)(jmax - jmin)); for (j = jmin; j < jmax; j++) { @@ -4484,11 +4901,13 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, } case (P_TOTAL_ENERGY_FLUX): { - value = min + 1.0*dy*(double)(j - jmin); - amp = 0.7*color_amplitude_linear(value, 1.0, 1); - while (amp > 1.0) amp -= 2.0; - while (amp < -1.0) amp += 2.0; - amp_to_rgb(0.5*(1.0 + amp), rgb); +// value = min + 1.0*dy*(double)(j - jmin); +// amp = 0.7*color_amplitude_linear(value, 1.0, 1); +// while (amp > 1.0) amp -= 2.0; +// while (amp < -1.0) amp += 2.0; +// amp_to_rgb(0.5*(1.0 + amp), rgb); + value = dy_phase*(double)(j - jmin); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; } case (P_PHASE): @@ -4505,19 +4924,55 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, amp_to_rgb(0.5*(1.0 + amp), rgb); break; } - case (50): /* to fix: put #define in correct file */ + case (Z_EULER_VORTICITY): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } - case (51): + case (Z_EULER_LOG_VORTICITY): { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } - case (52): /* to fix: put #define in correct file */ + case (Z_EULER_VORTICITY_ASYM): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_LPRESSURE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_PRESSURE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_DENSITY): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULER_SPEED): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULERC_VORTICITY): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + default: { value = min + 1.0*dy*(double)(j - jmin); color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); diff --git a/sub_wave_3d.c b/sub_wave_3d.c index 9a1d923..47edb7a 100644 --- a/sub_wave_3d.c +++ b/sub_wave_3d.c @@ -961,9 +961,9 @@ void draw_billiard_3d_front(int fade, double fade_value) void compute_energy_field(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY]) /* computes cosine of angle between normal vector and vector light */ { - int i, j; + int i, j, k; static int first = 1; - double energy, logenergy; + double energy, logenergy, gx, gy, arg, mod, sum; // printf("computing energy field\n"); // printf("COMPUTE_MEAN_ENERGY = %i\n", COMPUTE_MEAN_ENERGY); @@ -1008,6 +1008,20 @@ void compute_energy_field(double phi[NX*NY], double psi[NX*NY], short int xy_in[ if (logenergy > LOG_ENERGY_FLOOR) wave[i*NY+j].log_mean_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*logenergy; else wave[i*NY+j].mean_energy = LOG_SHIFT + PLOT_SCALE_LOG_ENERGY*LOG_ENERGY_FLOOR; } + + if (COMPUTE_ENERGY_FLUX) + { + compute_energy_flux_mod(phi, psi, xy_in, i, j, &gx, &gy, &arg, &mod); + wave[i*NY+j].flux_direction = arg/DPI; + + /* compute time-averaged flux intensity */ + wave[i*NY+j].flux_int_table[wave[i*NY+j].flux_counter] = mod*FLUX_SCALE; + sum = 0.0; + for (k = 0; k < FLUX_WINDOW; k++) sum += wave[i*NY+j].flux_int_table[k]; + wave[i*NY+j].flux_intensity = sum/(double)FLUX_WINDOW; + wave[i*NY+j].flux_counter++; + if (wave[i*NY+j].flux_counter == FLUX_WINDOW) wave[i*NY+j].flux_counter = 0; + } } else if (first) { @@ -1016,6 +1030,8 @@ void compute_energy_field(double phi[NX*NY], double psi[NX*NY], short int xy_in[ wave[i*NY+j].log_total_energy = LOG_ENERGY_FLOOR; wave[i*NY+j].mean_energy = 0.0; wave[i*NY+j].log_mean_energy = LOG_ENERGY_FLOOR; + wave[i*NY+j].flux_intensity = 0.0; + wave[i*NY+j].flux_direction = 0.0; } } @@ -1099,9 +1115,12 @@ void compute_light_angle(short int xy_in[NX*NY], t_wave wave[NX*NY], int movie) } -void compute_field_color(double value, int cplot, int palette, double rgb[3]) +void compute_field_color(double value, double value2, int cplot, int palette, double rgb[3]) /* compute the color depending on the field value and color palette */ +/* value2 is only used for flux representation */ { + int k; + switch (cplot) { case (P_3D_AMP_ANGLE): { @@ -1110,7 +1129,7 @@ void compute_field_color(double value, int cplot, int palette, double rgb[3]) } case (P_3D_ENERGY): { - if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, VSCALE_ENERGY*value, 1.0, 0, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); break; } @@ -1121,7 +1140,7 @@ void compute_field_color(double value, int cplot, int palette, double rgb[3]) } case (P_3D_TOTAL_ENERGY): { - if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, VSCALE_ENERGY*value, 1.0, 0, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); break; } @@ -1132,7 +1151,7 @@ void compute_field_color(double value, int cplot, int palette, double rgb[3]) } case (P_3D_MEAN_ENERGY): { - if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); + if (COLOR_PALETTE >= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, VSCALE_ENERGY*value, 1.0, 0, rgb); else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 0, rgb); break; } @@ -1146,16 +1165,27 @@ void compute_field_color(double value, int cplot, int palette, double rgb[3]) amp_to_rgb_palette(value, rgb, palette); break; } + case (P_3D_FLUX_INTENSITY): + { + color_scheme_asym_palette(COLOR_SCHEME, palette, value*FLUX_SCALE, 1.0, 0, rgb); + break; + } + case (P_3D_FLUX_DIRECTION): + { + amp_to_rgb_palette(value, rgb, palette); + for (k=0; k<3; k++) rgb[k] *= tanh(value2*FLUX_CSCALE); + break; + } } } -double compute_interpolated_colors_wave(int i, int j, short int xy_in[NX*NY], t_wave wave[NX*NY], double palette, int cplot, - double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3], - int fade, double fade_value, int movie) +double compute_interpolated_colors_wave(int i, int j, short int xy_in[NX*NY], t_wave wave[NX*NY], + double palette, int cplot, double rgb_e[3], double rgb_w[3], double rgb_n[3], double rgb_s[3], int fade, double fade_value, int movie) { int k; double cw, ce, cn, cs, c_sw, c_se, c_nw, c_ne, c_mid, ca, z_mid; + double cw2, ce2, cn2, cs2; double *z_sw, *z_se, *z_nw, *z_ne; z_sw = wave[i*NY+j].p_zfield[movie]; @@ -1177,10 +1207,26 @@ double compute_interpolated_colors_wave(int i, int j, short int xy_in[NX*NY], t_ cs = (c_sw + c_se + c_mid)/3.0; cn = (c_nw + c_ne + c_mid)/3.0; - compute_field_color(ce, cplot, palette, rgb_e); - compute_field_color(cw, cplot, palette, rgb_w); - compute_field_color(cn, cplot, palette, rgb_n); - compute_field_color(cs, cplot, palette, rgb_s); + /* data for second color parameter */ + if (CHANGE_LUMINOSITY) + { + c_sw = *wave[i*NY+j].p_cfield[movie+2]; + c_se = *wave[(i+1)*NY+j].p_cfield[movie+2]; + c_nw = *wave[i*NY+j+1].p_cfield[movie+2]; + c_ne = *wave[(i+1)*NY+j+1].p_cfield[movie+2]; + + c_mid = 0.25*(c_sw + c_se + c_nw + c_ne); + + cw2 = (c_sw + c_nw + c_mid)/3.0; + ce2 = (c_se + c_ne + c_mid)/3.0; + cs2 = (c_sw + c_se + c_mid)/3.0; + cn2 = (c_nw + c_ne + c_mid)/3.0; + } + + compute_field_color(ce, ce2, cplot, palette, rgb_e); + compute_field_color(cw, cw2, cplot, palette, rgb_w); + compute_field_color(cn, cn2, cplot, palette, rgb_n); + compute_field_color(cs, cs2, cplot, palette, rgb_s); if (SHADE_3D) { @@ -1229,9 +1275,9 @@ void compute_wave_fields(double phi[NX*NY], double psi[NX*NY], short int xy_in[N void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc[NX*NY], double tgamma[NX*NY]) /* initialise fields for wave speed and dissipation */ { - int i, j, k; - double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB; - double u, v, u1, x, y, xy[2], norm2, speed, r2, c; + int i, j, k, inlens; + double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1; + double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1; if (VARIABLE_IOR) { @@ -1327,6 +1373,47 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc } break; } + case (IOR_EXPLO_LENSING): + { + salpha = DPI/(double)NPOLY; +// lambda1 = LAMBDA; +// mu1 = LAMBDA; + lambda1 = 0.5*LAMBDA; + mu1 = 0.5*LAMBDA; + h = lambda1*tan(PI/(double)NPOLY); + if (h < mu1) ll = sqrt(mu1*mu1 - h*h); + else ll = 0.0; + +// #pragma omp parallel for private(i,j) + for (i=0; i= COL_TURBO) color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + else color_scheme_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb); + break; + } + case (P_3D_FLUX_DIRECTION): + { + value = 2.0*dy_phase*(double)(j - jmin); + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } case (Z_EULER_VORTICITY): { value = min + 1.0*dy*(double)(j - jmin); @@ -1788,6 +1925,20 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } + case (Z_EULER_LPRESSURE): + { + value = min + 1.0*dy*(double)(j - jmin); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + case (Z_EULERC_VORTICITY): + { + value = min + 1.0*dy*(double)(j - jmin); + printf("Palette value %.3lg\n", value); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; + } + } if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value; glColor3f(rgb[0], rgb[1], rgb[2]); @@ -1837,3 +1988,18 @@ void print_speed_3d(double speed, int fade, double fade_value) write_text(xlefttext + 0.28, y, message); } +void init_wave_fields(t_wave wave[NX*NY]) +/* initialize some auxiliary fields */ +{ + int i, k; + + #pragma omp parallel for private(i) + for (i=0; i #include -#define MOVIE 0 /* set to 1 to generate movie */ +#define MOVIE 1 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ -#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ +#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ +#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ /* General geometrical parameters */ -#define WINWIDTH 1920 /* window width */ -#define WINHEIGHT 1000 /* window height */ -#define NX 2500 /* number of grid points on x axis */ -#define NY 1250 /* number of grid points on y axis */ -// #define NX 2700 /* number of grid points on x axis */ -// #define NY 1350 /* number of grid points on y axis */ +// #define WINWIDTH 1920 /* window width */ +// #define WINHEIGHT 1150 /* window height */ +// // // // #define NX 2500 /* number of grid points on x axis */ +// // // // #define NY 1250 /* number of grid points on y axis */ +// #define NX 1800 /* number of grid points on x axis */ +// #define NY 1800 /* number of grid points on y axis */ +// // #define NX 2700 /* number of grid points on x axis */ +// // #define NY 1350 /* number of grid points on y axis */ +// // #define YMIN -1.041666667 +// // #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ +// +// #define XMIN -2.0 +// #define XMAX 2.0 /* x interval */ // #define YMIN -1.041666667 // #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ -#define XMIN -2.0 -#define XMAX 2.0 /* x interval */ -#define YMIN -1.041666667 -#define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ +#define HIGHRES 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 1280 /* number of grid points on x axis */ -// // #define NY 720 /* number of grid points on y axis */ +// #define NX 1280 /* number of grid points on x axis */ +#define NX 720 /* number of grid points on x axis */ +#define NY 720 /* number of grid points on y axis */ // #define NX 2560 /* number of grid points on x axis */ +// #define NX 1440 /* number of grid points on x axis */ // #define NY 1440 /* number of grid points on y axis */ -// -// // #define NX 640 /* number of grid points on x axis */ -// // #define NY 360 /* number of grid points on y axis */ -// + +// #define NX 360 /* 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 XMIN -1.1 +#define XMAX 1.1 /* x interval */ +#define YMIN -1.1 +#define YMAX 1.1 /* y interval */ #define JULIA_SCALE 0.8 /* scaling for Julia sets */ /* Choice of the billiard table */ -#define B_DOMAIN 57 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 32 /* choice of domain shape, see list in global_pdes.c */ +// #define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 2 /* pattern of circles or polygons, see list in global_pdes.c */ @@ -93,18 +103,18 @@ #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 2 /* choice of index of refraction, see list in global_pdes.c */ +#define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */ +#define IOR 3 /* 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.25 /* parameter controlling the dimensions of domain */ -#define MU 0.0 /* parameter controlling the dimensions of domain */ -#define NPOLY 3 /* number of sides of polygon */ -#define APOLY 2.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define LAMBDA 0.75 /* parameter controlling the dimensions of domain */ +#define MU 0.25 /* 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 7 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ @@ -130,7 +140,7 @@ /* Physical parameters of wave equation */ -#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ #define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ #define OSCILLATION_SCHEDULE 3 /* oscillation schedule, see list in global_pdes.c */ @@ -140,8 +150,9 @@ #define AMPLITUDE 0.8 /* amplitude of periodic excitation */ #define ACHIRP 0.2 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.05 /* Courant number */ -#define COURANTB 0.1 /* Courant number in medium B */ +#define COURANT 0.1 /* Courant number */ +// #define COURANTB 0.0 /* Courant number in medium B */ +#define COURANTB 0.04658753 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMAB 0.0 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ @@ -155,8 +166,8 @@ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ #define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ -#define OSCILLATING_SOURCE_PERIOD 100 /* period of oscillating source */ -#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 100 /* period of oscillating source */ +#define ALTERNATE_OSCILLATING_SOURCE 0 /* set to 1 to alternate sign of oscillating source */ /* Boundary conditions, see list in global_pdes.c */ @@ -167,16 +178,16 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1400 /* number of frames of movie */ -// #define NSTEPS 250 /* number of frames of movie */ -#define NVID 5 /* number of iterations between images displayed on screen */ +#define NSTEPS 2800 /* number of frames of movie */ +// #define NSTEPS 200 /* number of frames of movie */ +#define NVID 6 /* 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 */ #define PRINT_SPEED 0 /* set to 1 to print speed of moving source */ #define PAUSE 100 /* number of frames after which to pause */ -#define PSLEEP 2 /* sleep time during pause */ +#define PSLEEP 3 /* sleep time during pause */ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ #define MID_FRAMES 100 /* number of still frames between parts of two-part movie */ @@ -186,31 +197,41 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 0.25 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */ +// #define INITIAL_AMP 0.25 /* amplitude of initial condition */ +// #define INITIAL_VARIANCE 0.001 /* variance of initial condition */ +// #define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ +#define INITIAL_AMP 0.3 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0004 /* variance of initial condition */ #define INITIAL_WAVELENGTH 0.02 /* wavelength of initial condition */ +// #define INITIAL_AMP 0.35 /* amplitude of initial condition */ +// #define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ +// #define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -#define ZPLOT 103 /* wave height */ -#define CPLOT 103 /* color scheme */ +// #define ZPLOT 103 /* wave height */ +// #define CPLOT 103 /* color scheme */ +#define ZPLOT 104 /* wave height */ +#define CPLOT 104 /* color scheme */ -#define ZPLOT_B 107 -#define CPLOT_B 107 /* plot type for second movie */ +#define ZPLOT_B 108 +#define CPLOT_B 108 /* plot type for second movie */ +#define CHANGE_LUMINOSITY 1 /* set to 1 to let luminosity depend on energy flux intensity */ +#define FLUX_WINDOW 30 /* size of averaging window of flux intensity */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of plot */ #define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ #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 1 /* set to 1 to draw front of boundary after drawing wave */ +#define DRAW_BILLIARD 0 /* set to 1 to draw boundary */ +#define DRAW_BILLIARD_FRONT 0 /* 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.01 /* vertical scaling in energy plot */ -#define PLOT_SCALE_LOG_ENERGY 0.25 /* vertical scaling in log energy plot */ +#define PLOT_SCALE_LOG_ENERGY 0.2 /* vertical scaling in log energy plot */ /* 3D representation */ @@ -220,30 +241,32 @@ #define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */ #define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ -#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ +#define ROTATE_ANGLE 540.0 /* total angle of rotation during simulation */ +// #define ROTATE_ANGLE 45.0 /* total angle of rotation during simulation */ /* Color schemes */ -#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 11 /* 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 */ #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 1.0 /* sensitivity of color on wave amplitude */ +#define SLOPE 0.5 /* sensitivity of color on wave amplitude */ #define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ -#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 VSCALE_ENERGY 30.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 100.0 /* scaling factor for energy representation */ +#define E_SCALE 100.0 /* scaling factor for energy representation */ #define LOG_SCALE 0.75 /* scaling factor for energy log representation */ #define LOG_SHIFT 0.5 /* shift of colors on log scale */ #define LOG_ENERGY_FLOOR -10.0 /* floor value for log of (total) energy */ #define LOG_MEAN_ENERGY_SHIFT 1.0 /* additional shift for log of mean energy */ -#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */ +#define FLUX_SCALE 20.0 /* scaling factor for energy flux representation */ +#define FLUX_CSCALE 500.0 /* scaling factor for color in energy flux representation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -253,15 +276,15 @@ #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 */ +#define NXMAZE 8 /* width of maze */ +#define NYMAZE 32 /* height of maze */ +#define MAZE_MAX_NGBH 5 /* max number of neighbours of maze cell */ +#define RAND_SHIFT 5 /* 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.5 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 3.0 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE 6.0 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 6.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 */ @@ -283,11 +306,11 @@ 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, 8.0}; /* location of observer for REP_PROJ_3D representation */ +double observer[3] = {8.0, 8.0, 12.0}; /* location of observer for REP_PROJ_3D representation */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ -#define Z_SCALING_FACTOR 0.3 /* overall scaling factor of z axis for REP_PROJ_3D representation */ -#define XY_SCALING_FACTOR 2.4 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define Z_SCALING_FACTOR 0.25 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define XY_SCALING_FACTOR 2.0 /* 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 */ @@ -892,9 +915,12 @@ void evolve_wave(double phi[NX*NY], double psi[NX*NY], double tmp[NX*NY], short void draw_color_bar_palette(int plot, double range, int palette, int fade, double fade_value) { - double width = 0.14; + double width = 0.1; +// double width = 0.14; // double width = 0.2; + width *= (double)NX/(double)WINWIDTH; + if (ROTATE_COLOR_SCHEME) draw_color_scheme_palette_3d(-1.0, -0.8, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); else @@ -926,8 +952,9 @@ void viewpoint_schedule(int i) 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, *psi, *tmp, *total_energy, *color_scale, *tc, *tcc, *tgamma; + double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle, lambda1; + double *phi, *psi, *tmp, *color_scale, *tc, *tcc, *tgamma; +// double *total_energy; short int *xy_in; int i, j, s, sample_left[2], sample_right[2], period = 0, fade; static int counter = 0; @@ -946,7 +973,7 @@ void animation() phi = (double *)malloc(NX*NY*sizeof(double)); psi = (double *)malloc(NX*NY*sizeof(double)); tmp = (double *)malloc(NX*NY*sizeof(double)); - total_energy = (double *)malloc(NX*NY*sizeof(double)); +// total_energy = (double *)malloc(NX*NY*sizeof(double)); color_scale = (double *)malloc(NX*NY*sizeof(double)); tc = (double *)malloc(NX*NY*sizeof(double)); tcc = (double *)malloc(NX*NY*sizeof(double)); @@ -985,6 +1012,10 @@ void animation() npolyrect = init_polyrect(polyrect); for (i=0; i= INITIAL_TIME) save_frame(); +// if (i >= INITIAL_TIME) save_frame_counter(i); +// if (i >= INITIAL_TIME) save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); else printf("Initial phase time %i of %i\n", i, INITIAL_TIME); if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) @@ -1135,8 +1178,10 @@ void animation() if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0); glutSwapBuffers(); save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter); +// save_frame_counter(i); counter++; } + else if (NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); /* it seems that saving too many files too fast can cause trouble with the file system */ /* so this is to make a pause from time to time - parameter PAUSE may need adjusting */ @@ -1172,9 +1217,15 @@ void animation() draw_wave_3d(0, phi, psi, xy_in, wave, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value); if (PRINT_SPEED) print_speed_3d(speed, 1, fade_value); - glutSwapBuffers(); + if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers(); save_frame_counter(NSTEPS + i + 1); } + + if ((ROTATE_VIEW)&&(ROTATE_VIEW_WHILE_FADE)) + { + viewpoint_schedule(NSTEPS - INITIAL_TIME); + reset_view = 1; + } draw_wave_3d(1, phi, psi, xy_in, wave, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, 1); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0); @@ -1201,6 +1252,11 @@ void animation() if (!FADE) for (i=0; i #include -#define MOVIE 0 /* set to 1 to generate movie */ +#define MOVIE 1 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ -#define SAVE_MEMORY 0 /* set to 1 to save memory when writing tiff images */ +#define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ +#define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ /* 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 x axis */ +#define WINHEIGHT 1150 /* window height */ +// // #define NX 1920 /* number of grid points on x 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 */ - +#define NY 2300 /* 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 XMAX 2.0 /* x interval */ +#define YMIN -1.197916667 +#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ +// #define YMIN -1.041666667 +// #define YMAX 1.041666667 /* y interval for 9/16 aspect ratio */ +// #define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ // #define WINWIDTH 1280 /* window width */ @@ -82,7 +85,7 @@ /* Choice of the billiard table */ -#define B_DOMAIN 57 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 999 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 1 /* pattern of circles or polygons, see list in global_pdes.c */ @@ -94,8 +97,8 @@ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.25 /* parameter controlling the dimensions of domain */ -#define MU 0.0 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.5 /* parameter controlling the dimensions of domain */ +#define MU 0.5 /* 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 */ @@ -122,7 +125,7 @@ /* Physical parameters of wave equation */ -#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ #define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ #define OSCILLATION_SCHEDULE 1 /* oscillation schedule, see list in global_pdes.c */ @@ -131,8 +134,8 @@ #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.1 /* Courant number in medium B */ +#define COURANT 0.08 /* Courant number */ +#define COURANTB 0.04658753 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMAB 0.0 /* damping factor in wave equation */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ @@ -146,8 +149,8 @@ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ #define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ -#define OSCILLATING_SOURCE_PERIOD 50 /* period of oscillating source */ -#define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 15 /* period of oscillating source */ +#define ALTERNATE_OSCILLATING_SOURCE 0 /* set to 1 to alternate sign of oscillating source */ /* Boundary conditions, see list in global_pdes.c */ @@ -155,9 +158,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1700 /* number of frames of movie */ -// #define NSTEPS 3500 /* number of frames of movie */ -#define NVID 12 /* number of iterations between images displayed on screen */ +#define NSTEPS 2500 /* number of frames of movie */ +// #define NSTEPS 500 /* number of frames of movie */ +#define NVID 12 /* number of iterations between images displayed on screen */ +// #define NVID 9 /* 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 INITIAL_TIME 400 /* time after which to start saving frames */ @@ -174,35 +178,38 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 0.35 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.0002 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.01 /* wavelength of initial condition */ +#define INITIAL_AMP 0.1 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0003 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.015 /* wavelength of initial condition */ +// #define INITIAL_VARIANCE 0.00015 /* variance of initial condition */ +// #define INITIAL_WAVELENGTH 0.0075 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -#define PLOT 0 -// #define PLOT 1 +#define PLOT 7 +// #define PLOT 7 -#define PLOT_B 5 /* plot type for second movie */ +#define PLOT_B 6 /* plot type for second movie */ /* Color schemes */ -#define COLOR_PALETTE 13 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 18 /* Color palette, see list in global_pdes.c */ +// #define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 11 /* 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 1.0 /* sensitivity of color on wave amplitude */ +#define SLOPE 0.75 /* sensitivity of color on wave amplitude */ #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 100.0 /* scaling factor for energy representation */ -#define LOG_SCALE 0.25 /* scaling factor for energy log representation */ -#define LOG_SHIFT 0.0 /* shift of colors on log scale */ -#define FLUX_SCALE 1.0e4 /* scaling factor for enegy flux represtnation */ +#define E_SCALE 140.0 /* scaling factor for energy representation */ +#define LOG_SCALE 1.0 /* scaling factor for energy log representation */ +#define LOG_SHIFT 3.5 /* shift of colors on log scale */ +#define FLUX_SCALE 2.0e3 /* scaling factor for enegy flux represtnation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -213,16 +220,16 @@ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 7.0 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE 1.5 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #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 NXMAZE 8 /* width of maze */ +#define NYMAZE 32 /* height of maze */ +#define MAZE_MAX_NGBH 5 /* max number of neighbours of maze cell */ +#define RAND_SHIFT 0 /* seed of random number generator */ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ /* for compatibility with sub_wave and sub_maze */ @@ -231,7 +238,6 @@ #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 */ @@ -571,10 +577,10 @@ 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 time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle, x1, sign1; double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX], *total_flux; short int *xy_in[NX]; - int i, j, s, sample_left[2], sample_right[2], period = 0, fade; + int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0; static int counter = 0; long int wave_value; @@ -608,7 +614,12 @@ void animation() npolyrect = init_polyrect(polyrect); for (i=0; i XMIN)&&(x1 < XMAX)) + { + add_circular_wave(sign1, x1, y, phi, psi, xy_in); + printf("Adding wave at (%.2lg, %.2lg)\n", x1, y); + } + sign1 = -sign1; + } + source_counter++; + if (source_counter == 4) + { + source_counter = 0; + sign = -sign; + } +// for (j=0; j max) gradientx = max; + else if (gradientx < -max) gradientx = -max; + if (gradienty > max) gradienty = max; + else if (gradienty < -max) gradienty = -max; + + current_mod = velocity*module2(gradientx, gradienty); + if (current_mod > 1.0e-10) + { + current_arg = argument(gradientx,gradienty); + if (current_arg < 0.0) current_arg += DPI; + if (current_arg >= DPI) current_arg -= DPI; + } + else current_arg = PI; + *gx = velocity*gradientx; + *gy = velocity*gradienty; + } + else + { + current_mod = 0.0; + current_arg = PI; + *gx = 0.0; + *gy = 0.0; + } + + *module = current_mod; + *arg = current_arg; +} + + double compute_phase(double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], int i, int j) { double velocity, angle;