diff --git a/global_3d.c b/global_3d.c index 7237b39..5880a2b 100644 --- a/global_3d.c +++ b/global_3d.c @@ -12,6 +12,7 @@ #define E_EULER_INCOMP 6 /* incompressible Euler equation */ #define E_EULER_COMP 7 /* compressible Euler equation */ #define E_SHALLOW_WATER 8 /* shallow water equation */ +#define E_KURAMOTO 9 /* Kuramoto-type equation (nearest neighbor coupling) */ /* Choice of potential */ diff --git a/global_ljones.c b/global_ljones.c index 36be704..77c522b 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -13,15 +13,17 @@ #define NMAXCIRCLES 200000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define MAXNEIGH 20 /* max number of neighbours kept in memory */ -#define NMAXOBSTACLES 1000 /* max number of obstacles */ -#define NMAXSEGMENTS 1000 /* max number of repelling segments */ -#define NMAXGROUPS 50 /* max number of groups of segments */ +#define NMAXOBSTACLES 5000 /* max number of obstacles */ +#define NMAXSEGMENTS 1200 /* max number of repelling segments */ +#define NMAXGROUPS 100 /* max number of groups of segments */ #define NMAXCOLLISIONS 200000 /* max number of collisions */ #define NMAXPARTNERS 30 /* max number of partners in molecule */ #define NMAXPARTNERMOLECULES 30 /* max number of partners of a molecule */ #define NMAXPARTINCLUSTER 500 /* max number of particles in cluster */ #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 C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -46,6 +48,7 @@ /* pattern of additional obstacles */ #define O_CORNERS 0 /* obstacles in the corners (for Boy b.c.) */ #define O_GALTON_BOARD 1 /* Galton board pattern */ +#define O_GALTON_BOARD_WIDE 11 /* wide 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 */ @@ -58,6 +61,22 @@ #define O_SIEVE_B 81 /* obstacles form a sieve, v2 */ #define O_SIEVE_LONG 82 /* obstacles form a long sieve */ #define O_SIEVE_LONG_B 83 /* obstacles form a long sieve, version with varying spacing */ +#define O_SIEVE_VIDEO1400 84 /* sieve for vide #1400 */ +#define O_SQUARE 9 /* square lattice */ +#define O_SQUARE_TWOMASSES 91 /* square lattice with alternating masses */ +#define O_POISSON_DISC 20 /* Poisson disc sample */ + +/* type of obstacle pinning */ +#define OP_NPINNED 0 /* obstacles with NMAX_OBSTACLE_PINNED neighbours or less are pinned*/ +#define OP_LEFT 1 /* obstacles on the left are pinned */ +#define OP_LEFTTOPBOT 2 /* obstacles on the left, top or bottom are pinned */ +#define OP_CORNERS 3 /* obstacles in corners are pinned */ +#define OP_LEFTCORNERS 4 /* obstacles in left half and corners are pinned */ +#define OP_BNDRY_STEP 5 /* obstacles on boundary are pinned when at distance BDRY_PINNING_STEP */ + +/* type of obstacle recoupling */ +#define OR_FIXED_DIST 0 /* resample within fixed distance */ +#define OR_BREAK_MAX 1 /* break bonds when max length is reached */ /* pattern of additional repelling segments */ #define S_RECTANGLE 0 /* segments forming a rectangle */ @@ -95,6 +114,7 @@ #define S_COANDA_SHORT 26 /* shorter wall for Coanda effect */ #define S_CYLINDER 27 /* walls at top and bottom, for cylindrical b.c. */ #define S_TREE 28 /* Christmas tree(s) */ +#define S_TREES_B 281 /* Christmas trees, variant 2 */ #define S_CONE 29 /* cone */ #define S_CONVEYOR_BELT 30 /* conveyor belt */ #define S_TWO_CONVEYOR_BELTS 31 /* two angled conveyor belts */ @@ -105,8 +125,11 @@ #define S_CONVEYOR_SIEVE 35 /* conveyor belts for polygon sieve */ #define S_CONVEYOR_SIEVE_B 351 /* conveyor belts for polygon sieve, v2 with backward top conveyor */ #define S_CONVEYOR_SIEVE_LONG 352 /* conveyor belts for long polygon sieve */ +#define S_CONVEYORS_1400 353 /* conveyors for video #1400 */ #define S_MASS_SPECTROMETER 36 /* bins for mass spectrometer */ #define S_WIND_FORCE 361 /* bins for sorting by wind force */ +#define S_BINS_GALTON 362 /* bins for Galton board */ +#define S_BINS_GALTON_WIDE 363 /* more bins for Galton board */ /* particle interaction */ @@ -151,6 +174,7 @@ #define BC_ABSORBING 20 /* "no-return" boundary conditions outside BC area */ #define BC_REFLECT_ABS 21 /* reflecting on lower boundary, and "no-return" boundary conditions outside BC area */ #define BC_REFLECT_ABS_BOTTOM 22 /* absorbing on lower boundary, and reflecting elsewhere */ +#define BC_REFLECT_ABS_RIGHT 23 /* absorbing on right boundary, and reflecting elsewhere */ /* Regions for partial thermostat couplings */ @@ -302,6 +326,12 @@ #define IP_X 0 /* color depends on x coordinate of initial position */ #define IP_Y 1 /* color depends on y coordinate of initial position */ +/* Space dependence of electric field */ + +#define EF_CONST 0 /* constant magnetic field */ +#define EF_SQUARE 1 /* magnetic field concentrated in square */ +#define EF_LEFT 2 /* magnetic field concentrated region x < YMIN */ + /* Space dependence of magnetic field */ #define BF_CONST 0 /* constant magnetic field */ @@ -336,6 +366,16 @@ #define BG_CHARGE 2 /* background color depends on charge density */ #define BG_EKIN 3 /* background color depends on kinetic energy */ #define BG_FORCE 4 /* background color depends on total force */ +#define BG_EOBSTACLES 5 /* background color depends on obstacle energy */ +#define BG_EKIN_OBSTACLES 6 /* background color depends on kinetic energy plus obstacle energy */ +#define BG_DIR_OBSTACLES 7 /* background color depends on direction of velocity of obstacles */ +#define BG_POS_OBSTACLES 8 /* background color depends on displacement of obstacles */ + +/* Obstacle color schemes */ + +#define OC_ENERGY 0 /* obstacle color depends on energy */ +#define OC_DIRECTION 1 /* obstacle color depends on direction */ +#define OC_DIRECTION_ENERGY 2 /* obstacle hue depends on direction, luminosity on energy */ /* Particle add regions */ @@ -447,6 +487,8 @@ typedef struct { int number; /* total number of particles in cell */ int particles[HASHMAX]; /* numbers of particles in cell */ + int nobs; /* number of obstacles in cell */ + int obstacle; /* obstacle number */ int nneighb; /* number of neighbouring cells */ int neighbour[9]; /* numbers of neighbouring cells */ double x1, y1, x2, y2; /* coordinates of hashcell corners */ @@ -475,16 +517,47 @@ typedef struct typedef struct { double xc, yc, radius; /* center and radius of circle */ + double vx, vy; /* speed for option OSCILLATING_OBSTACLES */ + double fx, fy; /* force for option OSCILLATING_OBSTACLES */ double xc0, yc0; /* center of oscillation for option RATTLE_OBSTACLES */ + double energy; /* total energy (kinetic + potential) of obstacle */ short int active; /* circle is active */ double charge; /* charge of obstacle, for EM simulations */ + double mass; /* mass for option OSCILLATING_OBSTACLES */ double omega0, omega; /* speed of rotation */ double angle; /* angle of obstacle */ short int oscillate; /* has value 1 if the obstacles oscillates over time */ int period; /* oscillation period */ double amplitude, phase; /* amplitude and phase of oscillation */ + int nneighb; /* number of neighbours, for option OSCILLATING_OBSTACLES */ + int neighb[NMAX_OBSTACLE_NEIGHBOURS]; /* list of neighour numbers */ + double eqdist[NMAX_OBSTACLE_NEIGHBOURS]; /* equilibrium distance to neighbours */ + short int pinned; /* has value 1 if particle is pinned to a fixed position */ + short int chessboard; /* has value 1 on a chessboard, for some arrangements */ } t_obstacle; +typedef struct +{ + int i[3]; /* indices of obstacles forming triangle */ + int facet; /* index of facet triangle belongs to */ + short int acute; /* value depends on central angle */ + short int infacet; /* has value 1 if triangle belongs to a facet */ + double area0; /* initial area */ +} t_otriangle; + +typedef struct +{ + int ntriangles; /* number of triangles belonging to facet */ + int triangle[NMAX_TRIANGLES_PER_FACET]; /* indices of triangles belonging to facet */ + double area0; /* initial area */ +} t_ofacet; + + +typedef struct +{ + double xc, yc; /* point in R^2, for Poisson disc generation */ +} t_point; + typedef struct { double x1, x2, y1, y2; /* extremities of segment */ @@ -586,6 +659,6 @@ typedef struct -int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0; +int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0, n_otriangles = 0, n_ofacets = 0; FILE *lj_log; diff --git a/global_pdes.c b/global_pdes.c index 3abed13..e49f323 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -109,6 +109,12 @@ #define D_TWOCIRCLES 79 /* two circles of different size */ #define D_POLYCIRCLES 791 /* one large circle and NPOLY small ones */ #define D_POLYCIRCLES_ANGLED 792 /* variant of D_POLYCIRCLES with angled small cavities */ +#define D_TWOCIRCLES_ELLIPSE 793 /* ellipse with two adjacent circular cavities */ +#define D_CIRCLES_ELLIPSE 794 /* ellipse with NPOLY adjacent circles */ + +#define D_CARDIOID 90 /* cardioid */ +#define D_NEPHROID 91 /* nephroid */ +#define D_EPICYCLOID 92 /* epicycloid */ /* for wave_sphere.c */ @@ -126,6 +132,16 @@ #define D_SPHERE_MAZE_SPIRAL 101 /* circular maze on the sphere with slanted walls */ #define D_SPHERE_MAZE_WAVE 102 /* circular maze on the sphere with wavy walls */ +/* image files for D_IMAGE option */ + +#define IM_PHOTONS 0 /* word "photons" */ +#define IM_GUITAR 1 /* guitar shape */ +#define IM_GUITAR_NOHOLE 11 /* guitar shape without sound hole or neck */ +#define IM_EGUITAR 2 /* electric guitar shape */ +#define IM_HAPPY_NEW_YEAR 3 /* "Happy New Year 2025" */ +#define IM_DICKSON 4 /* Dickson fjord and others nearbys in Greenland */ +#define IM_DICKSON_ZOOM 5 /* zoom on Dickson fjord in Greenland */ + #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 NMAXSOURCES 30 /* maximal number of sources */ @@ -466,3 +482,6 @@ double julia_x = 0.37468, julia_y = 0.21115; /* parameters for Julia sets */ double wave_source_x[NMAXSOURCES], wave_source_y[NMAXSOURCES]; /* position of wave source */ double focus_x = XMAX; /* position of focus */ double michelson_position = 0.0; /* position of mirror in Michelson interferometer */ + +#define NEPICYCLOID 2000 /* number of sampled values for epicycloid */ +double epicycloid_phi_to_r[NEPICYCLOID]; /* angles for epicycloid */ diff --git a/heat.c b/heat.c index db20b5d..9b3db9d 100644 --- a/heat.c +++ b/heat.c @@ -67,6 +67,7 @@ #define B_DOMAIN 27 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 0 /* pattern of circles, see list in global_pdes.c */ +#define IMAGE_FILE 5 /* for option D_IMAGE */ #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 */ diff --git a/lennardjones.c b/lennardjones.c index 92792bb..e9fc290 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -58,42 +58,57 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -2.1 -#define INITXMAX -2.0 /* x interval for initial condition */ -#define INITYMIN -0.7 -#define INITYMAX -0.6 /* y interval for initial condition */ +#define INITXMIN -1.0 +#define INITXMAX -1.9 /* x interval for initial condition */ +#define INITYMIN -0.1 +#define INITYMAX 0.1 /* 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.1 -#define ADDXMAX -2.0 /* x interval for adding particles */ -#define ADDYMIN -0.6 -#define ADDYMAX -0.5 /* y interval for adding particles */ +#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 ADDRMIN 4.75 #define ADDRMAX 6.0 /* r interval for adding particles */ -#define BCXMIN -2.2 -#define BCXMAX 2.0 /* x interval for boundary condition */ +#define BCXMIN -2.5 +#define BCXMAX 2.5 /* x interval for boundary condition */ #define BCYMIN -1.125 -#define BCYMAX 1.4 /* y interval for boundary condition */ +#define BCYMAX 1.125 /* y interval for boundary condition */ -#define OBSXMIN -2.0 +#define OBSXMIN -1.8 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ +#define OBSYMIN -1.125 +#define OBSYMAX 1.125 /* x interval for motion of obstacle */ #define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_ljones.c */ #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 0 /* set to 1 do add fixed circular obstacles */ -#define OBSTACLE_PATTERN 82 /* pattern of obstacles, see list in global_ljones.c */ -#define RATTLE_OBSTACLES 1 /* set to 1 to rattle obstacles (for pattern O_SIEVE_B) */ +#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 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 NMAX_OBSTACLE_NEIGHBOURS 8 /* max number of obstacle neighbours */ +#define NMAX_OBSTACLE_PINNED 10 /* 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 */ +#define OBSTACLE_RECOUPLE_TYPE 1 /* algorithm for recoupling, see OR_ in global_ljones */ +#define OBSTACLE_RECOUPLE_TIME 200 /* time between obstacle recouplings */ +#define UNCOUPLE_MAXLENGTH 2.0 /* length at which bonds decouple */ +#define COUPLE_MINLENGTH 0.5 /* length at which bonds decouple */ -#define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */ -#define SEGMENT_PATTERN 361 /* 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 15 /* pattern of repelling segments, see list in global_ljones.c */ #define ROCKET_SHAPE 3 /* shape of rocket combustion chamber, see list in global_ljones.c */ #define ROCKET_SHAPE_B 3 /* shape of second rocket */ #define NOZZLE_SHAPE 6 /* shape of nozzle, see list in global_ljones.c */ @@ -124,7 +139,7 @@ #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.02 /* parameter controlling radius of particles */ +#define MU 0.007 /* parameter controlling radius of particles */ #define MU_B 0.03 /* parameter controlling radius of particles of second type */ #define NPOLY 3 /* number of sides of polygon */ #define APOLY 0.075 /* angle by which to turn polygon, in units of Pi/2 */ @@ -140,8 +155,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 10 -#define NOBSY 5 /* obstacles for O_HEX obstacle pattern */ +#define NOBSX 20 +#define NOBSY 10 /* obstacles for O_HEX obstacle pattern */ #define NTREES 15 /* number of trees in S_TREES */ #define X_SHOOTER -0.2 @@ -150,9 +165,9 @@ #define Y_TARGET 0.7 /* shooter and target positions in laser fight */ /* Parameters for length and speed of simulation */ - -#define NSTEPS 4800 /* number of frames of movie */ -#define NVID 120 /* number of iterations between images displayed on screen */ + +#define NSTEPS 4500 /* number of frames of movie */ +#define NVID 4500 /* 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 OBSTACLE_INITIAL_TIME 0 /* time after which to start moving obstacle */ @@ -169,29 +184,34 @@ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 1 +#define BOUNDARY_COND 23 /* Plot type, see list in global_ljones.c */ -#define PLOT 23 -#define PLOT_B 13 /* plot type for second movie */ +#define PLOT 13 +#define PLOT_B 12 /* plot type for second movie */ /* Background color depending on particle properties */ #define COLOR_BACKGROUND 0 /* set to 1 to color background */ -#define BG_COLOR 2 /* type of background coloring, see list in global_ljones.c */ -#define BG_COLOR_B 0 /* type of background coloring, see list in global_ljones.c */ +#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 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 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 */ -#define N_PARTICLE_COLORS 300 /* number of colors for P_NUMBER color scheme */ +#define N_PARTICLE_COLORS 5 /* number of colors for P_NUMBER color scheme */ #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 */ /* Color schemes */ @@ -242,8 +262,12 @@ #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_EMIN 10.0 /* energy of particle with coolest color */ -#define PARTICLE_EMAX 2000.0 /* energy of particle with hottest color */ +#define PARTICLE_EMIN 100.0 /* energy of particle with coolest color */ +#define PARTICLE_EMAX 3000.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_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 */ #define HUE_TYPE2 320.0 /* hue of particles of type 2 */ @@ -254,26 +278,26 @@ #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 RANDOM_RADIUS 1 /* set to 1 for random particle radius */ -#define RANDOM_RADIUS_MIN 0.25 /* min of random particle radius (default 0.75) */ -#define RANDOM_RADIUS_RANGE 1.5 /* range of random particle radius (default 0.5) */ -#define ADAPT_MASS_TO_RADIUS 0 /* set to positive value to for mass prop to power of radius */ +#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) */ +#define RANDOM_RADIUS_RANGE 1.0 /* range of random particle radius (default 0.5) */ +#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 3.0e-6 /* time step for particle displacement */ +#define DT_PARTICLE 2.0e-7 /* 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_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) */ -#define DAMPING 150.0 /* damping coefficient of particles */ +#define DAMPING 0.0 /* damping coefficient of particles */ #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_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 800.0 /* initial velocity range */ +#define V_INITIAL 50.0 /* initial velocity range */ #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 */ @@ -289,7 +313,7 @@ #define KSPRING_BOUNDARY 2.0e11 /* confining harmonic potential outside simulation region */ #define KSPRING_OBSTACLE 1.0e9 /* harmonic potential of obstacles */ #define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ -#define GRAVITY 5000.0 /* gravity acting on all particles */ +#define GRAVITY 0.0 /* gravity acting on all particles */ #define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */ #define CIRCULAR_GRAVITY 0 /* set to 1 to have gravity directed to center */ #define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */ @@ -300,10 +324,11 @@ #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 0 /* set to 1 to add an electric field */ -#define EFIELD 30000.0 /* value of electric field */ -#define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ -#define BFIELD 225.0 /* value of magnetic field */ +#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 CHARGE_B 1.0 /* charge of particles of second type */ #define INCREASE_E 0 /* set to 1 to increase electric field */ @@ -312,14 +337,20 @@ #define BFIELD_FACTOR 1.0 /* factor by which to increase magnetic field */ #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 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 */ -#define BFIELD_REGION 1 /* space-dependent magnetic field (0 for constant) */ +#define EFIELD_REGION 0 /* space-dependent magnetic field (0 for constant) */ +#define BFIELD_REGION 0 /* space-dependent magnetic field (0 for constant) */ -#define ADD_WIND 1 /* set to 1 to add a "wind" friction force */ +#define ADD_WIND 0 /* set to 1 to add a "wind" friction force */ #define WIND_FORCE 1.35e6 /* force of wind */ #define WIND_YMIN -0.6 /* min altitude of region with wind */ -#define ROTATION 1 /* set to 1 to include rotation of particles */ +#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 1.0e5 /* force constant in angular dynamics */ @@ -355,7 +386,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.018 /* radius of obstacle for circle boundary conditions */ +#define OBSTACLE_RADIUS 0.03 /* 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 */ @@ -382,17 +413,22 @@ #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 15 /* time interval between adding further particles */ +#define ADD_PERIOD 150 /* 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 SAFETY_FACTOR 4.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ #define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */ -#define N_TRACER_PARTICLES 500 /* number of tracer particles */ -#define TRAJECTORY_LENGTH 4800 /* length of recorded trajectory */ -#define TRACER_LUM_FACTOR 5.0 /* controls luminosity decrease of trajectories with time */ +#define N_TRACER_PARTICLES 100 /* 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 TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */ -#define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */ +#define TRAJECTORY_WIDTH 2 /* width of tracer particle trajectory */ + +#define TRACK_PARTICLE 0 /* set to 1 to track a given particle */ +#define TRACKED_PARTICLE 2 /* number of tracked particle */ +#define TRACK_INITIAL_TIME 900 /* time when starting to track */ #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) */ @@ -416,9 +452,11 @@ #define SHOW_SEGMENTS_PRESSURE 0 /* set to 1 to show (averaged) pressure acting on segments */ #define SEGMENT_PMAX 7.5e7 /* pressure of segment with hottest color */ #define P_AVRG_FACTOR 0.02 /* factor in computation of mean pressure */ +#define INACTIVATE_SEGMENTS_UNDER_PRESSURE 0 /* set to 1 to inactivate segment groups when limit pressure is reached */ +#define SEGMENT_P_INACTIVATE 6.0e7 /* pressure at which to inactivate group */ #define MOVE_SEGMENT_GROUPS 0 /* set to 1 to group segments into moving units */ -#define SEGMENT_GROUP_MASS 500.0 /* mass of segment group */ +#define SEGMENT_GROUP_MASS 500.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 0 /* set to 1 for groups of segments to repel each other */ @@ -478,7 +516,7 @@ #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.015 /* width of wall for BC_RECTANGLE_WALL b.c. */ +#define WALL_WIDTH 0.025 /* 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 */ @@ -517,19 +555,19 @@ #define PARTICLE_MASS_D 2.0 /* mass or partner particle */ #define CHARGE_D 1.0 /* charge of partner particle */ -#define NXMAZE 12 /* width of maze */ -#define NYMAZE 12 /* height of maze */ +#define NXMAZE 16 /* width of maze */ +#define NYMAZE 16 /* height of maze */ #define MAZE_MAX_NGBH 4 /* max number of neighbours of maze cell */ -#define RAND_SHIFT 4 /* seed of random number generator */ -#define MAZE_XSHIFT 0.5 /* horizontal shift of maze */ -#define MAZE_WIDTH 0.01 /* width of maze walls */ +#define RAND_SHIFT 5 /* seed of random number generator */ +#define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ +#define MAZE_WIDTH 0.015 /* width of maze walls */ #define FLOOR_FORCE 1 /* set to 1 to limit force on particle to FMAX */ #define FMAX 1.0e9 /* maximal force */ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 40 /* size of hashgrid in x direction */ +#define HASHX 29 /* 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 */ @@ -549,10 +587,11 @@ #define COUNT_PARTNER_TYPE ((RD_REACTION == CHEM_H2O_H_OH)||(RD_REACTION == CHEM_2H2O_H3O_OH)) #define PAIR_FORCE ((PAIR_PARTICLES)||((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_AGGREGATION)||(RD_REACTION == CHEM_AGGREGATION_CHARGE)||(RD_REACTION == CHEM_AGGREGATION_NNEIGH)||(RD_REACTION == CHEM_POLYGON_AGGREGATION)))) #define COMPUTE_PAIR_TORQUE (KTORQUE_PAIR_ANGLE != 0.0) -#define ADD_CONVEYOR_FORCE ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_BELT)||(SEGMENT_PATTERN == S_TWO_CONVEYOR_BELTS)||(SEGMENT_PATTERN == S_PERIODIC_CONVEYORS)||(SEGMENT_PATTERN == S_TEST_CONVEYORS)||(SEGMENT_PATTERN == S_CONVEYOR_SHOVELS)||(SEGMENT_PATTERN == S_CONVEYOR_MIXED)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_B)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_LONG))) -#define MOVE_CONVEYOR_BELT ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_SHOVELS)||(SEGMENT_PATTERN == S_CONVEYOR_MIXED)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_LONG))) -#define ROTATE_OBSTACLES ((ADD_FIXED_OBSTACLES)&&((OBSTACLE_PATTERN == O_SIEVE)||(OBSTACLE_PATTERN == O_SIEVE_B)||(OBSTACLE_PATTERN == O_SIEVE_LONG))) +#define ADD_CONVEYOR_FORCE ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_BELT)||(SEGMENT_PATTERN == S_TWO_CONVEYOR_BELTS)||(SEGMENT_PATTERN == S_PERIODIC_CONVEYORS)||(SEGMENT_PATTERN == S_TEST_CONVEYORS)||(SEGMENT_PATTERN == S_CONVEYOR_SHOVELS)||(SEGMENT_PATTERN == S_CONVEYOR_MIXED)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_B)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_LONG)||(SEGMENT_PATTERN == S_CONVEYORS_1400))) +#define MOVE_CONVEYOR_BELT ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_SHOVELS)||(SEGMENT_PATTERN == S_CONVEYOR_MIXED)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_LONG)||(SEGMENT_PATTERN == S_CONVEYORS_1400))) +#define ROTATE_OBSTACLES ((ADD_FIXED_OBSTACLES)&&((OBSTACLE_PATTERN == O_SIEVE)||(OBSTACLE_PATTERN == O_SIEVE_B)||(OBSTACLE_PATTERN == O_SIEVE_LONG)||(OBSTACLE_PATTERN == O_SIEVE_VIDEO1400))) #define POLYGON_INTERACTION ((INTERACTION == I_POLYGON)||(INTERACTION == I_POLYGON_ALIGN)) +#define TRACK ((TRACK_SEGMENT_GROUPS)||(TRACK_PARTICLE)) double xshift = 0.0; /* x shift of shown window */ double xspeed = 0.0; /* x speed of obstacle */ @@ -981,7 +1020,8 @@ double radius_schedule(int i) return(1.0 + (MU_RATIO - 1.0)*(double)i/(double)NSTEPS); } -double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], +double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], + 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 *ncoupled, int initial_phase) @@ -993,11 +1033,12 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA if (initial_phase) damping = INITIAL_DAMPING; else damping = DAMPING; -// printf("Evolving particles\n"); +// printf("Evolving %i particles\n", ncircles); - #pragma omp parallel for private(j,xi,totalenergy,a,move) +// #pragma omp parallel for private(j,xi,totalenergy,a,move) for (j=0; j 4.0*OBSTACLE_COUPLING_DIST) dist = 4.0*OBSTACLE_COUPLING_DIST; + if (dist > 5.0*obstacle[i].eqdist[j]) dist = 5.0*obstacle[i].eqdist[j]; + if (dist < 0.2*OBSTACLE_RADIUS) dist = 0.2*OBSTACLE_RADIUS; + + accel = KSPRING_OBSTACLE_COUPLE*(obstacle[i].eqdist[j] - dist)/(dist*obstacle[i].mass); + if (OBSTACLE_HARDCORE) + { + if (dist < 2.0*OBSTACLE_RADIUS) + accel += KSPRING_OBSTACLE_HARDCORE*(2.0*OBSTACLE_RADIUS - dist)/(obstacle[i].mass); + } + obstacle[i].vx += DT_PARTICLE*(obstacle[i].xc - obstacle[m].xc)*accel; + obstacle[i].vy += DT_PARTICLE*(obstacle[i].yc - obstacle[m].yc)*accel; + obstacle[i].energy += 0.5*KSPRING_OBSTACLE_COUPLE*(obstacle[i].eqdist[j] - dist)*(obstacle[i].eqdist[j] - dist); + } +// printf("2\n"); + if (obstacle[i].pinned) + { + obstacle[i].vx -= DT_PARTICLE*KSPRING_OBSTACLE_OSC*(obstacle[i].xc - obstacle[i].xc0)/obstacle[i].mass; + obstacle[i].vy -= DT_PARTICLE*KSPRING_OBSTACLE_OSC*(obstacle[i].yc - obstacle[i].yc0)/obstacle[i].mass; + obstacle[i].energy += KSPRING_OBSTACLE_OSC*((obstacle[i].xc - obstacle[i].xc0)*(obstacle[i].xc - obstacle[i].xc0) + (obstacle[i].yc - obstacle[i].yc0)*(obstacle[i].yc - obstacle[i].yc0)); + } + } + for (i = 0; i < nobstacles; i++) + { + obstacle[i].xc += DT_PARTICLE*obstacle[i].vx; + obstacle[i].yc += DT_PARTICLE*obstacle[i].vy; + } + for (i = 0; i < nobstacles; i++) + { + ekin = obstacle[i].mass*(obstacle[i].vx*obstacle[i].vx + obstacle[i].vy*obstacle[i].vy); + obstacle[i].energy += ekin; + } + } + else + { + for (i = 0; i < nobstacles; i++) + { + obstacle[i].vx += DT_PARTICLE*(obstacle[i].fx - KSPRING_OBSTACLE_OSC*(obstacle[i].xc - obstacle[i].xc0))/obstacle[i].mass; + obstacle[i].vy += DT_PARTICLE*(obstacle[i].fy - KSPRING_OBSTACLE_OSC*(obstacle[i].yc - obstacle[i].yc0))/obstacle[i].mass; + obstacle[i].xc += DT_PARTICLE*obstacle[i].vx; + obstacle[i].yc += DT_PARTICLE*obstacle[i].vy; + } + for (i = 0; i < nobstacles; i++) + { + ekin = obstacle[i].mass*(obstacle[i].vx*obstacle[i].vx + obstacle[i].vy*obstacle[i].vy); + epot = KSPRING_OBSTACLE_OSC*((obstacle[i].xc - obstacle[i].xc0)*(obstacle[i].xc - obstacle[i].xc0) + (obstacle[i].yc - obstacle[i].yc0)*(obstacle[i].yc - obstacle[i].yc0)); + etot = ekin + epot; + obstacle[i].energy = etot; + } + } +// printf("Evolved %i obstacles\n", nobstacles); +} 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, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, - angle, theta, sum, alpha, bfield; + a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, angle, theta, sum, alpha, bfield, track_x0, track_y0, efield, efieldy; double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds; double *cqx, *cqy, *cpx, *cpy, *cqangle, *cpangle; int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, cl, min_nb, max_nb, close, wrapx = 0, wrapy = 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, ncollisions = 0, ncoupled = 1, np, belt, obs; + group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt, obs, reset_obs; int *particle_numbers; static int imin, imax; static short int first = 1; @@ -1589,6 +1705,8 @@ void animation() t_molecule *molecule; t_lj_parameters params; t_belt *conveyor_belt; + t_otriangle *otriangle; + t_ofacet *ofacet; char message[100]; ratioc = 1.0 - ratio; @@ -1616,7 +1734,7 @@ void animation() } if (TRACER_PARTICLE) - trajectory = (t_tracer *)malloc(TRAJECTORY_LENGTH*N_TRACER_PARTICLES*sizeof(t_tracer)); + trajectory = (t_tracer *)malloc(TRAJECTORY_LENGTH*N_TRACER_PARTICLES*TRACER_STEPS*sizeof(t_tracer)); if (PAIR_PARTICLES) molecule = (t_molecule *)malloc(NMAXCIRCLES*sizeof(t_molecule)); /* molecules */ @@ -1641,6 +1759,12 @@ void animation() cpangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); } + if (FILL_OBSTACLE_TRIANGLES) + { + otriangle = (t_otriangle *)malloc(NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES*sizeof(t_otriangle)); + ofacet = (t_ofacet *)malloc(NMAXOBSTACLES*sizeof(t_ofacet)); + } + if (REACTION_DIFFUSION) { collisions = (t_collision *)malloc(2*NMAXCOLLISIONS*sizeof(t_collision)); @@ -1655,7 +1779,7 @@ void animation() lj_log = fopen("lj_logfile.txt", "w"); - if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle); + if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle, otriangle, ofacet); if (ADD_FIXED_SEGMENTS) init_segment_config(segment, conveyor_belt); if ((MOVE_SEGMENT_GROUPS)&&(ADD_FIXED_SEGMENTS)) @@ -1705,7 +1829,7 @@ void animation() // } // sleep(1); - update_hashgrid(particle, hashgrid, 1); + update_hashgrid(particle, obstacle, hashgrid, 1); printf("Updated hashgrid\n"); compute_relative_positions(particle, hashgrid); printf("Computed relative positions\n"); @@ -1787,6 +1911,8 @@ void animation() pright = 0.0; if (RECORD_PRESSURES) for (j=0; j= 2)) repair_cluster(cl, particle, cluster, 1000/NVID, 1); + /* evolution of lid coordinate */ if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(params.fboundary); if (BOUNDARY_COND == BC_RECTANGLE_WALL) @@ -2020,11 +2159,43 @@ void animation() } } + if (OSCILLATE_OBSTACLES) evolve_obstacles(obstacle); + if ((MOVE_SEGMENT_GROUPS)&&(i > INITIAL_TIME + SEGMENT_DEACTIVATION_TIME)) evolve_segment_groups(segment, i, segment_group); -// if ((MOVE_SEGMENT_GROUPS)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segment_groups(segment, i, segment_group); + /* update tracer particle trajectory */ + if ((TRACER_PARTICLE)&&(i > INITIAL_TIME)&&(n%(NVID/TRACER_STEPS)==0)) + { + for (j=0; j= TRAJECTORY_LENGTH*TRACER_STEPS) traj_position = 0; + traj_length++; + if (traj_length >= TRAJECTORY_LENGTH*TRACER_STEPS) + traj_length = TRAJECTORY_LENGTH*TRACER_STEPS - 1; + } } /* end of for (n=0; n INITIAL_TIME)) - { - 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 INITIAL_TIME)) +// { +// 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 ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) { @@ -2215,7 +2400,7 @@ void animation() /* compute force on segments */ if (PRINT_SEGMENTS_FORCE) compute_segments_force(¶ms, segment); - update_hashgrid(particle, hashgrid, 1); + update_hashgrid(particle, obstacle, hashgrid, 1); if ((i > INITIAL_TIME + WALL_TIME)&&(PRINT_ENTROPY)) { @@ -2236,7 +2421,7 @@ void animation() draw_frame(i, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, wall, pressure, pleft, pright, particle_numbers, 1, params, particle, cluster, - collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n); + collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet); if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers(); @@ -2268,7 +2453,7 @@ void animation() { draw_frame(i, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, - collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n); + collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet); glutSwapBuffers(); save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); counter++; @@ -2309,7 +2494,7 @@ void animation() blank(); draw_frame(NSTEPS, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, - collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n); + collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet); } if (DOUBLE_MOVIE) for (i=0; i max) max = n; } + + if (OSCILLATE_OBSTACLES) + { + for (i=0; i 0)&&(npoints < nmax)) + { + /* randomly select an active circle */ + i = rand()%(npoints); + while (!active_poisson[i]) i = rand()%(npoints); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= dpoisson*dpoisson); + /* new circle is in domain */ + far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin); + } + if (far) /* accept new circle */ + { + printf("New npoint at (%.3f,%.3f) accepted\n", x, y); + point[npoints].xc = x; + point[npoints].yc = y; + active_poisson[npoints] = 1; + npoints++; + n_p_active++; + naccepted++; + } + } + if (naccepted == 0) /* inactivate circle i */ + { + active_poisson[i] = 0; + n_p_active--; + } + printf("%i active npoints\n", n_p_active); + } + + printf("Generated %i points\n", npoints); + return(npoints); +} + void init_particle_config(t_particle particles[NMAXCIRCLES]) /* initialise particle configuration */ { @@ -887,6 +951,7 @@ void init_particle_config(t_particle particles[NMAXCIRCLES]) } case (C_POISSON_DISC): { + /* TODO: use generate_poisson_discs */ printf("Generating Poisson disc sample\n"); /* generate first circle */ // particles[0].xc = LAMBDA*(2.0*(double)rand()/RAND_MAX - 1.0); @@ -1125,7 +1190,7 @@ void add_particle_config(t_particle particles[NMAXCIRCLES], double xmin, double } break; } - case (C_HEX): /* TODO */ + case (C_HEX): { dx = (xmax - xmin)/((double)NGRIDX); dy = (ymax - ymin)/((double)NGRIDY); @@ -1253,12 +1318,488 @@ void add_obstacle(double x, double y, double radius, t_obstacle obstacle[NMAXOBS else printf("Warning: NMAXOBSTACLES should be increased\n"); } +void sort_angles(double *angle, int n, int *table) +/* writes to table positions of angles in increasing order */ +{ + int i, j, temp; + double min, temp_angle[NMAX_OBSTACLE_NEIGHBOURS]; + + for (i=0; i angle[table[j+1]]) + { + temp = table[j]; + table[j] = table[j+1]; + table[j+1] = temp; + } + } + } + + /* TEST */ + for (i=0; i 0.0); +} + +int old_triangle_overlap(t_otriangle triangle1, t_otriangle triangle2, t_obstacle obstacle[NMAXOBSTACLES]) +/* returns 1 if triangles overlap */ +{ + int j, jplus, jminus; + + /* find first vertex equal to triangle1.i[0] */ + j = 0; + while (triangle2.i[j] != triangle1.i[0]) j++; + + if (j < 3) + { + jplus = (j+1)%3; + jminus = (j+5)%3; + if (triangle1.i[1] == triangle2.i[jplus]) + return(triangle_vertex_overlap(0, 1, 2, j, jplus, jminus, triangle1, triangle2, obstacle)); + else if (triangle1.i[2] == triangle2.i[jminus]) + return(triangle_vertex_overlap(2, 0, 1, j, jminus, jplus, triangle1, triangle2, obstacle)); + else return(0); + } + else + { + if ((triangle1.i[1] == triangle2.i[1])&&(triangle1.i[2] == triangle2.i[2])) + return(triangle_vertex_overlap(1, 2, 0, 1, 2, 0, triangle1, triangle2, obstacle)); + else return(0); + } +} + +int triangle_overlap(t_otriangle triangle1, t_otriangle triangle2) +/* returns 1 if triangles overlap */ +{ + int j, jplus, jminus; + +// printf("Testing triangles (%i, %i, %i) and (%i, %i, %i)\n", triangle1.i[0], triangle1.i[1], triangle1.i[2], triangle2.i[0], triangle2.i[1], triangle2.i[2]); + + /* find first vertex equal to triangle1.i[0] */ + j = 0; + while ((j<3)&&(triangle2.i[j] != triangle1.i[0])) j++; + + if (j < 3) + { + jplus = (j+1)%3; + jminus = (j+5)%3; + if (triangle1.i[1] == triangle2.i[jplus]) + return(1); + else if (triangle1.i[2] == triangle2.i[jminus]) + return(1); + else return(0); + } + else + { + if ((triangle1.i[1] == triangle2.i[1])&&(triangle1.i[2] == triangle2.i[2])) + return(1); + else return(0); + } +} + +int init_obstacle_facets(t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES]) +/* group triangles into facets, for option COUPLE_OBSTACLES */ +{ + int i, j, k, n, ik, f, nt; + + printf("Initializing obstacle facets from %i triangles\n", n_otriangles); +// printf("1\n"); + +// printf("2\n"); + for (i=0; i NMAXOBSTACLES) + { + printf("NMAXOBSTACLES should be at least %i\n", n_ofacets); + exit(1); + } + + return(n_ofacets); +} + +double otriangle_area(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle triangle) +{ + int s1, s2, s3; + + s1 = triangle.i[0]; + s2 = triangle.i[1]; + s3 = triangle.i[2]; + return(vabs(triangle_area(obstacle[s1].xc, obstacle[s1].yc, obstacle[s2].xc, obstacle[s2].yc, obstacle[s3].xc, obstacle[s3].yc))); +} + + +int init_obstacle_coupling(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES]) +/* initialize list of neighbours, for option COUPLE_OBSTACLES */ +{ + int i, j, n, m, table[NMAX_OBSTACLE_NEIGHBOURS], neigh_table[NMAX_OBSTACLE_NEIGHBOURS]; + double dist, bdistx, bdisty, angle[NMAX_OBSTACLE_NEIGHBOURS], dist_table[NMAX_OBSTACLE_NEIGHBOURS]; + short int near_boundary; + + printf("Initializing obstacle neighbours\n"); + n_otriangles = 0; + + for (i=0; i= BCXMAX - bdistx)||(obstacle[i].yc <= BCYMIN + bdisty)||(obstacle[i].yc >= BCYMAX - bdisty)); + switch (OBSTACLE_PINNING_TYPE) { + case (OP_NPINNED): + { + obstacle[i].pinned = (obstacle[i].nneighb <= NMAX_OBSTACLE_PINNED); + break; + } + case (OP_LEFT): + { + obstacle[i].pinned = (obstacle[i].xc <= BCXMIN + bdistx); + break; + } + case (OP_LEFTTOPBOT): + { + obstacle[i].pinned = near_boundary; + break; + } + case (OP_CORNERS): + { + obstacle[i].pinned = (((obstacle[i].xc <= BCXMIN + bdistx)||(obstacle[i].xc >= BCXMAX - bdistx))&&((obstacle[i].yc <= BCYMIN + bdisty)||(obstacle[i].yc >= BCYMAX - bdisty))); + break; + } + case (OP_LEFTCORNERS): + { + obstacle[i].pinned = ((near_boundary)&&((obstacle[i].xc < 0.0)||(obstacle[i].xc > BCXMAX - bdistx))); + break; + } + case (OP_BNDRY_STEP): + { + obstacle[i].pinned = ((near_boundary)&&(obstacle[i].chessboard == 0)); + } + } + + } + + printf("Found %i facets\n", n_ofacets); +// for (i=0; i OBSXMAX) obstacle[n].xc += (OBSXMAX - OBSXMIN); + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + n++; + } + nobstacles = n; + break; + } + case (O_SQUARE): + { + dx = (OBSXMAX - OBSXMIN)/(double)NOBSX; + dy = (OBSYMAX - OBSYMIN)/(double)NOBSY; + n = 0; + + if ((ADD_FIXED_SEGMENTS)&&(SEGMENT_PATTERN == S_CYLINDER)) + /* pseudo-cylindrical boundary conditions (e.g. for Hall effect) */ + { + jmin = 1; + jmax = NOBSY-1; + } + else + { + jmin = 0; + jmax = NOBSY; + } + + for (i=0; i OBSXMAX) obstacle[n].xc += (OBSXMAX - OBSXMIN); + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + obstacle[n].chessboard = (i+j)%BDRY_PINNING_STEP; + n++; + } + nobstacles = n; + break; + } + case (O_SQUARE_TWOMASSES): { dx = (XMAX - XMIN)/(double)NOBSX; dy = (YMAX - YMIN)/(double)NOBSY; @@ -1413,15 +2039,18 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) { obstacle[n].xc = XMIN + ((double)i + 0.25)*dx; obstacle[n].yc = YMIN + ((double)j + 0.5)*dy; - if (j%2 == 1) obstacle[n].xc += 0.5*dx; if (obstacle[n].xc > XMAX) obstacle[n].xc += (XMAX - XMIN); obstacle[n].radius = OBSTACLE_RADIUS; + if ((i+j)%2 == 0) + { + obstacle[n].radius *= sqrt(2.0); + obstacle[n].mass *= 2.0; + } obstacle[n].active = 1; + obstacle[n].chessboard = (i+j)%BDRY_PINNING_STEP; n++; } nobstacles = n; -// printf("Added %i obstacles\n", nobstacles); -// for (n=0; n XMAX) obstacle[n].xc += (XMAX - XMIN); + obstacle[n].radius = OBSTACLE_RADIUS; + obstacle[n].active = 1; + } + free(point); + break; + } + default: { printf("Function init_obstacle_config not defined for this pattern \n"); @@ -1624,13 +2300,16 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES]) obstacle[n].omega0 = obstacle[n].omega; obstacle[n].xc0 = obstacle[n].xc; obstacle[n].yc0 = obstacle[n].yc; + obstacle[n].nneighb = 0; } + + if (COUPLE_OBSTACLES) init_obstacle_coupling(obstacle, otriangle, ofacet); } void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, int center, t_segment segment[NMAXSEGMENTS], int group) /* add four segments forming a rectangle, specified by two adjacent corners and width */ { - double tx, ty, ux, uy, norm, x3, y3, x4, y4; + double tx, ty, ux, uy, norm, x3, y3, x4, y4, angle; int i, n = nsegments; tx = x2 - x1; @@ -1671,9 +2350,25 @@ void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, d segment[n+3].y1 = y4; segment[n+3].x2 = x1; segment[n+3].y2 = y1; + + segment[n].angle1 = -PID; + segment[n].angle2 = 0.0; + segment[n+1].angle1 = PI; + segment[n+1].angle2 = 1.5*PI; + + segment[n+2].angle1 = PID; + segment[n+2].angle2 = PI; + + segment[n+3].angle1 = 0.0; + segment[n+3].angle2 = PID; + + angle = argument(tx, ty) + PI; + for (i=0; i<4; i++) { + segment[n+i].angle1 += angle; + segment[n+i].angle2 += angle; segment[n+i].concave = 1; segment[n+i].group = group; } @@ -3440,6 +4135,43 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[N break; } + case (S_TREES_B): + { + for (i=0; i PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN; - else if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX; + hue = SEGMENT_HUE_MIN + (SEGMENT_HUE_MAX - SEGMENT_HUE_MIN)*pressure/SEGMENT_PMAX; + if (hue > SEGMENT_HUE_MIN) hue = SEGMENT_HUE_MIN; + else if (hue < SEGMENT_HUE_MAX) hue = SEGMENT_HUE_MAX; // hsl_to_rgb_turbo(hue, 0.9, lum, rgb); @@ -7055,7 +7843,7 @@ void draw_one_triangle(t_particle particle, int same_table[9*HASHMAX], int p, in x = particle.xc + (double)p*(BCXMAX - BCXMIN); y = particle.yc + (double)q*(BCYMAX - BCYMIN); - if (TRACK_SEGMENT_GROUPS) + if (TRACK) { x -= xtrack; y -= ytrack; @@ -7542,7 +8330,7 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, if (CENTER_VIEW_ON_OBSTACLE) xc1 = xc - xshift; else xc1 = xc; - if (TRACK_SEGMENT_GROUPS) + if (TRACK) { xc1 -= xtrack; yc1 = yc - ytrack; @@ -7663,7 +8451,7 @@ void draw_collisions(t_collision *collisions, int ncollisions) if (CENTER_VIEW_ON_OBSTACLE) x1 = x - xshift; else x1 = x; - if (TRACK_SEGMENT_GROUPS) + if (TRACK) { x1 -= xtrack; y1 = y - ytrack; @@ -7686,72 +8474,155 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], glLineWidth(TRAJECTORY_WIDTH); printf("drawing trajectory\n"); - for (j=0; j= 0; i--) + for (j=0; j 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) + draw_line(x1, y1, x2, y2); + } + else + { +// for (i = traj_length-2; i >= traj_position + 1; i--) + for (i = traj_position + 1; i < traj_length-1; i++) + for (j=0; j 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) + draw_line(x1, y1, x2, y2); + } + for (i=0; i < traj_position-1; i++) +// for (i = traj_position-2; i >= 0; i--) + for (j=0; j 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) + draw_line(x1, y1, x2, y2); + } + } +} + +void old_draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, int *tracer_n, int plot) +/* draw tracer particle trajectory */ +{ + int i, j, time, p, width; + double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum; + +// blank(); + glLineWidth(TRAJECTORY_WIDTH); + printf("drawing trajectory\n"); + + for (j=0; j= 0; i--) { - x1 = trajectory[j*TRAJECTORY_LENGTH + i].xc; - x2 = trajectory[j*TRAJECTORY_LENGTH + i+1].xc; - y1 = trajectory[j*TRAJECTORY_LENGTH + i].yc; - y2 = trajectory[j*TRAJECTORY_LENGTH + i+1].yc; + x1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].xc; + x2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].xc; + y1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].yc; + y2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].yc; time = traj_length - i; - lum = 0.8 - TRACER_LUM_FACTOR*(double)time/(double)TRAJECTORY_LENGTH; + 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 (module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN)) draw_line(x1, y1, x2, y2); + if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) + draw_line(x1, y1, x2, y2); // printf("tracer = %i, (x1, y1) = (%.3lg, %.3lg), (x2, y2) = (%.3lg, %.3lg)\n", j, x1, y1, x2, y2); } else { - for (i = traj_position + 1; i < traj_length-1; i++) +// for (i = traj_position + 1; i < traj_length-1; i++) + for (i = traj_length-2; i >= traj_position + 1; i--) { - x1 = trajectory[j*TRAJECTORY_LENGTH + i].xc; - x2 = trajectory[j*TRAJECTORY_LENGTH + i+1].xc; - y1 = trajectory[j*TRAJECTORY_LENGTH + i].yc; - y2 = trajectory[j*TRAJECTORY_LENGTH + i+1].yc; + x1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].xc; + x2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].xc; + y1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].yc; + y2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].yc; time = traj_position + traj_length - i; - lum = 1.0 - (double)time/(double)TRAJECTORY_LENGTH; + 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 (module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN)) draw_line(x1, y1, x2, y2); + if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) + draw_line(x1, y1, x2, y2); } - for (i=0; i < traj_position-1; i++) +// for (i=0; i < traj_position-1; i++) + for (i = traj_position-2; i >= 0; i--) { - x1 = trajectory[j*TRAJECTORY_LENGTH + i].xc; - x2 = trajectory[j*TRAJECTORY_LENGTH + i+1].xc; - y1 = trajectory[j*TRAJECTORY_LENGTH + i].yc; - y2 = trajectory[j*TRAJECTORY_LENGTH + i+1].yc; + x1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].xc; + x2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].xc; + y1 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i].yc; + y2 = trajectory[j*TRAJECTORY_LENGTH*TRACER_STEPS + i+1].yc; time = traj_position - i; - lum = 1.0 - (double)time/(double)TRAJECTORY_LENGTH; + 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 (module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN)) draw_line(x1, y1, x2, y2); + if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) + draw_line(x1, y1, x2, y2); } } } } - -void color_background(t_particle particle[NMAXCIRCLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY]) +void color_background(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY]) /* color background according to particle properties */ { - int i, j, k, n, p, q, m, nnb, number, avrg_fact; - double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue; + int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs; + double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue, valx, valy, lum; static int first = 1, counter = 0; p1 = 0.75; @@ -7851,6 +8722,119 @@ void color_background(t_particle particle[NMAXCIRCLES], int bg_color, t_hashgrid hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); break; } + case (BG_EOBSTACLES): + { + value = 0.0; + nnb = hashgrid[n].nneighb; + for (q=0; q ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; + if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; + hue = p1*oldhue + p2*hue; + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); + break; + } + case (BG_EKIN_OBSTACLES): + { + value = 0.0; + nnb = hashgrid[n].nneighb; + for (q=0; q ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; + if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; + hue = p1*oldhue + p2*hue; + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); + break; + } + case (BG_DIR_OBSTACLES): + { + valx = 0.0; + valy = 0.0; + nnb = hashgrid[n].nneighb; + for (q=0; q DPI) hue -= DPI; + if (hue < 0.0) hue += DPI; + hue = pp1*oldhue + pp2*hue; + lum = module2(valx, valy)/OBSTACLE_VMAX; + if (lum > 1.0) lum = 1.0; + hsl_to_rgb_palette(hue*360.0/DPI, 0.9, 0.5*lum, rgb, COLOR_PALETTE_DIRECTION); + break; + } + case (BG_POS_OBSTACLES): + { + valx = 0.0; + valy = 0.0; + nnb = hashgrid[n].nneighb; +// for (q=0; q DPI) hue -= DPI; + if (hue < 0.0) hue += DPI; + hue = pp1*oldhue + pp2*hue; + lum = module2(valx, valy); + if (lum > 1.0) lum = 1.0; +// lum = 1.0; + hsl_to_rgb_palette(hue*360.0/DPI, 0.9, 0.5*lum, rgb, COLOR_PALETTE_DIRECTION); + break; + } case (BG_FORCE): { value = 0.0; @@ -7902,12 +8886,15 @@ void draw_particles(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRC 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]; + printf("Drawing particles\n"); // if (!TRACER_PARTICLE) blank(); if (plot == P_NOPARTICLE) blank(); - if ((COLOR_BACKGROUND)&&(bg_color > 0)) color_background(particle, bg_color, hashgrid); - else if (!TRACER_PARTICLE) blank(); +// if ((COLOR_BACKGROUND)&&(bg_color > 0)) color_background(particle, bg_color, hashgrid); +// else + if ((!TRACER_PARTICLE)&&(!COLOR_BACKGROUND)&&(!DRAW_OBSTACLE_LINKS)&&(!FILL_OBSTACLE_TRIANGLES)) blank(); +// printf("Drawing particles 1\n"); glColor3f(1.0, 1.0, 1.0); /* show region of partial thermostat */ @@ -7986,6 +8973,8 @@ void draw_particles(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRC } } } +// printf("Drawing particles 2\n"); + /* draw "altitude lines" */ if (ALTITUDE_LINES) draw_altitude_lines(); @@ -8195,12 +9184,254 @@ void draw_particles(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRC } } +void obstacle_hue(t_obstacle obstacle, double rgb[3]) +/* compute obstacle color hue */ +{ + int k; + double etot, angle, hue, lum; + + switch (OBSTACLE_COLOR) { + case (OC_ENERGY): + { + etot = obstacle.energy; + hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*etot/OBSTACLE_EMAX; + if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; + if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN); + break; + } + case (OC_DIRECTION): + { + hue = argument(obstacle.vx, obstacle.vy); + if (hue < 0.0) hue += DPI; + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIRECTION); + break; + } + case (OC_DIRECTION_ENERGY): + { + hue = argument(obstacle.vx, obstacle.vy); + etot = obstacle.energy; + if (hue < 0.0) hue += DPI; + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; + lum = etot/OBSTACLE_EMAX; + if (lum > 0.5) lum = 0.5; + hsl_to_rgb_palette(hue, 0.9, lum, rgb, COLOR_PALETTE_DIRECTION); + break; + } + } +} + +int area_to_rgb(t_obstacle obstacle[NMAXOBSTACLES], int i, int n1, int n2, double rgb[3]) +/* computes rgb code of triangle according to area */ +/* returns 1 if angle of triangle is less than PI/2 */ +/* OLD VERSION */ +{ + int p, q, k, n, counter; + double area, angle1, angle2; + static int first; + static double mean_area; + + if (first) + { + /* compute mean area of triangles */ + if (FILL_OBSTACLE_TRIANGLES) + { + counter = 0; + mean_area = 0.0; + for (p = 0; p < nobstacles; p++) + { + n = obstacle[p].nneighb; + for (q = 0; q < n - 1; q++) + { + n1 = obstacle[p].neighb[q]; + n2 = obstacle[p].neighb[q+1]; + mean_area += vabs(triangle_area(obstacle[p].xc, obstacle[p].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc)); + counter++; + } + n1 = obstacle[p].neighb[n-1]; + n2 = obstacle[p].neighb[0]; + mean_area += vabs(triangle_area(obstacle[p].xc, obstacle[p].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc)); + counter++; + } + mean_area *= 1.0/(double)counter; + } + first = 0; + } + + area = vabs(triangle_area(obstacle[i].xc, obstacle[i].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc)); + area = OBSTACLE_AREA_SHADE_FACTOR*(area - mean_area); + if (area > 1.0) area = 1.0; + if (area < 0.0) area = 0.0; + for (k=0; k<3; k++) rgb[k] = area; + + angle1 = argument(obstacle[n1].xc - obstacle[i].xc, obstacle[n1].yc - obstacle[i].yc); + angle2 = argument(obstacle[n2].xc - obstacle[i].xc, obstacle[n2].yc - obstacle[i].yc); + if (angle2 < angle1) angle2 += DPI; + + return(angle2 - angle1 < 1.5*PID); +} + +void old_draw_obstacle_triangles(t_obstacle obstacle[NMAXOBSTACLES]) +/* draw the triangles between obstacles, for option FILL_OBSTACLE_TRIANGLES */ +{ + int i, j, k, n, n1, n2, neigh; + short int acute; + double rgb[3], angle1, angle2; + + if (FILL_OBSTACLE_TRIANGLES) for (i = 0; i < nobstacles; i++) + { + n = obstacle[i].nneighb; + for (j = 0; j < n - 1; j++) + { + n1 = obstacle[i].neighb[j]; + n2 = obstacle[i].neighb[j+1]; + acute = area_to_rgb(obstacle, i, n1, n2, rgb); + if (acute) draw_colored_triangle(obstacle[i].xc, obstacle[i].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc, rgb); + } + n1 = obstacle[i].neighb[n-1]; + n2 = obstacle[i].neighb[0]; + acute = area_to_rgb(obstacle, i, n1, n2, rgb); + if (acute) draw_colored_triangle(obstacle[i].xc, obstacle[i].yc, obstacle[n1].xc, obstacle[n1].yc, obstacle[n2].xc, obstacle[n2].yc, rgb); + } + glColor3f(1.0, 1.0, 1.0); + for (i = 0; i < nobstacles; i++) + for (j = 0; j < obstacle[i].nneighb; j++) + { + neigh = obstacle[i].neighb[j]; + draw_line(obstacle[i].xc, obstacle[i].yc, obstacle[neigh].xc, obstacle[neigh].yc); + } +} + +void facet_area_to_rgb(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle triangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet facet[NMAXOBSTACLES], int nfacet, double rgb[3]) +/* computes rgb code of triangle according to area */ +/* returns 1 if angle of triangle is less than PI/2 */ +{ + int p, q, k, n, counter; + double area, area0; + static int first; + static double mean_area; + + if (first) + { + /* compute mean area of triangles */ + if (FILL_OBSTACLE_TRIANGLES) + { + counter = 0; + mean_area = 0.0; + for (p = 0; p < n_otriangles; p++) + { + mean_area += otriangle_area(obstacle, triangle[p]); + counter++; + } + mean_area *= 1.0/(double)counter; + + for (p = 0; p < n_ofacets; p++) + { + area0 = 0.0; + for (k = 0; k < facet[p].ntriangles; k++) + area0 += triangle[facet[p].triangle[k]].area0; + facet[p].area0 = area0/(double)facet[p].ntriangles; + } + } + first = 0; + } + + area = 0.0; + for (k = 0; k < facet[nfacet].ntriangles; k++) + { + p = facet[nfacet].triangle[k]; + area += otriangle_area(obstacle, triangle[p]); + } + area *= 1.0/(double)facet[nfacet].ntriangles; + + area = OBSTACLE_AREA_SHADE_FACTOR*(area - facet[nfacet].area0); + if (area > 1.0) area = 1.0; + if (area < 0.0) area = 0.0; + for (k=0; k<3; k++) rgb[k] = area; +} + +double triangle_area_to_rgb(t_obstacle obstacle[NMAXOBSTACLES], t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], int triangle, double rgb[3]) +{ + double area; + int k, p; + + area = otriangle_area(obstacle, otriangle[triangle]); + area = 0.6 + OBSTACLE_AREA_SHADE_FACTOR*(area - otriangle[triangle].area0); + if (area > 0.9) area = 0.9; + if (area < 0.1) area = 0.1; + for (k=0; k<3; k++) rgb[k] = area; + + return(area); +} + +void draw_obstacle_triangles(t_obstacle obstacle[NMAXOBSTACLES], + t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], + t_ofacet ofacet[NMAXOBSTACLES]) +/* draw the triangles between obstacles, for option FILL_OBSTACLE_TRIANGLES */ +{ + int i, j, k, n, n1, n2, neigh, nf, f, t, s1, s2, s3; + short int acute; + double rgb[3], pvect; + + printf("drawing %i facets with %i triangles\n", n_ofacets, n_otriangles); + if (FILL_OBSTACLE_TRIANGLES) + { + /* determine acute angles */ + for (t = 0; t < n_otriangles; t++) + { + s1 = otriangle[t].i[0]; + s2 = otriangle[t].i[1]; + s3 = otriangle[t].i[2]; + + pvect = (obstacle[s2].xc - obstacle[s1].xc)*(obstacle[s3].yc - obstacle[s1].yc); + pvect -= (obstacle[s2].yc - obstacle[s1].yc)*(obstacle[s3].xc - obstacle[s1].xc); + + otriangle[t].acute = (pvect > 0.0); + } + + if (SHADE_OBSTACLE_FACETS) for (i = 0; i < n_ofacets; i++) + { + facet_area_to_rgb(obstacle, otriangle, ofacet, i, rgb); + nf = ofacet[i].ntriangles; + for (f = 0; f < nf; f++) + { + t = ofacet[i].triangle[f]; + if (otriangle[t].acute) + { + s1 = otriangle[t].i[0]; + s2 = otriangle[t].i[1]; + s3 = otriangle[t].i[2]; + + draw_colored_triangle(obstacle[s1].xc, obstacle[s1].yc, obstacle[s2].xc, obstacle[s2].yc, obstacle[s3].xc, obstacle[s3].yc, rgb); + } + } + } + else for (t = 0; t < n_otriangles; t++) if (otriangle[t].acute) + { + triangle_area_to_rgb(obstacle, otriangle, t, rgb); + s1 = otriangle[t].i[0]; + s2 = otriangle[t].i[1]; + s3 = otriangle[t].i[2]; + + draw_colored_triangle(obstacle[s1].xc, obstacle[s1].yc, obstacle[s2].xc, obstacle[s2].yc, obstacle[s3].xc, obstacle[s3].yc, rgb); + } + } + + glColor3f(1.0, 1.0, 1.0); + for (i = 0; i < nobstacles; i++) + for (j = 0; j < obstacle[i].nneighb; j++) + { + neigh = obstacle[i].neighb[j]; + draw_line(obstacle[i].xc, obstacle[i].yc, obstacle[neigh].xc, obstacle[neigh].yc); + } +} void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], t_belt *conveyor_belt, int wall) /* draw the container, for certain boundary conditions */ { - int i, j; - double rgb[3], x, y, r, phi, angle, dx, dy, ybin, x1, x2, y1, y2, h, bangle, blength, pos0, pos, bsegment, tx, ty, bwidth, ca, sa; + int i, j, k, n, n1, n2, neigh; + double rgb[3], x, y, r, phi, angle, dx, dy, ybin, x1, x2, y1, y2, h, bangle, blength, pos0, pos, bsegment, tx, ty, bwidth, ca, sa, ekin, epot, etot, hue, area; char message[100]; /* draw fixed obstacles */ @@ -8208,9 +9439,16 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] { glLineWidth(CONTAINER_WIDTH); if (CHARGE_OBSTACLES) hsl_to_rgb(30.0, 0.1, 0.5, rgb); - else hsl_to_rgb(300.0, 0.1, 0.5, rgb); + else if (!OSCILLATE_OBSTACLES) hsl_to_rgb(300.0, 0.1, 0.5, rgb); for (i = 0; i < nobstacles; i++) + { + + if (OSCILLATE_OBSTACLES) + { + obstacle_hue(obstacle[i], rgb); + } draw_colored_circle_precomp(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius, rgb); + } glColor3f(1.0, 1.0, 1.0); for (i = 0; i < nobstacles; i++) { @@ -8256,14 +9494,14 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] dx = r*cos(bangle); dy = r*sin(bangle); - x1 = conveyor_belt[i].x1; - y1 = conveyor_belt[i].y1; + x1 = conveyor_belt[i].x1 - xtrack; + y1 = conveyor_belt[i].y1 - ytrack; draw_line(x1-dx, y1-dy, x1+dx, y1+dy); draw_line(x1+dy, y1-dx, x1-dy, y1+dx); draw_circle_precomp(x1, y1, r); - x2 = conveyor_belt[i].x2; - y2 = conveyor_belt[i].y2; + x2 = conveyor_belt[i].x2 - xtrack; + y2 = conveyor_belt[i].y2 - ytrack; draw_line(x2-dx, y2-dy, x2+dx, y2+dy); draw_line(x2+dy, y2-dx, x2-dy, y2+dx); draw_circle_precomp(x2, y2, r); @@ -9231,10 +10469,13 @@ void print_particles_speeds(t_particle particle[NMAXCIRCLES]) write_text(xrighttext + 0.1, y, message); } -double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall, double krepel) +double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall, double krepel, int reset) { - int i, k, corner; + int i, s, k, corner, group, ngroups; double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, y1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin, charge, d2, speed, pangle, distance1, sangle; + double group_pressure[NMAXGROUPS]; + int segments_per_group[NMAXGROUPS]; + short int inactivate_group[NMAXGROUPS], inactivate_segment; static int first = 1; static double seqangle; @@ -9244,6 +10485,12 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl first = 0; } + if ((OSCILLATE_OBSTACLES)&&(reset)) for (i=0; i ngroups) ngroups = group; + } + for (group=1; group 0) + { +// group_pressure[group] *= 1.0/(double)segments_per_group[group]; + if (group_pressure[group] > SEGMENT_P_INACTIVATE*(double)segments_per_group[group]) inactivate_group[group] = 1; + else inactivate_group[group] = 0; + } + for (i=0; i 0)&&(inactivate_group[group])) segment[i].active = 0; + } + } switch(BOUNDARY_COND){ case (BC_SCREEN): @@ -9829,6 +11110,30 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl else if (x < BCXMIN + padding) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - x + padding); } + return(0.0); + } + case (BC_REFLECT_ABS_RIGHT): + { + /* add harmonic force outside screen */ + padding = MU; + x = particle[j].xc; + y = particle[j].yc; + + if (x > BCXMAX + padding) + { + particle[j].active = 0; + particle[j].vx = 0.0; + particle[j].vy = 0.0; + particle[j].xc = BCXMAX + 2.0*padding; + particle[j].yc = BCYMIN - 2.0*padding; + } + else + { + if (y > BCYMAX - padding) particle[j].fy -= KSPRING_BOUNDARY*(y - BCYMAX + padding); + else if (y < BCYMIN + padding) particle[j].fy += KSPRING_BOUNDARY*(BCYMIN - y + padding); + if (x < BCXMIN + padding) particle[j].fx += KSPRING_BOUNDARY*(BCXMIN - x + padding); + } + return(0.0); } } @@ -10714,8 +12019,14 @@ int tracer_n[N_TRACER_PARTICLES]) } } - add_particle(x, y, 0.0, V_INITIAL, PARTICLE_MASS, type, particle); + 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); + + printf("Particle number %i\n", ncircles-1); + +// add_particle(x, y, V_INITIAL, 0.25*V_INITIAL*(double)rand()/RAND_MAX, PARTICLE_MASS, type, particle); // add_particle(x, y, V_INITIAL*(double)rand()/RAND_MAX, 0.0, PARTICLE_MASS, type, particle); // if (y > 0.0) @@ -10921,6 +12232,16 @@ int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin, t } +int partial_efield(double x, double y) +/* */ +{ + switch (EFIELD_REGION) { + case (EF_CONST): return(1); + case (EF_SQUARE): return(vabs(x) < YMAX); + case (EF_LEFT): return(x < YMIN); + } +} + int partial_bfield(double x, double y) /* */ { @@ -15303,9 +16624,15 @@ void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_positio t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], t_group_data *group_speeds, t_group_segments *segment_group, - t_belt *conveyor_belt, int *tracer_n) + t_belt *conveyor_belt, int *tracer_n, + t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES]) /* draw a movie frame */ { + printf("Drawing frame\n"); + if ((COLOR_BACKGROUND)&&(bg_color > 0)) color_background(particle, obstacle, bg_color, hashgrid); +// else if (!TRACER_PARTICLE) blank(); + if (DRAW_OBSTACLE_LINKS) draw_obstacle_triangles(obstacle, otriangle, ofacet); +// if (DRAW_OBSTACLE_LINKS) old_draw_obstacle_triangles(obstacle); if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length, particle, cluster, tracer_n, plot); draw_particles(particle, cluster, plot, params.beta, collisions, ncollisions, bg_color, hashgrid, params); draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, conveyor_belt, wall); @@ -15322,5 +16649,6 @@ void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_positio } else if (PRINT_PARTICLE_SPEEDS) print_particles_speeds(particle); else if (PRINT_SEGMENTS_SPEEDS) print_segment_group_speeds(segment_group); +// printf("5\n"); } diff --git a/sub_rde.c b/sub_rde.c index 03c54f7..9b94346 100644 --- a/sub_rde.c +++ b/sub_rde.c @@ -4877,6 +4877,98 @@ void compute_laplacian_rde(double phi_in[NX*NY], double phi_out[NX*NY], short in } } +void compute_kuramoto_int_planar(double phi_in[NX*NY], double phi_out[NX*NY], short int xy_in[NX*NY]) +/* computes Kuramoto interaction of phi_in and stores it in phi_out - case with whole rectangular domain */ +{ + int i, j, iplus, iminus, jplus, jminus; + double phi; + + #pragma omp parallel for private(i,j) + for (i=1; i nmaxpixels) { - printf("Image too large, increase nmaxpixels in init_xyin_from_image()\n"); + printf("Image too large, increase nmaxpixels in init_xyin_from_image() to %i\n", nx*ny); exit(0); } @@ -2679,17 +2806,24 @@ int init_xyin_from_image(short int * xy_in[NX]) for (i=0; i nx-1) ii = nx-1; if (ii < 0) ii = 0; if (jj > ny-1) jj = ny-1; if (jj < 0) jj = 0; k = 3*(jj*nx+ii); rgbtot = rgb_values[k] + rgb_values[k+1] + rgb_values[k+2]; - if (rgbtot > 3*maxrgb/2) xy_in[i][j] = 1; + if (rgbtot*sign > rgbmax*sign) xy_in[i][j] = 1; else xy_in[i][j] = 0; + + if (ii == 0) xy_in[i][j] = 0; +// if (ii == nx-1) xy_in[i][j] = 0; } + + /* to avoid reflection on left boundary */ + if (IMAGE_FILE == IM_HAPPY_NEW_YEAR) + for (j=0; j -MU)&&(h < 2.0*MU)&&(r*r - h*h < WALL_WIDTH*WALL_WIDTH)) return(1); + } + return(0); + } + case (D_CARDIOID): + { + y1 = y - MU; + r = x*x + y1*y1; + return(r*r + 4.0*LAMBDA*y1*r - 4.0*LAMBDA*LAMBDA*x*x < 0.0); + } + case (D_NEPHROID): + { + return(108*ipow(LAMBDA,4)*x*x > ipow(x*x + y*y - 4.0*LAMBDA*LAMBDA,3)); + } + case (D_EPICYCLOID): + { + if (first) + { + init_epicyloid(); + first = 0; + } + phi = argument(x, y) - APOLY*PID; + if (phi < 0.0) phi += DPI; + if (phi >= DPI) phi -= DPI; + a = DPI/(double)NPOLY; + phi += PID; + while (phi > a) phi -= a; + if (phi > 0.5*a) phi = a - phi; + k = (int)(phi*(double)NEPICYCLOID/a); + return(module2(x,y) < epicycloid_phi_to_r[k]*LAMBDA/(double)(NPOLY)); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -4074,7 +4302,7 @@ void draw_rc_hyp() void draw_billiard(int fade, double fade_value) /* draws the billiard boundary */ { - double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, alpha2, dphi, omega, z, l, width, a, b, c, d, r1, r2, ymax, height, xmin, xmax, ca, sa, xshift, x5, x6, f, fp, xratio, w; + double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, alpha2, dphi, omega, z, l, width, a, b, c, d, r1, r2, ymax, height, xmin, xmax, ca, sa, xshift, x5, x6, f, fp, xratio, w, nx, ny, x3, y3, psi, psi0, psi2; int i, j, k, k1, k2, mr2, ntiles; static int first = 1, nsides; static double h, hh, sqr3, ll, salpha, arcangle; @@ -5916,6 +6144,122 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound } break; } + case (D_TWOCIRCLES_ELLIPSE): + { + a = 0.9*XMAX - 3.0*MU; + b = a/LAMBDA; + alpha2 = asin(WALL_WIDTH/MU); + alpha = asin(WALL_WIDTH/b); + x1 = a + 2.0*MU; + draw_circle_arc(x1, 0.0, MU, -PI + alpha2, DPI - 2.0*alpha2, NSEG); + draw_line(x1 - MU*cos(alpha2), WALL_WIDTH, a*cos(alpha), WALL_WIDTH); + draw_ellipse_arc(0.0, 0.0, a, b, alpha, PI - 2.0*alpha, NSEG); + draw_line(-a*cos(alpha), WALL_WIDTH, -x1 + MU*cos(alpha2), WALL_WIDTH); + draw_circle_arc(-x1, 0.0, MU, alpha2, DPI - 2.0*alpha2, NSEG); + draw_line(-x1 + MU*cos(alpha2), -WALL_WIDTH, -a*cos(alpha), -WALL_WIDTH); + draw_ellipse_arc(0.0, 0.0, a, b, PI + alpha, PI - 2.0*alpha, NSEG); + draw_line(a*cos(alpha), -WALL_WIDTH, x1 - MU*cos(alpha2), -WALL_WIDTH); + break; + } + case (D_CIRCLES_ELLIPSE): + { + a = JULIA_SCALE*(0.9*XMAX - 3.0*MU); + b = a/LAMBDA; + alpha2 = asin(WALL_WIDTH/MU); + for (k=0; k0) draw_ellipse_arc(0.0, 0.0, a, b, psi, psi2 - psi, NSEG); + + draw_line(x2 - JULIA_SCALE*MU*cos(alpha + alpha2), y2 - JULIA_SCALE*MU*sin(alpha + alpha2), x3, y3); + + draw_circle_arc(x2, y2, JULIA_SCALE*MU, alpha - PI + alpha2, DPI - 2.0*alpha2, NSEG); + + x3 = x1 - JULIA_SCALE*WALL_WIDTH*ny; + y3 = y1 + JULIA_SCALE*WALL_WIDTH*nx; + draw_line(x2 - JULIA_SCALE*MU*cos(alpha - alpha2), y2 - JULIA_SCALE*MU*sin(alpha - alpha2), x3, y3); + + psi = argument(x3/a, y3/b); + if (k == 0) psi0 = psi2; + } + if (psi0 < psi) psi0 += DPI; + draw_ellipse_arc(0.0, 0.0, a, b, psi, psi0 - psi, NSEG); + + /* draw foci */ + if (FOCI) + { + if (fade) glColor3f(0.3*fade_value, 0.3*fade_value, 0.3*fade_value); + else glColor3f(0.3, 0.3, 0.3); + x0 = sqrt(a*a - b*b); + + glLineWidth(2); + glEnable(GL_LINE_SMOOTH); + + draw_circle(x0, 0.0, 0.01, NSEG); + draw_circle(-x0, 0.0, 0.01, NSEG); + } + + break; + } + case (D_CARDIOID): + { + glBegin(GL_LINE_LOOP); + for (i=0; i nmaxpixels) + { + printf("Image too large, increase nmaxpixels in init_xyin_from_image() to %i\n", nx*ny); + exit(0); + } + + /* read rgb values */ + for (j=0; j nx-1) ii = nx-1; + if (ii < 0) ii = 0; + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + k = 3*(jj*nx+ii); + rgbtot = rgb_values[k] + rgb_values[k+1] + rgb_values[k+2]; + if (rgbtot*sign > rgbmax*sign) xy_in[i*NY+j] = 1; + else xy_in[i*NY+j] = 0; + + if (ii == 0) xy_in[i*NY+j] = 0; +// if (ii == nx-1) xy_in[i*NY+j] = 0; + } + + fclose(image_file); + free(rgb_values); + return(1); +} diff --git a/wave_billiard.c b/wave_billiard.c index 520c4f0..d77ca49 100644 --- a/wave_billiard.c +++ b/wave_billiard.c @@ -44,7 +44,7 @@ #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 when writing tiff images */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -55,29 +55,34 @@ /* General geometrical parameters */ -#define WINWIDTH 1920 /* window width */ +// #define WINWIDTH 1920 /* window width */ +#define WINWIDTH 1150 /* window width */ #define WINHEIGHT 1150 /* window height */ -#define NX 3840 /* number of grid points on x axis */ +// #define NX 3840 /* number of grid points on x axis */ +#define NX 2300 /* 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 YMIN -0.997916667 -#define YMAX 1.397916667 /* y interval for 9/16 aspect ratio */ +// #define XMIN -2.0 +// #define XMAX 2.0 /* x interval */ +#define XMIN -1.197916667 +#define XMAX 1.197916667 /* x interval */ +#define YMIN -1.197916667 +#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ #define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ #define HRES 1 /* dummy, only used by rde.c */ -#define JULIA_SCALE 1.0 /* scaling for Julia sets */ +#define JULIA_SCALE 0.8 /* scaling for Julia sets and some other domains */ /* Choice of the billiard table */ -#define B_DOMAIN 792 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 92 /* 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 IMAGE_FILE 5 /* for option D_IMAGE */ #define COMPARISON 0 /* set to 1 to compare two different patterns (beta) */ -#define B_DOMAIN_B 20 /* second domain shape, for comparisons */ +#define B_DOMAIN_B 20 /* second domain shape, for comparisons */ #define CIRCLE_PATTERN_B 0 /* second pattern of circles or polygons */ #define P_PERCOL 0.15 /* probability of having a circle in C_RAND_PERCOL arrangement */ @@ -85,11 +90,11 @@ #define PDISC_FACTOR 3.5 /* controls density of Poisson disc process (default: 3.25) */ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.85 /* parameter controlling the dimensions of domain */ -#define MU 0.18 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.8 /* parameter controlling the dimensions of domain */ +#define MU 0.75 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */ -#define NPOLY 3 /* number of sides of polygon */ -#define APOLY -0.17 /* angle by which to turn polygon, in units of Pi/2 */ +#define NPOLY 7 /* 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 */ @@ -113,7 +118,6 @@ /* You can add more billiard tables by adapting the functions */ /* xy_in_billiard and draw_billiard below */ - /* Physical parameters of wave equation */ #define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ @@ -128,10 +132,8 @@ #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.12 /* Courant number */ -// #define COURANT 0.07 /* Courant number */ -#define COURANTB 0.08 /* Courant number in medium B */ -// #define GAMMA 1.0e-7 /* damping factor in wave equation */ +#define COURANT 0.08 /* 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 */ #define GAMMA_SIDES 1.0e-4 /* damping factor on boundary */ @@ -146,10 +148,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 16 /* period of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 49 /* period of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ -#define N_SOURCES 3 /* number of sources, for option draw_sources */ +#define N_SOURCES 7 /* 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 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 */ @@ -164,7 +167,7 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2300 /* number of frames of movie */ +#define NSTEPS 1350 /* number of frames of movie */ #define NVID 8 /* 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 */ @@ -182,20 +185,20 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 1.5 /* amplitude of initial condition */ -#define INITIAL_VARIANCE 0.00001 /* variance of initial condition */ +#define INITIAL_AMP 1.4 /* amplitude of initial condition */ +#define INITIAL_VARIANCE 0.00005 /* variance of initial condition */ #define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ /* Plot type, see list in global_pdes.c */ -#define PLOT 0 +#define PLOT 8 -#define PLOT_B 8 /* plot type for second movie */ +#define PLOT_B 0 /* plot type for second movie */ /* Color schemes */ -#define COLOR_PALETTE 13 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ +#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 BLACK 1 /* background */ @@ -206,17 +209,17 @@ #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.0 /* additional shift for wave amplitude */ -#define VSCALE_AMPLITUDE 0.5 /* additional scaling factor for wave amplitude */ -#define E_SCALE 10.0 /* scaling factor for energy representation */ +#define VSHIFT_AMPLITUDE -0.5 /* 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 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.8 /* controls time window size in P_AVERAGE_ENERGY scheme */ +#define AVRG_E_FACTOR 0.9 /* 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 0 /* set to 1 to add pseudo-3d shading effect */ -#define SHADE_SCALE_2D 0.25 /* lower value increases sensitivity of shading */ +#define SHADE_2D 1 /* set to 1 to add pseudo-3d shading effect */ +#define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */ #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 */ @@ -225,10 +228,10 @@ #define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ -#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 0.2 /* scale of color scheme bar for 2nd part */ -#define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ +#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 CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */ #define CIRC_COLORBAR_B 0 /* set to 1 to draw circular color scheme */ @@ -594,7 +597,7 @@ 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; + 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 *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]; @@ -829,17 +832,17 @@ void animation() /* add oscillating waves */ for (source = 0; source < N_SOURCES; source++) { - angle = APOLY*PID + (double)source*DPI/(double)NPOLY; - wave_source_x[source] = LAMBDA*cos(angle) + 2.0*MU*cos(angle - 0.5*PID); - wave_source_y[source] = LAMBDA*sin(angle) + 2.0*MU*sin(angle - 0.5*PID); - source_shift[source] = source*OSCILLATING_SOURCE_PERIOD/N_SOURCES; - if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == source_shift[source])) + angle = ((double)source)*DPI/(double)NPOLY + APOLY*PID - PID; + wave_source_x[source] = 1.0*LAMBDA*cos(angle); + wave_source_y[source] = 1.0*LAMBDA*sin(angle); + source_shift[source] = 0; +// source_shift[source] = OSCILLATING_SOURCE_PERIOD*source/N_SOURCES; + if ((ADD_OSCILLATING_SOURCE)&&(i%(OSCILLATING_SOURCE_PERIOD) == source_shift[source])&&(i