diff --git a/colors_waves.c b/colors_waves.c index dc0ad4c..0b658f0 100644 --- a/colors_waves.c +++ b/colors_waves.c @@ -63,6 +63,14 @@ void color_scheme(int scheme, double value, double scale, int time, double rgb[3 amp_to_rgb(amplitude, rgb); break; } + case (C_BASIC_LINEAR): + { + amplitude = color_amplitude_linear(value, scale, time); + if (amplitude > 1.0) amplitude -= 1.0; + else if (amplitude < 0.0) amplitude += 1.0; + amp_to_rgb(amplitude, rgb); + break; + } } } @@ -186,6 +194,17 @@ void color_scheme_palette(int scheme, int palette, double value, double scale, i amp_to_rgb_palette(amplitude, rgb, palette); break; } + case (C_BASIC_LINEAR): + { + if (value > 1.0) value -= 1.0; + else if (value < 0.0) value += 1.0; + amp_to_rgb_palette(value, rgb, palette); + break; +// if (value > 2.0) value -= 2.0; +// else if (value < 0.0) value += 2.0; +// amp_to_rgb_palette(value, rgb, palette); +// break; + } } } diff --git a/drop_billiard.c b/drop_billiard.c index e82ec1f..5a143f1 100644 --- a/drop_billiard.c +++ b/drop_billiard.c @@ -49,6 +49,7 @@ #define POLYLINE_PATTERN 1 /* pattern of polyline */ #define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ +#define NABSCIRCLES 10 /* number of absorbing circles */ #define NMAXCIRCLES 1000 /* total number of circles (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 1000 /* total number of sides of polygonal line */ @@ -70,6 +71,9 @@ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NPOLY 4 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define LAMBDA_B 1.0 /* parameter controlling shape of domain (for P_POLYRING) */ +#define NPOLY_B 100000 /* number of sides of second polygon */ +#define APOLY_B 1.0 /* angle by which to turn second polygon, in units of Pi/2 */ #define DRAW_BILLIARD 1 /* set to 1 to draw billiard */ #define DRAW_CONSTRUCTION_LINES 1 /* set to 1 to draw additional construction lines for billiard */ #define PERIODIC_BC 0 /* set to 1 to enforce periodic boundary conditions when drawing particles */ diff --git a/global_ljones.c b/global_ljones.c index 77c522b..6e96b10 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -23,7 +23,7 @@ #define NMAXBELTS 10 /* max number of conveyor belts */ #define NMAXSHOVELS 50 /* max number of shovels */ #define NMAX_TRIANGLES_PER_OBSTACLE 10 /* max number of triangles per obstacle */ -#define NMAX_TRIANGLES_PER_FACET 5 /* max number of triangles per facet between obstacles */ +#define NMAX_TRIANGLES_PER_FACET 4 /* max number of triangles per facet between obstacles */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -388,7 +388,8 @@ #define C_HUE 1 /* color scheme modifies hue */ #define C_PHASE 2 /* color scheme shows phase */ #define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma, Twilight) */ -#define C_ONEDIM_LINEAR 4 /* use preset 1d color scheme with linear scale */ +#define C_ONEDIM_LINEAR 4 /* use preset 1d color scheme with linear scale */ +#define C_BASIC_LINEAR 5 /* use preset 1d color scheme with linear scale */ /* Color palettes */ diff --git a/global_particles.c b/global_particles.c index 319be34..d0cda13 100644 --- a/global_particles.c +++ b/global_particles.c @@ -92,6 +92,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define P_TOKA_PRIME 6 /* Tokarsky room made of 86 triangles */ #define P_TREE 7 /* pine tree */ #define P_TOKA_NONSELF 8 /* Tokarsky non-self-unilluminable room */ +#define P_ISOCELES_TRIANGLE 9 /* rectangle isoceles triangle */ #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 */ @@ -99,6 +100,7 @@ double x_shooter = -0.2, y_shooter = -0.6, x_target = 0.4, y_target = 0.7; #define P_MAZE_CIRC_SCATTERER 14 /* circular maze with scatterers */ #define P_MAZE_HEX 15 /* hexagonal maze */ #define P_MAZE_OCT 16 /* maze with octagonal and square cells */ +#define P_STAR 17 /* star-shaped domain */ /* Color palettes */ diff --git a/global_pdes.c b/global_pdes.c index e49f323..05f26aa 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -115,6 +115,11 @@ #define D_CARDIOID 90 /* cardioid */ #define D_NEPHROID 91 /* nephroid */ #define D_EPICYCLOID 92 /* epicycloid */ +#define D_CIRCLE_LATTICE 93 /* lattice of connected circles */ +#define D_CIRCLE_LATTICE_RANDOM 931 /* lattice of connected circles with random channel widths */ +#define D_CIRCLE_LATTICE_HEX 94 /* hex lattice of connected circles */ +#define D_CIRCLE_LATTICE_HONEY 95 /* honeycomb lattice of connected circles */ +#define D_CIRCLE_LATTICE_POISSON 96 /* Poisson disc process of connected circles */ /* for wave_sphere.c */ @@ -144,6 +149,7 @@ #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 NMAXLINES 1000 /* maximal number of lines (for bounadries) */ #define NMAXSOURCES 30 /* maximal number of sources */ #define C_SQUARE 0 /* square grid of circles */ @@ -304,7 +310,9 @@ #define C_HUE 1 /* color scheme modifies hue */ #define C_PHASE 2 /* color scheme shows phase */ #define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma) */ -#define C_ONEDIM_LINEAR 4 /* use preset 1d color scheme with linear scale */ +#define C_ONEDIM_LINEAR 4 /* use preset 1d color scheme with linear scale */ +#define C_BASIC_LINEAR 5 /* use preset 1d color scheme with linear scale and slope 1, range is assumed to be [0,1] */ + /* Color palettes */ @@ -330,6 +338,7 @@ /* plot types used by rde */ #define Z_AMPLITUDE 0 /* amplitude of first field */ +#define Z_ANGLE 10 /* angle, for Kuramoto model */ #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 */ @@ -403,6 +412,11 @@ typedef struct double posi, posj; /* (i,j) coordinates of vertex */ } t_vertex; +// typedef struct +// { +// t_vertex z1, z2; /* extremities of line */ +// } t_line; + typedef struct { double x1, y1, x2, y2; /* (x,y) coordinates of vertices */ @@ -457,6 +471,7 @@ 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 */ +// int nlines = NMAXLINES; /* actual number of lines */ short int input_signal[NSTEPS]; /* time-dependent source signal */ @@ -466,6 +481,7 @@ 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 */ +// t_line line[NMAXLINES]; /* data of lines */ /* the same for comparisons between different domains */ int ncircles_b = NMAXCIRCLES; /* actual number of circles, can be decreased e.g. for random patterns */ diff --git a/global_perc.c b/global_perc.c index e2f0e7e..5216202 100644 --- a/global_perc.c +++ b/global_perc.c @@ -42,6 +42,7 @@ #define C_PHASE 2 /* color scheme shows phase */ #define C_ONEDIM 3 /* use preset 1d color scheme (for Turbo, Viridis, Magma, Inferno, Plasma) */ #define C_ONEDIM_LINEAR 4 /* use preset 1d color scheme with linear scale */ +#define C_BASIC_LINEAR 5 /* use preset 1d color scheme with linear scale */ /* Color palettes */ diff --git a/heat.c b/heat.c index 9b3db9d..960211e 100644 --- a/heat.c +++ b/heat.c @@ -222,6 +222,8 @@ #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ #define WALL_WIDTH 0.1 /* width of wall separating lenses */ +#define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ +#define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */ #define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define INITIAL_TIME 50 /* time after which to start saving frames */ #define OSCIL_YMAX 0.35 /* defines oscillation range */ diff --git a/lennardjones.c b/lennardjones.c index e9fc290..1867002 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -37,7 +37,7 @@ #include #define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 1 /* set to 1 to save memory while saving frames */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -58,25 +58,25 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.0 +#define INITXMIN -1.95 #define INITXMAX -1.9 /* x interval for initial condition */ -#define INITYMIN -0.1 -#define INITYMAX 0.1 /* y interval for initial condition */ +#define INITYMIN 0.9 +#define INITYMAX 1.0 /* y interval for initial condition */ #define THERMOXMIN -1.25 #define THERMOXMAX 1.25 /* x interval for initial condition */ #define THERMOYMIN 0.0 #define THERMOYMAX 0.75 /* y interval for initial condition */ -#define ADDXMIN -2.2 -#define ADDXMAX -2.0 /* x interval for adding particles */ -#define ADDYMIN -1.0 -#define ADDYMAX 1.0 /* y interval for adding particles */ +#define ADDXMIN -1.95 +#define ADDXMAX -1.9 /* x interval for adding particles */ +#define ADDYMIN 0.0 +#define ADDYMAX 0.6 /* y interval for adding particles */ #define ADDRMIN 4.75 #define ADDRMAX 6.0 /* r interval for adding particles */ -#define BCXMIN -2.5 -#define BCXMAX 2.5 /* x interval for boundary condition */ +#define BCXMIN -2.0 +#define BCXMAX 2.0 /* x interval for boundary condition */ #define BCYMIN -1.125 #define BCYMAX 1.125 /* y interval for boundary condition */ @@ -90,15 +90,15 @@ #define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */ #define CIRCLE_PATTERN_B 0 /* pattern of circles for additional particles */ -#define ADD_FIXED_OBSTACLES 1 /* set to 1 do add fixed circular obstacles */ -#define OBSTACLE_PATTERN 20 /* pattern of obstacles, see list in global_ljones.c */ +#define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */ +#define OBSTACLE_PATTERN 91 /* pattern of obstacles, see list in global_ljones.c */ #define RATTLE_OBSTACLES 0 /* set to 1 to rattle obstacles (for pattern O_SIEVE_B) */ #define OSCILLATE_OBSTACLES 1 /* set to 1 to make obstacles oscillate */ #define COUPLE_OBSTACLES 1 /* set to 1 to couple obstacles to neighbours */ -#define OBSTACLE_PISC_DISTANCE 0.12 /* minimal distance in Poisson disc process for obstacles, controls density of obstacles */ -#define OBSTACLE_COUPLING_DIST 0.18 /* max distance of coupled obstacles */ +#define OBSTACLE_PISC_DISTANCE 0.08 /* minimal distance in Poisson disc process for obstacles, controls density of obstacles */ +#define OBSTACLE_COUPLING_DIST 0.12 /* max distance of coupled obstacles */ #define NMAX_OBSTACLE_NEIGHBOURS 8 /* max number of obstacle neighbours */ -#define NMAX_OBSTACLE_PINNED 10 /* max number of neighbours to be pinned */ +#define NMAX_OBSTACLE_PINNED 3 /* max number of neighbours to be pinned */ #define OBSTACLE_PINNING_TYPE 0 /* type of obstacle pinning, see OP_ in global_ljones */ #define BDRY_PINNING_STEP 4 /* interval between pinned obstacles on boundary */ #define RECOUPLE_OBSTACLES 0 /* set to 1 to reset obstacle coupling */ @@ -126,7 +126,7 @@ #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 12 /* 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 4.0 /* angular frequency of spin-spin interaction */ #define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */ @@ -134,15 +134,16 @@ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define PDISC_DISTANCE 1.0 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 5.0 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ -#define MU 0.007 /* parameter controlling radius of particles */ +#define MU 0.02 /* parameter controlling radius of particles */ #define MU_B 0.03 /* parameter controlling radius of particles of second type */ +#define MU_ADD 0.02 /* parameter controlling radius of added particles */ #define NPOLY 3 /* number of sides of polygon */ -#define APOLY 0.075 /* angle by which to turn polygon, in units of Pi/2 */ +#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define AWEDGE 0.5 /* opening angle of wedge, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ @@ -155,8 +156,8 @@ #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 NOBSX 20 -#define NOBSY 10 /* obstacles for O_HEX obstacle pattern */ +#define NOBSX 40 +#define NOBSY 24 /* obstacles for O_HEX obstacle pattern */ #define NTREES 15 /* number of trees in S_TREES */ #define X_SHOOTER -0.2 @@ -166,10 +167,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 4500 /* number of frames of movie */ -#define NVID 4500 /* number of iterations between images displayed on screen */ +#define NSTEPS 2200 /* number of frames of movie */ +#define NVID 150 /* number of iterations between images displayed on screen */ #define NSEG 25 /* number of segments of boundary of circles */ -#define INITIAL_TIME 0 /* time after which to start saving frames */ +#define INITIAL_TIME 30 /* time after which to start saving frames */ #define OBSTACLE_INITIAL_TIME 0 /* time after which to start moving obstacle */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ @@ -184,25 +185,26 @@ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 23 +#define BOUNDARY_COND 1 /* Plot type, see list in global_ljones.c */ -#define PLOT 13 -#define PLOT_B 12 /* plot type for second movie */ +#define PLOT 17 +// #define PLOT 23 +#define PLOT_B 17 /* plot type for second movie */ /* Background color depending on particle properties */ -#define COLOR_BACKGROUND 0 /* set to 1 to color background */ -#define BG_COLOR 7 /* type of background coloring, see list in global_ljones.c */ -#define BG_COLOR_B 5 /* type of background coloring, see list in global_ljones.c */ +#define COLOR_BACKGROUND 1 /* set to 1 to color background */ +#define BG_COLOR 4 /* type of background coloring, see list in global_ljones.c */ +#define BG_COLOR_B 2 /* type of background coloring, see list in global_ljones.c */ #define OBSTACLE_COLOR 0 /* type of obstacle, see OC_ in global_ljones.c */ #define DRAW_BONDS 0 /* set to 1 to draw bonds between neighbours */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ #define FILL_TRIANGLES 0 /* set to 1 to fill triangles between neighbours */ #define DRAW_CLUSTER_LINKS 0 /* set to 1 to draw links between particles in cluster */ -#define DRAW_OBSTACLE_LINKS 1 /* set to 1 to draw links between interacting obstacles */ +#define DRAW_OBSTACLE_LINKS 0 /* set to 1 to draw links between interacting obstacles */ #define FILL_OBSTACLE_TRIANGLES 0 /* set to 1 to fill triangles between interacting obstacles */ #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 */ @@ -210,8 +212,8 @@ #define INITIAL_POS_TYPE 0 /* type of initial position dependence */ #define ERATIO 0.995 /* ratio for time-averaging in P_EMEAN color scheme */ #define DRATIO 0.999 /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */ -#define OBSTACLE_AREA_SHADE_FACTOR 160.0 /* controls sensitivity of triangle shade for option FILL_OBSTACLE_TRIANGLES */ -#define SHADE_OBSTACLE_FACETS 0 /* set to 1 to shade facets instead of triangles */ +#define OBSTACLE_AREA_SHADE_FACTOR 80.0 /* controls sensitivity of triangle shade for option FILL_OBSTACLE_TRIANGLES */ +#define SHADE_OBSTACLE_FACETS 1 /* set to 1 to shade facets instead of triangles */ /* Color schemes */ @@ -263,10 +265,10 @@ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ #define PARTICLE_EMIN 100.0 /* energy of particle with coolest color */ -#define PARTICLE_EMAX 3000.0 /* energy of particle with hottest color */ +#define PARTICLE_EMAX 100000.0 /* energy of particle with hottest color */ #define SEGMENT_HUE_MIN 275.0 /* color of original segment */ #define SEGMENT_HUE_MAX 30.0 /* color of saturated segment */ -#define OBSTACLE_EMAX 200.0 /* energy of obstacle with hottest color */ +#define OBSTACLE_EMAX 150.0 /* energy of obstacle with hottest color */ #define OBSTACLE_VMAX 4.0 /* speed of obstacle with largest luminosity */ #define HUE_TYPE0 320.0 /* hue of particles of type 0 */ #define HUE_TYPE1 60.0 /* hue of particles of type 1 */ @@ -276,7 +278,7 @@ #define HUE_TYPE5 60.0 /* hue of particles of type 5 */ #define HUE_TYPE6 130.0 /* hue of particles of type 6 */ #define HUE_TYPE7 150.0 /* hue of particles of type 7 */ -#define BG_FORCE_SLOPE 7.5e-8 /* contant in BG_FORCE backgound color scheme*/ +#define BG_FORCE_SLOPE 1.0e-6 /* contant in BG_FORCE backgound color scheme*/ #define RANDOM_RADIUS 0 /* set to 1 for random particle radius */ #define RANDOM_RADIUS_MIN 0.4 /* min of random particle radius (default 0.75) */ @@ -284,9 +286,9 @@ #define ADAPT_MASS_TO_RADIUS 1 /* set to positive value to for mass prop to power of radius */ #define ADAPT_DAMPING_TO_RADIUS 0.5 /* set to positive value to for friction prop to power of radius */ #define ADAPT_DAMPING_FACTOR 0.5 /* factor by which damping is adapted to radius */ -#define DT_PARTICLE 2.0e-7 /* time step for particle displacement */ +#define DT_PARTICLE 1.5e-6 /* time step for particle displacement */ #define KREPEL 50.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 2.5 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST 2.0 /* Lennard-Jones equilibrium distance */ #define EQUILIBRIUM_DIST_B 2.5 /* Lennard-Jones equilibrium distance for second type of particle */ #define SEGMENT_FORCE_EQR 1.0 /* equilibrium distance factor for force from segments (default 1.5) */ #define REPEL_RADIUS 25.0 /* radius in which repelling force acts (in units of particle radius) */ @@ -294,10 +296,12 @@ #define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */ #define DAMPING_ROT 5.0 /* damping coefficient for rotation of particles */ #define PARTICLE_MASS 2.0 /* mass of particle of radius MU */ -#define PARTICLE_MASS_B 2.0 /* mass of particle of radius MU_B */ +#define PARTICLE_MASS_B 2.0 /* mass of particle of radius MU_B */ +#define PARTICLE_ADD_MASS 2.0 /* mass of added particles */ #define PARTICLE_INERTIA_MOMENT 0.5 /* moment of inertia of particle */ #define PARTICLE_INERTIA_MOMENT_B 0.5 /* moment of inertia of second type of particle */ -#define V_INITIAL 50.0 /* initial velocity range */ +#define V_INITIAL 200.0 /* initial velocity range */ +#define V_INITIAL_ADD 50.0 /* initial velocity range for added particles */ #define OMEGA_INITIAL 100.0 /* initial angular velocity range */ #define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */ #define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */ @@ -324,13 +328,14 @@ #define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */ #define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */ -#define ADD_EFIELD 1 /* set to 1 to add an electric field */ -#define EFIELD 0.0 /* value of electric field */ -#define EFIELD_Y 1500.0 /* value of electric field */ -#define ADD_BFIELD 1 /* set to 1 to add a magnetic field */ -#define BFIELD 1200.0 /* value of magnetic field */ -#define CHARGE 1.0 /* charge of particles of first type */ +#define ADD_EFIELD 0 /* set to 1 to add an electric field */ +#define EFIELD -100.0 /* value of electric field */ +#define EFIELD_Y 0.0 /* value of electric field */ +#define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ +#define BFIELD 0.0 /* value of magnetic field */ +#define CHARGE 1.0 /* charge of particles of first type */ #define CHARGE_B 1.0 /* charge of particles of second type */ +#define CHARGE_ADD 1.0 /* charge of added particles */ #define INCREASE_E 0 /* set to 1 to increase electric field */ #define EFIELD_FACTOR 5000000.0 /* factor by which to increase electric field */ #define INCREASE_B 0 /* set to 1 to increase magnetic field */ @@ -338,8 +343,8 @@ #define CHARGE_OBSTACLES 0 /* set to 1 for obstacles to be charged */ #define OBSTACLE_CHARGE 3.0 /* charge of obstacles */ #define OBSTACLE_MASS 100.0 /* mass of obstacles, if oscillating */ -#define KSPRING_OBSTACLE_OSC 1.0e7 /* spring constant for oscillating obstacles */ -#define KSPRING_OBSTACLE_COUPLE 5.0e6 /* spring constant for coupled obstacles */ +#define KSPRING_OBSTACLE_OSC 5.0e5 /* spring constant for oscillating obstacles */ +#define KSPRING_OBSTACLE_COUPLE 2.0e5 /* spring constant for coupled obstacles */ #define OBSTACLE_HARDCORE 1 /* set to 1 to add "hard core" repulsion between obstacles */ #define KSPRING_OBSTACLE_HARDCORE 1.0e11 /* spring constant for obstacle hard core repulsion */ #define KCOULOMB_OBSTACLE 1000.0 /* Coulomb force constant for charged obstacles */ @@ -359,8 +364,8 @@ #define KTORQUE_DIFF 500.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 0 /* set to 1 to draw cross on particles of second type */ -#define DRAW_MINUS 0 /* set to 1 to draw cross on particles of negative charge */ +#define DRAW_CROSS 1 /* set to 1 to draw cross on particles of second type */ +#define DRAW_MINUS 1 /* set to 1 to draw cross on particles of negative charge */ #define SPIN_RANGE 10.0 /* range of spin-spin interaction */ #define SPIN_RANGE_B 10.0 /* range of spin-spin interaction for second type of particle */ #define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */ @@ -386,7 +391,7 @@ #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.03 /* radius of obstacle for circle boundary conditions */ +#define OBSTACLE_RADIUS 0.015 /* 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 */ @@ -413,16 +418,17 @@ #define ADD_PARTICLES 1 /* set to 1 to add particles */ #define ADD_REGION 0 /* shape of add regions, cf ADD_* in global_ljones */ #define ADD_TIME 0 /* time at which to add first particle */ -#define ADD_PERIOD 150 /* time interval between adding further particles */ +#define ADD_PERIOD 35 /* time interval between adding further particles */ #define N_ADD_PARTICLES 1 /* number of particles to add */ -#define FINAL_NOADD_PERIOD 500 /* final period where no particles are added */ +#define FINAL_NOADD_PERIOD 250 /* final period where no particles are added */ #define SAFETY_FACTOR 4.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ +#define ADD_ALTERNATE_CHARGE 1 /* set to 1 to randomly select sign of added charge */ #define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */ -#define N_TRACER_PARTICLES 100 /* number of tracer particles */ +#define N_TRACER_PARTICLES 200 /* number of tracer particles */ #define TRACER_STEPS 5 /* number of tracer steps recorded between images */ -#define TRAJECTORY_LENGTH 40000 /* length of recorded trajectory */ -#define TRACER_LUM_FACTOR 100.0 /* controls luminosity decrease of trajectories with time */ +#define TRAJECTORY_LENGTH 5000 /* length of recorded trajectory */ +#define TRACER_LUM_FACTOR 40.0 /* controls luminosity decrease of trajectories with time */ #define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */ #define TRAJECTORY_WIDTH 2 /* width of tracer particle trajectory */ @@ -567,7 +573,7 @@ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 29 /* size of hashgrid in x direction */ +#define HASHX 40 /* size of hashgrid in x direction */ #define HASHY 20 /* 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 */ @@ -2368,7 +2374,7 @@ void animation() // draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, conveyor_belt, wall); /* add a particle */ - if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) + if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i - INITIAL_TIME < NSTEPS - FINAL_NOADD_PERIOD)) { /* add enzymes */ if ((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR))) diff --git a/mangrove.c b/mangrove.c index aab03ee..d475dbd 100644 --- a/mangrove.c +++ b/mangrove.c @@ -266,6 +266,8 @@ #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ #define DRAW_WAVE_TIMESERIES 0 /* set to 1 to draw a time series of the wave */ #define WALL_WIDTH 0.1 /* width of wall separating lenses */ +#define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ +#define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */ #define OSCIL_YMAX 0.35 /* defines oscillation range */ #define MESSAGE_LDASH 14 /* length of dash for Morse code message */ #define MESSAGE_LDOT 8 /* length of dot for Morse code message */ diff --git a/particle_billiard.c b/particle_billiard.c index 7f109d0..6668d18 100644 --- a/particle_billiard.c +++ b/particle_billiard.c @@ -35,20 +35,20 @@ #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 INVERT_COUNTER 0 /* set to 1 to save frames in inverse order */ // #define WINWIDTH 1280 /* window width */ -#define WINWIDTH 720 /* window width */ -#define WINHEIGHT 720 /* window height */ +#define WINWIDTH 1200 /* window width */ +#define WINHEIGHT 1200 /* window height */ // #define XMIN -1.5 // #define XMAX 2.5 /* x interval */ -#define XMIN -1.125 -#define XMAX 1.125 /* x interval */ -#define YMIN -1.125 -#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ +#define XMIN -1.2 +#define XMAX 1.4 /* x interval */ +#define YMIN -1.4 +#define YMAX 1.2 /* y interval for 9/16 aspect ratio */ #define SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */ @@ -57,9 +57,10 @@ #define B_DOMAIN 30 /* choice of domain shape */ #define CIRCLE_PATTERN 1 /* pattern of circles */ -#define POLYLINE_PATTERN 10 /* pattern of polyline */ +#define POLYLINE_PATTERN 9 /* pattern of polyline */ -#define ABSORBING_CIRCLES 0 /* set to 1 for circular scatterers to be absorbing */ +#define ABSORBING_CIRCLES 1 /* set to 1 for circular scatterers to be absorbing */ +#define NABSCIRCLES 10 /* number of absorbing circles */ #define NMAXCIRCLES 100000 /* total number of circles (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 100000 /* total number of sides of polygonal line */ @@ -70,10 +71,13 @@ #define SDEPTH 1 /* Sierpinski gastket depth */ #define LAMBDA 1.5 /* parameter controlling shape of domain */ -#define MU 0.005 /* second parameter controlling shape of billiard */ +#define MU 0.02 /* 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 */ +#define LAMBDA_B 1.0 /* parameter controlling shape of domain (for P_POLYRING) */ +#define NPOLY_B 100000 /* number of sides of second polygon */ +#define APOLY_B 1.0 /* angle by which to turn second polygon, in units of Pi/2 */ #define PENROSE_RATIO 2.5 /* parameter controlling the shape of small ellipses in Penrose room */ #define DRAW_BILLIARD 1 /* set to 1 to draw billiard */ @@ -85,21 +89,21 @@ /* Simulation parameters */ -// #define NPART 10 /* number of particles */ -#define NPART 50000 /* number of particles */ +#define NPART 21 /* number of particles */ +// #define NPART 50000 /* 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 0 /* set to 1 to keep trails of the particles */ -#define HEATMAP 1 /* set to 1 to show heat map of particles */ -#define DRAW_FINAL_HEATMAP 1 /* set to 1 to show final heat map of particles */ +#define SHOWTRAILS 1 /* set to 1 to keep trails of the particles */ +#define HEATMAP 0 /* set to 1 to show heat map of particles */ +#define DRAW_FINAL_HEATMAP 0 /* set to 1 to show final heat map of particles */ #define DRAW_HEATMAP_HISTOGRAM 0 /* set to 1 to draw a histogram of particle distribution in heat map */ #define NBIN_FACTOR 6.0 /* constant controlling number of bins in histogram */ -#define DRAW_HEATMAP_PARTICLES 1 /* set to 1 to draw particles in heat map */ +#define DRAW_HEATMAP_PARTICLES 0 /* set to 1 to draw particles in heat map */ #define HEATMAP_MAX_PART_BY_CELL 50 /* set to positive value to draw only limited number of particles in cell */ -#define PLOT_HEATMAP_AVERAGE 1 /* 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 PLOT_HEATMAP_AVERAGE 0 /* set to 1 to plot average number of particles in heat map */ +#define SHOWZOOM 1 /* 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 0 /* 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 */ @@ -108,11 +112,11 @@ #define TEST_INITIAL_COND 0 /* set to 1 to allow only initial conditions that pass a test */ -#define NSTEPS 1300 /* number of frames of movie */ -#define TIME 3000 /* time between movie frames, for fluidity of real-time simulation */ +#define NSTEPS 6000 /* number of frames of movie */ +#define TIME 1000 /* time between movie frames, for fluidity of real-time simulation */ // #define DPHI 0.000002 /* integration step */ -#define DPHI 0.00002 /* integration step */ -#define NVID 25 /* number of iterations between images displayed on screen */ +#define DPHI 0.000005 /* integration step */ +#define NVID 50 /* number of iterations between images displayed on screen */ #define END_FRAMES 50 /* number of still frames at the end of the movie */ /* Decreasing TIME accelerates the animation and the movie */ @@ -132,7 +136,7 @@ #define RAINBOW_COLOR 1 /* 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.025 /* length of velocity vectors */ +#define LENGTH 0.01 /* 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 */ @@ -328,9 +332,14 @@ void draw_zoom(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTM glEnd(); /* draw billiard boundaries in zoom */ - glLineWidth(BILLIARD_WIDTH*2); + glLineWidth(BILLIARD_WIDTH*10); - if (y_target + width > 1.0) + if (POLYLINE_PATTERN == P_ISOCELES_TRIANGLE) + { + draw_line(shiftx, shifty, x1, shifty); + draw_line(shiftx, shifty, x1, y2); + } + else if (y_target + width > 1.0) { yb = shifty + 0.5*(1.0 - y_target)/width; glBegin(GL_LINE_STRIP); @@ -352,6 +361,7 @@ void draw_zoom(int color[NPARTMAX], double *configs[NPARTMAX], int active[NPARTM // glLineWidth(PARTICLE_WIDTH*2); + /* draw particles */ for (i=0; i +#include +#include +#include +#include +#include +#include /* Sam Leffler's libtiff library. */ +#include + +#define MOVIE 0 /* set to 1 to generate movie */ + +#define WINWIDTH 1760 /* window width */ +#define WINHEIGHT 990 /* window height */ + +#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 */ + +#define XMIN -1.2 +#define XMAX 2.8 /* x interval */ +#define YMIN -1.125 +#define YMAX 1.125 /* y interval for 9/16 aspect ratio */ + +#define SCALING_FACTOR 1.0 /* scaling factor of drawing, needed for flower billiards, otherwise set to 1.0 */ + +/* Choice of the billiard table, see global_particles.c */ + +#define B_DOMAIN 30 /* choice of domain shape */ + +#define CIRCLE_PATTERN 1 /* pattern of circles */ +#define POLYLINE_PATTERN 8 /* pattern of polyline */ + +#define ABSORBING_CIRCLES 1 /* set to 1 for circular scatterers to be absorbing */ +#define NABSCIRCLES 10 /* number of absorbing circles */ + +#define NMAXCIRCLES 50000 /* total number of circles (must be at least NCX*NCY for square grid) */ +#define NMAXPOLY 100004 /* total number of sides of polygonal line */ +#define NCX 9 /* number of circles in x direction */ +#define NCY 20 /* number of circles in y direction */ +#define NPOISSON 500 /* number of points for Poisson C_RAND_POISSON arrangement */ +#define NGOLDENSPIRAL 2000 /* max number of points for C_GOLDEN_SPIRAL arrandement */ +#define SDEPTH 2 /* Sierpinski gastket depth */ + +#define LAMBDA 1.0 /* parameter controlling shape of domain */ +#define MU 0.0075 /* second parameter controlling shape of billiard */ +#define FOCI 1 /* set to 1 to draw focal points of ellipse */ +#define NPOLY 5 /* number of sides of polygon */ +#define APOLY 0.2 /* angle by which to turn polygon, in units of Pi/2 */ +#define LAMBDA_B 1.0 /* parameter controlling shape of domain (for P_POLYRING) */ +#define NPOLY_B 100000 /* number of sides of second polygon */ +#define APOLY_B 1.0 /* angle by which to turn second polygon, in units of Pi/2 */ +#define DRAW_BILLIARD 1 /* set to 1 to draw billiard */ +#define DRAW_CONSTRUCTION_LINES 0 /* set to 1 to draw additional construction lines for billiard */ +#define PERIODIC_BC 0 /* set to 1 to enforce periodic boundary conditions when drawing particles */ +#define PENROSE_RATIO 2.5 /* parameter controlling the shape of small ellipses in Penrose room */ + +#define RESAMPLE 0 /* set to 1 if particles should be added when dispersion too large */ +#define DEBUG 0 /* draw trajectories, for debugging purposes */ + +/* Simulation parameters */ + +#define NPART 1 /* number of particles */ +#define NPARTMAX 100000 /* maximal number of particles after resampling */ +#define TRAJ_LENGTH 120 /* length of trajectory */ +#define PLOT_NMAX 100 /* max length on length plot */ +#define LMAX 0.01 /* minimal segment length triggering resampling */ +#define DMIN 0.02 /* minimal distance to boundary for triggering resampling */ +#define CYCLE 0 /* set to 1 for closed curve (start in all directions) */ +#define SHOWTRAILS 0 /* set to 1 to keep trails of the particles */ +#define HEATMAP 0 /* set to 1 to show heat map of particles */ +#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 TEST_ACTIVE 0 /* set to 1 to test whether particle is in billiard */ +#define PRINT_TRAJECTORY_LENGTH 1 /* set to 1 to print length of trajectory 0 */ +#define PRINT_TRAJECTORY_ANGLE 1 /* set to 1 to print angle of trajectory 0 */ +#define PRINT_TRAJECTORY_SLOPE 0 /* set to 1 to print slope of trajectory 0 */ +#define PRINT_TRAJECTORY_PERIOD 0 /* set to 1 to print period of trajectory 0 */ +#define DRAW_LENGTHS_PLOT 1 /* set to 1 to plot trajectory lengths */ +#define LENGTHS_LOG_SCALE 0 /* set to 1 to use log scale for plot of lengths */ +#define LENGTH_PLOT_POLAR 0 /* set to 1 to plot lengths in polar coordinates */ +#define MIN_ANGLE 0.0 /* range of angles of trajectory */ +#define MAX_ANGLE 135.0 /* range of angles of trajectory */ +#define TEST_INITIAL_COND 0 /* set to 1 to allow only initial conditions that pass a test */ + +#define SLOW_AT_LONG_TRAJ 0 /* set to 1 to slow down movie for long trajectories */ +#define ADD_SUCCESS_GALLERY 0 /* set to 1 to add gallery of successful trajectories at end of movie */ +#define SUCCESS_GALLERY_FRAMES 25 /* number of frames per success */ +#define EXIT_BOTH_WAYS 1 /* set to 1 to add exits to he left to succesful trajectories */ + +#define NSTEPS 6400 /* number of frames of movie */ +#define TIME 2500 /* time between movie frames, for fluidity of real-time simulation */ +#define DPHI 0.00001 /* integration step */ +#define NVID 150 /* number of iterations between images displayed on screen */ + +/* Decreasing TIME accelerates the animation and the movie */ +/* For constant speed of movie, TIME*DPHI should be kept constant */ +/* However, increasing DPHI too much deterioriates quality of simulation */ +/* NVID tells how often a picture is drawn in the animation, increase it for faster anim */ +/* For a good quality movie, take for instance TIME = 400, DPHI = 0.00005, NVID = 100 */ + +/* Colors and other graphical parameters */ + +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ + +#define NCOLORS 64 /* number of colors */ +#define COLORSHIFT 0 /* 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 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.03 /* 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 */ +#define GRAPH_HUE 200.0 /* color hue of graph */ +#define SUCCESS_HUE 30.0 /* color hue of success */ + +#define BLACK 1 /* set to 1 for black background */ +#define COLOR_OUTSIDE 0 /* set to 1 for colored outside */ +#define OUTER_COLOR 270.0 /* color outside billiard */ +#define PAINT_INT 0 /* set to 1 to paint interior in other color (for polygon/Reuleaux) */ +#define PAINT_EXT 1 /* set to 1 to paint exterior */ + +#define PAUSE 1000 /* number of frames after which to pause */ +#define PSLEEP 2 /* sleep time during pause */ +#define SLEEP1 1 /* initial sleeping time */ +#define SLEEP2 1 /* final sleeping time */ +#define END_FRAMES 250 /* number of frames at end of movie */ + +#define NXMAZE 25 /* width of maze */ +#define NYMAZE 25 /* height of maze */ +#define MAZE_MAX_NGBH 8 /* max number of neighbours of maze cell */ +#define RAND_SHIFT 10 /* seed of random number generator */ +#define MAZE_XSHIFT -0.7 /* horizontal shift of maze */ +#define MAZE_RANDOM_FACTOR 0.1 /* randomization factor for S_MAZE_RANDOM */ +#define MAZE_CORNER_RADIUS 0.4 /* radius of tounded corners in maze */ +#define CLOSE_MAZE 0 /* set to 1 to close maze exits */ + +#include "global_particles.c" +#include "sub_maze.c" +#include "sub_part_billiard.c" + + +/*********************/ +/* animation part */ +/*********************/ + +void init_boundary_config(double smin, double smax, double anglemin, double anglemax, double *configs[NPARTMAX]) +/* initialize configuration: drop on the boundary, beta version */ +/* WORKS FOR ELLIPSE, HAS TO BE ADAPTED TO GENERAL BILLIARD */ +{ + int i; + double ds, da, s, angle, theta, alpha, pos[2]; + + if (anglemin <= 0.0) anglemin = PI/((double)NPART); + if (anglemax >= PI) anglemax = PI*(1.0 - 1.0/((double)NPART)); + ds = (smax - smin)/((double)NPART); + da = (anglemax - anglemin)/((double)NPART); + for (i=0; i 1) dalpha = (angle2 - angle1)/((double)(NPART-1)); + else dalpha = 0.0; + for (i=0; i 1) dalpha = (angle2 - angle1)/((double)(particle2 - particle1-1)); + else dalpha = 0.0; + for (i=particle1; i=1; j--) + { + lum = 0.5*(1.0 - (double)j/(double)TRAJ_LENGTH); + + for (i=0; i 1.0) + { + yb = shifty + 0.5*(1.0 - y_target)/width; + glBegin(GL_LINE_STRIP); + glVertex2d(x1, yb); + glVertex2d(x2, yb); + glVertex2d(x2, yb + 0.02); + glVertex2d(x1, yb + 0.02); + glEnd(); + } + /* other boundaries not yet implemented */ + + /* draw target in zoom */ + glLineWidth(BILLIARD_WIDTH*2); + + if (shooter) glColor3f(1.0, 0.0, 0.0); + else glColor3f(0.0, 0.8, 0.2); + + glBegin(GL_LINE_LOOP); + tradius = zoomwidth*MU/width; + for (i=0; i<=NSEG; i++) + { + phi = (double)i*DPI/(double)NSEG; + x1 = shiftx + tradius*cos(phi); + y1 = shifty + tradius*sin(phi); + glVertex2d(x1, y1); + } + glEnd (); + +// glLineWidth(PARTICLE_WIDTH*2); +} + +void draw_zoom_area(double x_target, double y_target, double width) +/* draw rectangle marking zoomed area */ +{ + double x1, x2, y1, y2; + + glEnable(GL_LINE_SMOOTH); + glLineWidth(BILLIARD_WIDTH/2); + + x1 = x_target - width; + y1 = y_target - width; + x2 = x_target + width; + y2 = y_target + width; + + glColor3f(1.0, 1.0, 1.0); + glBegin(GL_LINE_LOOP); + glVertex2d(x1, y1); + glVertex2d(x2, y1); + glVertex2d(x2, y2); + glVertex2d(x1, y2); + glEnd(); +} + +void draw_trajectories_lshape(int color[NPARTMAX], double *trajectory[NPARTMAX], int traj_length[NPARTMAX]) +/* draw the trajectories */ +{ + int i, j; + double rgb[3], lum, x1, y1, x2, y2, x3, y3, x4, y4; + + for (j=TRAJ_LENGTH-2; j>=1; j--) + { + lum = 0.5*(1.0 - (double)j/(double)TRAJ_LENGTH); + + for (i=0; i= 0.0) x3 = 0.0; + else x3 = 1.0; + y3 = y1; + } + else if ((x1 == 0.0)||(x1 == 1.0)) + { + x3 = -1.0; + y3 = y1; + } + if (y1 == -1.0) + { + if (x1 >= 0.0) y3 = 0.0; + else y3 = 1.0; + x3 = x1; + } + else if ((y1 == 0.0)||(y1 == 1.0)) + { + y3 = -1.0; + x3 = x1; + } + + if (j < traj_length[i]) + { + rgb_color_scheme_lum(color[i], lum, rgb); + glColor3f(rgb[0], rgb[1], rgb[2]); +// glBegin(GL_LINE_STRIP); +// glVertex2d(x1, y1); +// glVertex2d(x2, y2); +// glEnd (); + glBegin(GL_LINE_STRIP); + glVertex2d(x2, y2); + glVertex2d(x3, y3); + glEnd (); + } + } + } + } + + +void draw_trajectories(int color[NPARTMAX], double *trajectory[NPARTMAX], int traj_length[NPARTMAX]) +/* draw the trajectories */ +{ + int i, j, col; + double rgb[3], lum, x1, y1, x2, y2; + + glutSwapBuffers(); + if (!SHOWTRAILS) blank(); + if (PAINT_INT) paint_billiard_interior(); + + if (SHOWZOOM) switch (POLYLINE_PATTERN) { + case (P_TOKA_PRIME): + { + draw_zoom(color, trajectory, traj_length, x_target, y_target, 0.05, 1.65, 0.75, 0.3, 0); + draw_zoom(color, trajectory, traj_length, x_shooter, y_shooter, 0.05, -1.65, 0.75, 0.3, 1); + break; + } + case (P_TOKA_NONSELF): + { + draw_zoom(color, trajectory, traj_length, 0.0, 0.0, 0.1, 1.65, 0.75, 0.3, 1); + break; + } + case (P_ISOCELES_TRIANGLE): + { + draw_zoom(color, trajectory, traj_length, 1.0, -1.0, 0.2, 0.6, 0.6, 0.4, 0); + break; + } + } + + glLineWidth(PARTICLE_WIDTH); + + glEnable(GL_LINE_SMOOTH); + + /* if domain is L shape with periodic boundary conditions, we have to take care of those */ + if (B_DOMAIN == D_LSHAPE) draw_trajectories_lshape(color, trajectory, traj_length); + + else for (j=TRAJ_LENGTH-1; j>=1; j--) + { + lum = 0.5*(1.0 - 0.5*(double)j/(double)TRAJ_LENGTH); + + for (i=0; i=1; j--) if (j < traj_length[i]) + { + x1 = trajectory[i][2*j]; + y1 = trajectory[i][2*j+1]; + x2 = trajectory[i][2*j-2]; + y2 = trajectory[i][2*j-1]; + + len = module2(x2 - x1, y2 - y1); + + for (k=0; k<(int)(len/pathstep); k++) + { + x = x1 + (x2 - x1)*(double)k*pathstep/len; + y = y1 + (y2 - y1)*(double)k*pathstep/len; + + 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]++; + } + +// if (HEATMAP_MAX_PART_BY_CELL > 0) +// { +// drawtable[i] = ((n == -2)||((n > 0)&&(heatmap_number[n] <= HEATMAP_MAX_PART_BY_CELL))); +// } +// else + drawtable[i] = (n > -1); + } + } + } + + } + + for (n=0; n= 100) shiftx = 0.1; + else if (j >= 10) shiftx = 0.0775; + else shiftx = 0.055; + write_text_fixedwidth(xmin - shiftx, y - 0.01, message); + } + + /* horizontal axis graduation */ + for (j=0; j iw) x -= c*(double)(i - iw); + if ((x > xmin)&&(x < xmax - 0.05)) + { + y = ymin - 0.025; + y1 = ymin + 0.025; + draw_line(x, y, x, y1); + if (x < xmax - 0.15) + { + sprintf(message, "%i", jmin + j); + if (j < 100) x += 0.01; + if (j < 10) x += 0.015; + write_text_fixedwidth(x - 0.04, ymin - 0.075, message); + } + } + } + + /* to reduce aliasing, plot a broader graph at smaller luminosity */ + glLineWidth(3); + if (TEST_INITIAL_COND) + { + glColor3f(rgb1_success[0], rgb1_success[1], rgb1_success[2]); + for (j=imin; jimin) + { + draw_line(x0, y0, x, y0); + draw_line(x, y0, x, y); + } + + x0 = x; + y0 = y; + } + + glLineWidth(2); + x0 = xmin; + y0 = ymin; + + if (TEST_INITIAL_COND) + { + glColor3f(rgb_success[0], rgb_success[1], rgb_success[2]); + for (j=imin; jimin) + { + draw_line(x0, y0, x, y0); + draw_line(x, y0, x, y); + } + + x0 = x; + y0 = y; + } + + draw_initial_condition_circle(x0, y0, 0.02, 30); + +// draw_initial_condition_circle(x0, y0, 0.02, traj_length_table[i]); + +} + +void plot_lengths_polar(int traj_length_table[NSTEPS+1], short int traj_test_table[NSTEPS+1], int i) +/* draw a plot of thajectory lengths */ +{ + int j, k; + double xc, yc, x0, y0, x1, y1, xmin = 0.5, xmax = 1.9, ymin = -0.7, ymax = 0.3, phi, dphi, r; + char message[50]; + static int first = 1; + static double rgb[3], rgb_success[3], logfactor; + + if (first) + { + rgb_color_scheme(GRAPH_HUE, rgb); + rgb_color_scheme(SUCCESS_HUE, rgb_success); + logfactor = 0.07; + } + + dphi = DPI/(double)NSTEPS; + xc = 0.5*(xmin + xmax); + yc = 0.5*(ymin + ymax); + + glColor3f(1.0, 1.0, 1.0); + glLineWidth(1); + + /* draw grid */ + r = 0.5*(xmax - xmin); + for (j=0; j<12; j++) + { + phi = (double)j*DPI/12.0; + draw_line(xc, yc, xc + r*cos(phi), yc + r*sin(phi)); + } + for (j=1; j<5; j++) + { + r = logfactor*(double)j*log(10.0); + draw_circle(xc, yc, r, NSEG); + } + + glLineWidth(2); + /* draw plot of succesful paths */ +// if (TEST_INITIAL_COND) +// { +// r = 0.5*(xmax - xmin); +// glColor3f(rgb_success[0], rgb_success[1], rgb_success[2]); +// for (j=0; j 1) dalpha = (angle2 - angle1)/((double)(NPART)); + else dalpha = 0.0; + for (i=0; i= 0.0)&&(traj_length[i] < TRAJ_LENGTH - 1)) + { +// c = vbilliard(configs[i]); + pos_billiard(configs[i], pos, &alpha); + c = vbilliard_xy(configs[i], alpha, pos); + j = traj_length[i]; + trajectory[i][2*j] = configs[i][4]; + trajectory[i][2*j+1] = configs[i][5]; + traj_length[i]++; + + if ((!loop)&&(i == 0)) + { + if (module2(configs[i][0] - s0, configs[i][1] - u0) > 0.001) period++; + else loop = 1; + } + } + if (i==0) c0 = c; + j = traj_length[i]; + trajectory[i][2*j] = configs[i][6]; + trajectory[i][2*j+1] = configs[i][7]; + traj_length[i]++; + + } + + *xmax = trajectory[0][2*(traj_length[0]-1)]; + +// printf("Period = %i\n", period); + if ((c0 == DUMMY_SIDE_ABS)||(c0 < 0)) return(0); + else return(period + 1); +// if ((c != DUMMY_SIDE_ABS)&&(c >= 0.0)) return (period); +// else return(0); +} + +void compute_trajectories_ellipse(double x, double y, double *configs[NPARTMAX], + double *trajectory[NPARTMAX], int traj_length[NPARTMAX]) +/* compute trajectories when starting in (x,y), with angles between angle1 and angle2 */ +{ + int i, j, c; + double dalpha, alpha, a, b, cc, dx, d; + double conf[2], pos[2]; + + cc = sqrt(LAMBDA*LAMBDA - 1.0); + dx = 4.0*cc/(double)NPART; + + for (i=0; i XMAX)||((EXIT_BOTH_WAYS)&&(xmax < XMIN))) + { + traj_test_table[i] = 1; + nsuccess++; + } + else traj_test_table[i] = 0; +// traj_test_table[i] = test_initial_condition(configs, active, newcolor); +// if (traj_test_table[i] >= 1) nsuccess++; + } + +// printf("Period = %i\n", period); + + draw_trajectories(color, trajectory, traj_length); + if (HEATMAP) draw_config_heatmap(trajectory, traj_length, configs, active, heatmap_number, heatmap_total, 0); +// else draw_trajectories(color, trajectory, traj_length); + + print_parameters(i, traj_length, traj_length_table, period, alpha); + if (DRAW_BILLIARD) draw_billiard(); + if (DRAW_LENGTHS_PLOT) plot_lengths(traj_length_table, traj_test_table, i+1); + + /* draw the shooter in red */ + rgb[0] = 1.0; rgb[1] = 0.0; rgb[2] = 0.0; + draw_colored_circle(x, y, 0.01, NSEG, rgb); +// draw_colored_circle(x_shooter, y_shooter, 0.01, NSEG, rgb); + + for (j=0; j= 10)) + { + npause = (int)(0.5*sqrt((double)traj_length[0])); + for (j=0; j0)&&(traj_test_table[i])) + { +// sleep(0.5); + /* additional buffer swap seems to be necessary on some OS */ + glutSwapBuffers(); + if (first == 0) + { +// gallery_counter = NSTEPS + END_FRAMES + success_counter*SUCCESS_GALLERY_FRAMES; + glutSwapBuffers(); + first = 1; + } + else + { + for (j=0; j= nend) nend = gallery_counter - NSTEPS; + for (i=0; i 0)) + x += NU_KURAMOTO*delta_phi[0][i*NY+j]; + + x = phi_in[0][i*NY+j] + intstep*x; + if (x > PI) + { + y = (x + PI)/DPI; + phi_out[0][i*NY+j] = x - DPI*(double)((int)y); + } + else if (x < -PI) + { + y = (-x + PI)/DPI; + phi_out[0][i*NY+j] = x + DPI*(double)((int)y); + } + else phi_out[0][i*NY+j] = x; + break; +// x = phi_in[0][i*NY+j]; +// phi_out[0][i*NY+j] = phi_in[0][i*NY+j] + intstep*(K_KURAMOTO*kuramoto_int[i*NY+j] + OMEGA_KURAMOTO); +// if (phi_out[0][i*NY+j] > PI) +// { +// x = (phi_out[0][i*NY+j] + PI)/DPI; +// phi_out[0][i*NY+j] -= DPI*(double)((int)x); +// } +// else if (phi_out[0][i*NY+j] < -PI) +// { +// x = (-phi_out[0][i*NY+j] + PI)/DPI; +// phi_out[0][i*NY+j] += DPI*(double)((int)x); +// } +// break; + } case (E_CAHN_HILLIARD): { /* TO DO */ @@ -1681,6 +1729,10 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) free(delta_v); } } + else if (RDE_EQUATION == E_KURAMOTO) + { + free(kuramoto_int); + } if (COMPUTE_PRESSURE) { @@ -2111,7 +2163,7 @@ void animation() xy_in = (short int *)malloc(NX*NY*sizeof(short int)); rde = (t_rde *)malloc(NX*NY*sizeof(t_rde)); - if (SPHERE) +// if (SPHERE) { wsphere = (t_wave_sphere *)malloc(NX*NY*sizeof(t_wave_sphere)); init_wave_sphere_rde(wsphere,1); @@ -2243,8 +2295,11 @@ void animation() // init_tidal_state(1, 0.0015, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); // init_linear_blob_sphere(0, 1.3*PI, 0.65*PI, 0.0, 0.0, 0.3, 0.1, 0.1, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); + + init_random(0.0, 0.4, phi, xy_in, wsphere); +// init_random_smoothed(0.0, 0.4, phi, xy_in, wsphere); - init_expanding_blob_sphere(0, (279.0/180.0)*PI, (115.0/180.0)*PI, 0.75, 0.04, 0.06, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); +// init_expanding_blob_sphere(0, (279.0/180.0)*PI, (115.0/180.0)*PI, 0.75, 0.04, 0.06, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); // add_gaussian_wave(-1.6, -0.5, 0.015, 0.25, SWATER_MIN_HEIGHT, phi, xy_in); @@ -2507,7 +2562,7 @@ void animation() free(phi_tmp[i]); } free(xy_in); - if (SPHERE) +// if (SPHERE) { free(wsphere); free(wsphere_hr); diff --git a/schrodinger.c b/schrodinger.c index eb550bf..5a6b960 100644 --- a/schrodinger.c +++ b/schrodinger.c @@ -202,6 +202,8 @@ #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ #define WALL_WIDTH 0.1 /* width of wall separating lenses */ +#define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ +#define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */ #define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define INITIAL_TIME 50 /* time after which to start saving frames */ #define OSCIL_YMAX 0.35 /* defines oscillation range */ diff --git a/sub_lj.c b/sub_lj.c index 089222e..af43730 100644 --- a/sub_lj.c +++ b/sub_lj.c @@ -2018,8 +2018,8 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle otrian } case (O_SQUARE_TWOMASSES): { - dx = (XMAX - XMIN)/(double)NOBSX; - dy = (YMAX - YMIN)/(double)NOBSY; + dx = (OBSXMAX - OBSXMIN)/(double)NOBSX; + dy = (OBSYMAX - OBSYMIN)/(double)NOBSY; n = 0; if ((ADD_FIXED_SEGMENTS)&&(SEGMENT_PATTERN == S_CYLINDER)) @@ -2037,9 +2037,9 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle otrian for (i=0; i XMAX) obstacle[n].xc += (XMAX - XMIN); + obstacle[n].xc = OBSXMIN + ((double)i + 0.25)*dx; + obstacle[n].yc = OBSYMIN + ((double)j + 0.5)*dy; + if (obstacle[n].xc > OBSXMAX) obstacle[n].xc += (OBSXMAX - OBSXMIN); obstacle[n].radius = OBSTACLE_RADIUS; if ((i+j)%2 == 0) { @@ -7241,13 +7241,14 @@ int add_particle(double x, double y, double vx, double vy, double mass, short in particle[i].xc = x; particle[i].yc = y; - particle[i].radius = MU; + particle[i].radius = MU_ADD; // particle[i].radius = MU*sqrt(mass); particle[i].active = 1; particle[i].neighb = 0; particle[i].diff_neighb = 0; particle[i].thermostat = 1; - particle[i].charge = CHARGE; + particle[i].charge = CHARGE_ADD; + if ((ADD_ALTERNATE_CHARGE)&&(rand()%2 == 1)) particle[i].charge *= -1.0; particle[i].energy = 0.0; particle[i].emean = 0.0; @@ -8357,7 +8358,7 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, } /* draw crosses/bars on charged particles */ - if (TWO_TYPES) +// if (TWO_TYPES) { if ((DRAW_CROSS)&&(particle.charge > 0.0)) { @@ -8468,7 +8469,7 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], /* draw tracer particle trajectory */ { int i, j, time, p, width; - double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum; + double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum, lum1; // blank(); glLineWidth(TRAJECTORY_WIDTH); @@ -8493,7 +8494,12 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], time = traj_length - i; lum = 0.9 - TRACER_LUM_FACTOR*(double)time/(double)(TRAJECTORY_LENGTH*TRACER_STEPS); if (lum < 0.0) lum = 0.0; - glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); + if (FILL_OBSTACLE_TRIANGLES) + { + lum1 = 0.5*(1.0-lum); + glColor3f(lum1 + lum*rgb[0], lum1 + lum*rgb[1], lum1 + lum*rgb[2]); + } + else glColor3f(lum*rgb[0], lum*rgb[1], lum*rgb[2]); if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); @@ -12022,7 +12028,7 @@ int tracer_n[N_TRACER_PARTICLES]) printf("Added a particle at (%.5lg, %.5lg)\n\n\n", x, y); fprintf(lj_log, "Added a particle at (%.5lg, %.5lg)\n\n\n", x, y); - add_particle(x, y, V_INITIAL, 0.05*V_INITIAL*(double)rand()/RAND_MAX, PARTICLE_MASS, type, particle); + add_particle(x, y, V_INITIAL_ADD, 0.05*V_INITIAL_ADD*(double)rand()/RAND_MAX, PARTICLE_ADD_MASS, type, particle); printf("Particle number %i\n", ncircles-1); diff --git a/sub_part_billiard.c b/sub_part_billiard.c index 02a594f..b74cdad 100644 --- a/sub_part_billiard.c +++ b/sub_part_billiard.c @@ -191,6 +191,7 @@ void rgb_color_scheme_minmax(int i, double rgb[3]) /* saturation = r, luminosity = 0.5 */ } + void rgb_color_scheme_density(int part_number, double rgb[3], double minprop) /* color scheme with specified color interval */ { @@ -722,6 +723,13 @@ void draw_billiard() /* draws the billiard boundary */ draw_circle(x0, 0.0, r, NSEG); draw_circle(-x0, 0.0, r, NSEG); } + + if (ABSORBING_CIRCLES) + { + rgb[0] = 0.7; rgb[1] = 0.7; rgb[2] = 0.7; + for (k=0; k DPI) config[1] -= DPI; + + /* experimental */ + if (ABSORBING_CIRCLES) for (i=0; i margin)&&(t > 2.0*MU)) /* there is an intersection with circle i */ + { + config[0] = DUMMY_ABSORBING; + config[1] = PI; + } + } config[2] = 0.0; /* running time */ config[3] = module2(x1-pos[0], y1-pos[1]); /* distance to collision */ @@ -6792,7 +6812,7 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES], int i, j, k, l, n, z, ii, jj, terni[SDEPTH], ternj[SDEPTH], quater[SDEPTH], cond, p, q, block, nblocks, ww; short int vkoch[NMAXCIRCLES], turnright; 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; + x1, y1, x2, y2, dx, dy, padding = 0.02, width = 0.01, xright, yright, xtop, ytop, rand_factor, rmin, rmax, dr, phi, dphi, rscat, omegab, swb; double *maze_coords; t_maze* maze; @@ -6867,14 +6887,17 @@ void init_polyline(t_segment polyline[NMAXPOLY], t_circle circles[NMAXCIRCLES], } case (P_POLYRING): { - nsides = 2*NPOLY; - ncircles = 0; + nsides = NPOLY + NPOLY_B; + if (ABSORBING_CIRCLES) ncircles = NABSCIRCLES; + else ncircles = 0; omega = DPI/(double)NPOLY; + omegab = DPI/(double)NPOLY_B; sw = sin(omega/2.0); + swb = sin(omegab/2.0); for (i=0; i= PI) angle -= DPI; + rde[i*NY+j].theta = angle; + } + else rde[i*NY+j].theta = 0.0; + } +} + double colors_rps(int type) { switch (type) { @@ -4888,7 +4908,7 @@ void compute_kuramoto_int_planar(double phi_in[NX*NY], double phi_out[NX*NY], sh for (j=1; j PI) value -= DPI; + } +// if (vabs(value) > PI) printf("Warning, color out of range\n"); + color_scheme_palette(C_BASIC_LINEAR, palette, value/PI, 1.0, 1, rgb); +// color_scheme_palette(C_ONEDIM_LINEAR, palette, value/PI, 1.0, 1, rgb); +// if (value > PI) value -= PI; +// if (value < 0.0) value += PI, +// amp_to_rgb_palette(value/PI, rgb, palette); + break; + } case (Z_RGB): { /* TO DO */ @@ -5757,6 +5893,18 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, // printf("[compute_rde_fields] zplot = %i\n", zplot); switch (RDE_EQUATION) { + case (E_KURAMOTO): + { + compute_theta_kuramoto(phi, xy_in, rde); + + /* experimental */ + if ((zplot == Z_VORTICITY)||(cplot == Z_VORTICITY)||(zplot == Z_VORTICITY_ABS)||(cplot == Z_VORTICITY_ABS)) + { + compute_gradient_theta(rde); + compute_curl(rde); + } + break; + } case (E_RPS): { if ((COMPUTE_THETA)||(COMPUTE_THETAZ)) @@ -5850,6 +5998,12 @@ void init_zfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, t_ for (i=0; i= wsphere[i*NY+j].altitude + FLOODING_VSHIFT_2D))); - else draw = wsphere[i*NY+j].indomain; +// else draw = wsphere[i*NY+j].indomain; + else draw = 1; - for (p=0; p NX) ii -= NX; @@ -6419,7 +6589,7 @@ void draw_wave_2d_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere ws first = 0; } -// printf("Drawing wave\n"); + printf("Drawing wave\n"); /* draw the field */ glBegin(GL_QUADS); for (i=0; i= wsphere[i*NY+j].altitude + FLOODING_VSHIFT))); - else draw = wsphere[i*NY+j].indomain; -// draw = 1; - for (p=0; p= dpoisson*dpoisson); - /* new circle is in domain */ + + /* take care of periodic boundary conditions */ + if (periodic_bc) for (p=-1; p<2; p++) for (q=-1; q<2; q++) + { + x1 = circles[k].xc + (double)p*(xmax - xmin); + y1 = circles[k].yc + (double)q*(ymax - ymin); + far = far*((x - x1)*(x - x1) + (y - y1)*(y - y1) >= dpoisson*dpoisson); + } +// else +// { + /* new circle is in domain */ +// far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin); +// } far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin); } if (far) /* accept new circle */ @@ -1168,7 +1180,7 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern } case (C_RINGS_POISSONDISC): { - ncircles = generate_poisson_discs(circles, YMIN, YMAX, YMIN, YMAX, PDISC_FACTOR*MU); + ncircles = generate_poisson_discs(circles, YMIN, YMAX, YMIN, YMAX, PDISC_FACTOR*MU, 0); for (i=0; i 0.0) + { + polyline[nlines].x = x1 + MU*cos(alpha1); + polyline[nlines].y = y1 + MU*sin(alpha1); + nlines++; + polyline[nlines].x = x1 + dx - MU*cos(alpha1); + polyline[nlines].y = y1 + MU*sin(alpha1); + nlines++; + + polyline[nlines].x = x1 + MU*cos(alpha1); + polyline[nlines].y = y1 - MU*sin(alpha1); + nlines++; + polyline[nlines].x = x1 + dx - MU*cos(alpha1); + polyline[nlines].y = y1 - MU*sin(alpha1); + nlines++; + + polyrect[nrect].x1 = x1 + MU*cos(alpha1); + polyrect[nrect].y1 = y1 - MU*sin(alpha1); + polyrect[nrect].x2 = x1 + dx - MU*cos(alpha1); + polyrect[nrect].y2 = y1 + MU*sin(alpha1); + nrect++; + } + + if (alpha2 > 0.0) + { + polyline[nlines].x = x1 + MU*sin(alpha2); + polyline[nlines].y = y1 + MU*cos(alpha2); + nlines++; + polyline[nlines].x = x1 + MU*sin(alpha2); + polyline[nlines].y = y1 + dy - MU*cos(alpha2); + nlines++; + + polyline[nlines].x = x1 - MU*sin(alpha2); + polyline[nlines].y = y1 + MU*cos(alpha2); + nlines++; + polyline[nlines].x = x1 - MU*sin(alpha2); + polyline[nlines].y = y1 + dy - MU*cos(alpha2); + nlines++; + + polyrect[nrect].x1 = x1 - MU*sin(alpha2); + polyrect[nrect].y1 = y1 + MU*cos(alpha2); + polyrect[nrect].x2 = x1 + MU*sin(alpha2); + polyrect[nrect].y2 = y1 + dy - MU*cos(alpha2); + nrect++; + } + } + } + + for (i=0; i XMAX - MU) + { + circles[ncirc].xc = circles[i].xc - XMAX + XMIN; + circles[ncirc].yc = circles[i].yc; + ncirc++; + printf("%i circles\n", ncirc); + } + if (circles[i].yc < YMIN + MU) + { + circles[ncirc].xc = circles[i].xc; + circles[ncirc].yc = circles[i].yc + YMAX - YMIN; + ncirc++; + printf("%i circles\n", ncirc); + } + if (circles[i].yc > YMAX - MU) + { + circles[ncirc].xc = circles[i].xc; + circles[ncirc].yc = circles[i].yc - YMAX + YMIN; + ncirc++; + printf("%i circles\n", ncirc); + } + } + + n_neighbours = (int *)malloc(ncirc*sizeof(int)); + neighbour_angles = (double *)malloc(ncirc*maxneigh*sizeof(double)); + + for (i=0; i neighbour_angles[i*maxneigh + k + 1]) + { + temp = neighbour_angles[i*maxneigh + k]; + neighbour_angles[i*maxneigh + k] = neighbour_angles[i*maxneigh + k + 1]; + neighbour_angles[i*maxneigh + k + 1] = temp; + } + } + } + + /* initialize arcs */ + for (j = 0; j < n; j++) + { + polyarc[narcs].xc = circles[i].xc; + polyarc[narcs].yc = circles[i].yc; + polyarc[narcs].r = MU; + polyarc[narcs].angle1 = neighbour_angles[i*maxneigh+j] + alpha; + if (j < n-1) + polyarc[narcs].dangle = neighbour_angles[i*maxneigh+j+1] - neighbour_angles[i*maxneigh+j] - 2.0*alpha; + else + polyarc[narcs].dangle = neighbour_angles[i*maxneigh] - neighbour_angles[i*maxneigh+j] + DPI - 2.0*alpha; + narcs++; + } + } + + 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); + break; + } case (Z_EULER_VORTICITY): { value = min + 1.0*dy*(double)(j); diff --git a/sub_wave_3d_rde.c b/sub_wave_3d_rde.c index 5b7d559..7248239 100644 --- a/sub_wave_3d_rde.c +++ b/sub_wave_3d_rde.c @@ -921,6 +921,7 @@ int ij_to_sphere(int i, int j, double r, t_wave_sphere wsphere[NX*NY], double xy if (use_wave_radius) { newr = wsphere[i*NY+j].radius; + newr += (RSCALE_B-1.0)*r; xyz[0] *= newr; xyz[1] *= newr; xyz[2] *= newr; @@ -958,6 +959,7 @@ int ij_to_sphere_hres(int i, int j, double r, t_wave_sphere wsphere[HRES*HRES*NX if (use_wave_radius) { newr = wsphere[i*HRES*NY+j].radius; + newr += (RSCALE_B-1.0)*r; xyz[0] *= newr; xyz[1] *= newr; xyz[2] *= newr; @@ -2745,6 +2747,12 @@ void draw_circular_color_scheme_palette_3d(double x1, double y1, double radius, value = dy_phase*(double)(j); color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; + } + case (Z_AMPLITUDE): + { + value = min + 1.0*dy*(double)(j); + color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); + break; } case (Z_POLAR): { diff --git a/wave_3d.c b/wave_3d.c index 2ba99ca..dae3dc2 100644 --- a/wave_3d.c +++ b/wave_3d.c @@ -43,7 +43,7 @@ #include #include -#define MOVIE 1 /* set to 1 to generate movie */ +#define MOVIE 0 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ #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 */ @@ -52,10 +52,8 @@ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1150 /* window height */ -// #define NX 1920 /* number of grid points on x axis */ -// #define NY 1150 /* number of grid points on y axis */ -#define NX 3000 /* number of grid points on x axis */ -#define NY 1600 /* number of grid points on y axis */ +#define NX 1920 /* number of grid points on x axis */ +#define NY 1150 /* number of grid points on y axis */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ @@ -68,36 +66,39 @@ /* Choice of the billiard table */ -#define B_DOMAIN 76 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 96 /* 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 */ #define COMPARISON 0 /* set to 1 to compare two different patterns */ #define B_DOMAIN_B 20 /* second domain shape, for comparisons */ #define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */ +#define IMAGE_FILE 5 /* for option D_IMAGE */ -#define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */ +#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ #define IOR 181 /* choice of index of refraction, see list in global_pdes.c */ #define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */ #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 1000 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define PDISC_FACTOR 3.25 /* controls density of Poisson disc process (default: 3.25) */ +#define PDISC_FACTOR 2.3 /* controls density of Poisson disc process (default: 3.25) */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ +#define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ -#define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ -#define MU 0.35 /* parameter controlling the dimensions of domain */ -#define NPOLY 6 /* number of sides of polygon */ +#define LAMBDA 1.25 /* parameter controlling the dimensions of domain */ +#define MU 0.065 /* parameter controlling the dimensions of domain */ +#define NPOLY 4 /* 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 2000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 20.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 30 /* number of grid point for grid of disks */ -#define NGRIDY 18 /* number of grid point for grid of disks */ -#define WALL_WIDTH 0.1 /* width of wall separating lenses */ +#define NGRIDX 15 /* number of grid point for grid of disks */ +#define NGRIDY 10 /* number of grid point for grid of disks */ +#define WALL_WIDTH 0.022 /* width of wall separating lenses */ +#define WALL_WIDTH_RND 0.5 /* proportion of width of width for random arrangements */ #define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define X_SHOOTER -0.2 @@ -145,8 +146,9 @@ /* 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 9 /* period of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 8 /* period of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ +#define MAX_PULSING_TIME 1500 /* max time for adding pulses */ #define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ #define WAVE_PACKET_SOURCE_TYPE 1 /* type of wave packet sources */ @@ -155,14 +157,14 @@ /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 2 +#define B_COND 1 #define PRECOMPUTE_BC 0 /* set to 1 to compute neighbours for Laplacian in advance */ /* Parameters for length and speed of simulation */ -#define NSTEPS 1800 /* number of frames of movie */ -#define NVID 15 /* number of iterations between images displayed on screen */ +#define NSTEPS 3200 /* number of frames of movie */ +#define NVID 3 /* 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 */ @@ -179,10 +181,9 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 1.0 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.00001 /* variance of initial condition */ -#define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ - +#define INITIAL_AMP 0.75 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */ +#define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ @@ -196,7 +197,7 @@ #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 NON_DIRICHLET_BC 1 /* 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 0 /* set to 1 to draw boundary */ #define DRAW_BILLIARD_FRONT 0 /* set to 1 to draw front of boundary after drawing wave */ @@ -214,29 +215,29 @@ #define REP_AXO_3D 0 /* linear projection (axonometry) */ #define REP_PROJ_3D 1 /* projection on plane orthogonal to observer line of sight */ -#define ROTATE_VIEW 0 /* set to 1 to rotate position of observer */ +#define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ /* Color schemes */ -#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 10 /* 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 VSCALE_AMPLITUDE 2.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define SLOPE 0.75 /* sensitivity of color on wave amplitude */ +#define VSCALE_AMPLITUDE 0.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ #define VSHIFT_AMPLITUDE 0.0 /* additional shift for wave amplitude */ -#define VSCALE_ENERGY 3.0 /* additional scaling factor for color scheme P_3D_ENERGY */ +#define VSCALE_ENERGY 1.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 LOG_SCALE 0.75 /* scaling factor for energy log representation */ -#define LOG_SHIFT 0.5 /* shift of colors on log scale */ +#define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ +#define E_SCALE 30.0 /* scaling factor for energy representation */ +#define LOG_SCALE 1.5 /* scaling factor for energy log representation */ +#define LOG_SHIFT 0.25 /* 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 4000.0 /* scaling factor for energy flux representation */ @@ -266,13 +267,13 @@ #define MAZE_WIDTH 0.02 /* half width of maze walls */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 3.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 2.0 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 0.6 /* 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 ADD_POTENTIAL 1 /* set to 1 to add potential to z coordinate */ +#define ADD_POTENTIAL 0 /* set to 1 to add potential to z coordinate */ #define POTENTIAL 10 #define POT_FACT 200.0 #define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ @@ -289,7 +290,7 @@ #define HRES 1 /* dummy, only used by rde.c */ #define SHADE_2D 0 /* set to 1 to add pseudo-3d shading effect */ #define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */ -#define N_SOURCES 1 /* number of sources, for option draw_sources */ +#define N_SOURCES 2 /* number of sources, for option draw_sources */ #define XYIN_INITIALISED (B_DOMAIN == D_IMAGE) /* end of constants only used by sub_wave and sub_maze */ @@ -304,14 +305,14 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; double light[3] = {0.816496581, -0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {8.0, 6.0, 6.0}; /* location of observer for REP_PROJ_3D representation */ +double observer[3] = {-8.0, -6.0, 6.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.15 /* 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 Z_SCALING_FACTOR 0.05 /* overall scaling factor of z axis for REP_PROJ_3D representation */ +#define XY_SCALING_FACTOR 1.8 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ #define XSHIFT_3D -0.1 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D 0.3 /* overall y shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.2 /* overall y shift for REP_PROJ_3D representation */ #include "global_pdes.c" /* constants and global variables */ @@ -950,11 +951,11 @@ 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, angle = 0.0, lambda1, y, x1, sign1, omega, phase_shift; + double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle = 0.0, lambda1, y, x1, sign1, omega, phase_shift, source_amp[N_SOURCES]; 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, source_counter = 0, k, p, q; + int i, j, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, k, p, q, source, source_periods[N_SOURCES]; static int counter = 0, first_source = 1; long int wave_value; t_wave *wave; @@ -1003,15 +1004,26 @@ void animation() } printf("Polygons initialized\n"); + if (XYIN_INITIALISED) init_xyin_from_image_3d(xy_in); + /* initialise polyline for von Koch and similar domains */ - npolyline = init_polyline(MDEPTH, polyline); - for (i=0; i #define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ #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 */ @@ -55,17 +55,13 @@ /* General geometrical parameters */ -// #define WINWIDTH 1920 /* window width */ -#define WINWIDTH 1150 /* window width */ +#define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1150 /* window height */ -// #define NX 3840 /* number of grid points on x axis */ -#define NX 2300 /* number of grid points on x axis */ +#define NX 3840 /* number of grid points on x axis */ #define NY 2300 /* number of grid points on y axis */ -// #define XMIN -2.0 -// #define XMAX 2.0 /* x interval */ -#define XMIN -1.197916667 -#define XMAX 1.197916667 /* x interval */ +#define XMIN -2.0 +#define XMAX 2.0 /* x interval */ #define YMIN -1.197916667 #define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ @@ -76,9 +72,9 @@ /* Choice of the billiard table */ -#define B_DOMAIN 92 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 96 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */ +#define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */ #define IMAGE_FILE 5 /* for option D_IMAGE */ #define COMPARISON 0 /* set to 1 to compare two different patterns (beta) */ @@ -87,22 +83,24 @@ #define P_PERCOL 0.15 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define PDISC_FACTOR 3.5 /* controls density of Poisson disc process (default: 3.25) */ +#define PDISC_FACTOR 2.3 /* controls density of Poisson disc process (default: 3.25) */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ +#define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ -#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ -#define MU 0.75 /* parameter controlling the dimensions of domain */ +#define LAMBDA 1.25 /* parameter controlling the dimensions of domain */ +#define MU 0.065 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */ -#define NPOLY 7 /* number of sides of polygon */ +#define NPOLY 6 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 6 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #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 60 /* number of grid point for grid of disks */ -#define NGRIDY 25 /* number of grid point for grid of disks */ -#define WALL_WIDTH 0.05 /* width of wall separating lenses */ +#define NGRIDX 15 /* number of grid point for grid of disks */ +#define NGRIDY 10 /* number of grid point for grid of disks */ +#define WALL_WIDTH 0.017 /* width of wall separating lenses */ +#define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */ #define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define X_SHOOTER -0.2 @@ -110,11 +108,11 @@ #define X_TARGET 0.4 #define Y_TARGET 0.7 /* shooter and target positions in laser fight */ -#define ISO_XSHIFT_LEFT -2.9 +#define ISO_XSHIFT_LEFT -2.9 #define ISO_XSHIFT_RIGHT 1.4 -#define ISO_YSHIFT_LEFT -0.15 -#define ISO_YSHIFT_RIGHT -0.15 -#define ISO_SCALE 0.5 /* coordinates for isospectral billiards */ +#define ISO_YSHIFT_LEFT -0.2 +#define ISO_YSHIFT_RIGHT 0.15 +#define ISO_SCALE 0.475 /* coordinates for isospectral billiards */ /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard below */ @@ -132,7 +130,7 @@ #define AMPLITUDE 0.5 /* amplitude of periodic excitation */ #define ACHIRP 0.25 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.08 /* Courant number */ +#define COURANT 0.1 /* Courant number */ #define COURANTB 0.025 /* Courant number in medium B */ #define GAMMA 0.0 /* damping factor in wave equation */ #define GAMMAB 0.0 /* damping factor in wave equation */ @@ -148,11 +146,11 @@ /* 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 49 /* period of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 8 /* period of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ -#define N_SOURCES 7 /* number of sources, for option draw_sources */ +#define N_SOURCES 2 /* number of sources, for option draw_sources */ #define ALTERNATE_SOURCE_PHASES 0 /* set to 1 to alternate initial phases of sources */ -#define MAX_PULSING_TIME 100 /* max time for adding pulses */ +#define MAX_PULSING_TIME 1500 /* max time for adding pulses */ #define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ #define WAVE_PACKET_SOURCE_TYPE 3 /* type of wave packet sources */ @@ -163,12 +161,12 @@ /* Boundary conditions, see list in global_pdes.c */ -#define B_COND 2 +#define B_COND 1 /* Parameters for length and speed of simulation */ -#define NSTEPS 1350 /* number of frames of movie */ -#define NVID 8 /* number of iterations between images displayed on screen */ +#define NSTEPS 3200 /* number of frames of movie */ +#define NVID 7 /* 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 */ @@ -185,20 +183,20 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 1.4 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.00005 /* variance of initial condition */ +#define INITIAL_AMP 0.75 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.0005 /* variance of initial condition */ #define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -#define PLOT 8 +#define PLOT 0 -#define PLOT_B 0 /* plot type for second movie */ +#define PLOT_B 8 /* plot type for second movie */ /* Color schemes */ -#define COLOR_PALETTE 12 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 10 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -209,13 +207,13 @@ #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 VSHIFT_AMPLITUDE -0.5 /* additional shift for wave amplitude */ +#define VSHIFT_AMPLITUDE 0.0 /* additional shift for wave amplitude */ #define VSCALE_AMPLITUDE 0.2 /* additional scaling factor for wave amplitude */ -#define E_SCALE 13.5 /* scaling factor for energy representation */ +#define E_SCALE 30.0 /* scaling factor for energy representation */ #define LOG_SCALE 0.75 /* scaling factor for energy log representation */ #define LOG_SHIFT 0.75 /* shift of colors on log scale */ #define FLUX_SCALE 250.0 /* scaling factor for energy flux represtnation */ -#define AVRG_E_FACTOR 0.9 /* controls time window size in P_AVERAGE_ENERGY scheme */ +#define AVRG_E_FACTOR 0.85 /* controls time window size in P_AVERAGE_ENERGY scheme */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */ #define SHADE_2D 1 /* set to 1 to add pseudo-3d shading effect */ @@ -230,8 +228,8 @@ #define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ #define COLORBAR_RANGE 1.5 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 0.3 /* scale of color scheme bar for 2nd part */ -#define ROTATE_COLOR_SCHEME 1 /* set to 1 to draw color scheme horizontally */ +#define COLORBAR_RANGE_B 0.12 /* scale of color scheme bar for 2nd part */ +#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */ #define CIRC_COLORBAR_B 0 /* set to 1 to draw circular color scheme */ @@ -588,7 +586,7 @@ void draw_color_bar_palette(int plot, double range, int palette, int circular, i // double width = 0.2; 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(-1.0, -0.85, XMAX - 0.1, -1.0, plot, -range, range, palette, fade, fade_value); else if (circular) draw_circular_color_scheme_palette_fade(XMAX - 2.0*width, YMIN + 2.0*width, 1.5*width, plot, -range, range, palette, fade, fade_value); else @@ -597,10 +595,10 @@ void draw_color_bar_palette(int plot, double range, int palette, int circular, i void animation() { - double time, scale, ratio, startleft[2], startright[2], sign[N_SOURCES], r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle = 0.0, x1, ior_angle = 0.0, omega, phase_shift, vshift, dsource, finv, source_amp, nx, ny, r; + double time, scale, ratio, startleft[2], startright[2], sign[N_SOURCES], r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle = 0.0, x1, ior_angle = 0.0, omega, phase_shift, vshift, dsource, finv, source_amp[N_SOURCES], nx, ny, r; double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *average_energy[NX], *color_scale[NX], *total_flux, *tcc_table[NX], *tgamma_table[NX], *fade_table; short int *xy_in[NX]; - int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q, first_source = 1, imin, imax, ij[2], source, source_period, source_shift[N_SOURCES]; + int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q, first_source = 1, imin, imax, ij[2], source, source_period, source_shift[N_SOURCES], source_periods[N_SOURCES]; // static int image_counter = 0; int image_counter = 0; long int wave_value; @@ -646,16 +644,29 @@ void animation() printf("Polygons initialized\n"); /* initialise polyline for von Koch and similar domains */ - npolyline = init_polyline(MDEPTH, polyline); - for (i=0; i