From 623d353390a4d02d44d9ecf20b3ca0d0a863a533 Mon Sep 17 00:00:00 2001 From: Nils Berglund <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Sat, 17 Aug 2024 12:04:42 +0200 Subject: [PATCH] Add files via upload --- global_3d.c | 10 +- global_ljones.c | 68 +- global_pdes.c | 7 +- heat.c | 3 + lennardjones.c | 471 +++++++--- mangrove.c | 5 + rde.c | 479 ++++++---- schrodinger.c | 3 + sub_lj.c | 2211 +++++++++++++++++++++++++++++++++++++++++++-- sub_rde.c | 424 ++++++++- sub_wave.c | 240 ++++- sub_wave_3d_rde.c | 145 ++- wave_3d.c | 6 +- wave_billiard.c | 122 +-- wave_common.c | 91 +- wave_comparison.c | 5 + wave_energy.c | 5 + wave_sphere.c | 4 + 18 files changed, 3808 insertions(+), 491 deletions(-) diff --git a/global_3d.c b/global_3d.c index 51109a6..3bdbba4 100644 --- a/global_3d.c +++ b/global_3d.c @@ -38,6 +38,8 @@ #define GF_WING 4 /* wing shape */ #define GF_COMPUTE_FROM_BC 5 /* compute force field as gradient of bc_field2 */ #define GF_EARTH 6 /* field depends on altitude on continents */ +#define GF_MARS 7 /* field depends on altitude on Mars */ +#define GF_VENUS 8 /* field depends on altitude on Venus */ /* Choice of water depth for shallow water equation */ @@ -45,10 +47,14 @@ #define SH_CIRCLES 2 /* shallow obstacle specified by CIRCLE_PATTERN */ #define SH_COAST 3 /* depth varying with x-coordinate */ #define SH_COAST_MONOTONE 4 /* depth decreasing with x-coordinate */ +#define SH_SPHERE_CUBE 5 /* cube embedded in sphere */ +#define SH_SPHERE_OCTAHEDRON 6 /* octahedron embedded in sphere */ +#define SH_SPHERE_DODECAHEDRON 7 /* dodecahedron embedded in sphere */ +#define SH_SPHERE_ICOSAHEDRON 8 /* icosahedron embedded in sphere */ /* Type of rotating viewpoint */ -#define VP_HORIZONTAL 0 /* rotate in a horizontal plane */ +#define VP_HORIZONTAL 0 /* rotate in a horizontal plane (constant latitude) */ #define VP_ORBIT 1 /* rotate in a plane containing the origin */ #define VP_ORBIT2 11 /* rotate in a plane specified by max latitude */ #define VP_POLAR 2 /* polar orbit */ @@ -83,7 +89,7 @@ #define PLANET ((B_DOMAIN == D_SPHERE_EARTH)||(B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)||(B_DOMAIN == D_SPHERE_VENUS)||(B_DOMAIN == D_SPHERE_MERCURY)) #define OTHER_PLANET ((B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)||(B_DOMAIN == D_SPHERE_VENUS)||(B_DOMAIN == D_SPHERE_MERCURY)) -#define RDE_PLANET ((ADAPT_STATE_TO_BC)&&(OBSTACLE_GEOMETRY == D_SPHERE_EARTH)) +#define RDE_PLANET ((ADAPT_STATE_TO_BC)&&((OBSTACLE_GEOMETRY == D_SPHERE_EARTH)||(OBSTACLE_GEOMETRY == D_SPHERE_MARS)||(OBSTACLE_GEOMETRY == D_SPHERE_VENUS))) #define NMAXCIRC_SPHERE 100 /* max number of circles on sphere */ #define NMAX_TRACER_PTS 20 /* max number of tracer points recorded per cell */ diff --git a/global_ljones.c b/global_ljones.c index e376dfb..4ed659a 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -11,14 +11,16 @@ #define D_CIRCLES 20 /* several circles */ #define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */ -#define NMAXCIRCLES 100000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ +#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 NMAXCOLLISIONS 200000 /* max number of collisions */ #define NMAXPARTNERS 30 /* max number of partners in molecule */ -#define NMAXPARTNERMOLECULES 10 /* max number of partners of a 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 C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -76,6 +78,7 @@ #define S_DAM_BRICKS 17 /* dam made of several bricks */ #define S_HLINE_HOLE 18 /* horizontal line with a hole in the bottom */ #define S_HLINE_HOLE_SPOKES 181 /* horizontal line with a hole in the bottom and extra spokes */ +#define S_HLINE_HOLE_SLOPED 182 /* slanted lines with a hole */ #define S_EXT_CIRCLE_RECT 19 /* particles outside a circle and a rectangle */ #define S_BIN_OPENING 20 /* bin containing particles opening at deactivation time */ #define S_BIN_LARGE 201 /* larger bin */ @@ -88,6 +91,10 @@ #define S_CYLINDER 27 /* walls at top and bottom, for cylindrical b.c. */ #define S_TREE 28 /* Christmas tree(s) */ #define S_CONE 29 /* cone */ +#define S_CONVEYOR_BELT 30 /* conveyor belt */ +#define S_TWO_CONVEYOR_BELTS 31 /* two angled conveyor belts */ +#define S_PERIODIC_CONVEYORS 32 /* one wrapping belt, and one short horizontal belt */ +#define S_TEST_CONVEYORS 321 /* test */ /* particle interaction */ @@ -108,6 +115,10 @@ #define I_COULOMB_IMAGINARY 14 /* Coulomb interaction with "imaginary charge" */ #define I_DNA_CHARGED 15 /* Coulomb-type interaction between end points of DNA nucleotides */ #define I_DNA_CHARGED_B 151 /* stronger Coulomb-type interaction between end points of DNA nucleotides */ +#define I_SEGMENT 16 /* harmonic interaction between segments */ +#define I_SEGMENT_CHARGED 161 /* harmonic interaction between segments and Coulomb interaction between ends*/ +#define I_POLYGON 17 /* harmonic interaction between regular polygons */ +#define I_POLYGON_ALIGN 171 /* harmonic interaction between polygons with an aligning torque */ /* Boundary conditions */ @@ -208,6 +219,9 @@ #define CHEM_DNA_BASE_SPLIT 254 /* aggregation/splitting of DNA molecules when base pairs don't match */ #define CHEM_DNA_ENZYME 255 /* aggregation/splitting of DNA molecules in presence of enzymes */ #define CHEM_DNA_ENZYME_REPAIR 256 /* aggregation/splitting of DNA molecules in presence of enzymes and additional repairing of bad connections */ +#define CHEM_POLYGON_AGGREGATION 26 /* aggregation of polygons */ +#define CHEM_POLYGON_CLUSTER 261 /* clustering of polygons into new clusters */ +#define CHEM_POLYGON_ONECLUSTER 262 /* clustering of polygons, with only one cluster allowed */ /* Initial conditions for chemical reactions */ @@ -254,12 +268,16 @@ #define P_NUMBER 12 /* colors depend on particle number */ #define P_EMEAN 13 /* averaged kinetic energy (with exponential damping) */ #define P_LOG_EMEAN 131 /* log of averaged kinetic energy (with exponential damping) */ +#define P_EMEAN_DENSITY 132 /* averaged kinetic energy divided by the cluster size */ #define P_DIRECT_EMEAN 14 /* averaged version of P_DIRECT_ENERGY */ #define P_NOPARTICLE 15 /* particles are not drawn (only the links between them) */ #define P_NPARTNERS 16 /* number of partners */ #define P_CHARGE 17 /* colors represent charge */ #define P_MOL_ANGLE 18 /* orientation of molecule defined by partners */ #define P_CLUSTER 19 /* colors depend on connected component */ +#define P_CLUSTER_SIZE 20 /* colors depend on size of connected component */ +#define P_CLUSTER_SELECTED 21 /* colors show which clusters are slected for growth */ +#define P_COLLISION 22 /* colors depend on number of collision/reaction */ /* Rotation schedules */ @@ -275,6 +293,8 @@ #define POLY_STAR 0 /* star-shaped graph (central molecule attracts outer ones) */ #define POLY_ALL 1 /* all-to-all coupling */ +#define POLY_STAR_CHARGED 11 /* star-shaped graph with charged molecules */ +#define POLY_POLYGON 12 /* polygonal shape */ #define POLY_WATER 2 /* star-shaped with a 120° separation between anions */ #define POLY_SOAP 3 /* polymers with all-to-all coupling and polar end */ #define POLY_SOAP_B 4 /* polymers with pairwise coupling and polar end */ @@ -287,6 +307,9 @@ #define POLY_DNA_ALT 71 /* simplified model for DNA with different short ends */ #define POLY_DNA_DOUBLE 72 /* simplified model for DNA with double ends for rigidity */ #define POLY_DNA_FLEX 73 /* simplified model for DNA with less backbone rigidity (beta) */ +#define POLY_KITE 8 /* kite for kites and darts quasicrystal */ +#define POLY_DART 81 /* dart for kites and darts quasicrystal */ +#define POLY_SEG_POLYGON 9 /* polygon of segments */ /* Background color schemes */ @@ -365,18 +388,42 @@ typedef struct int partner[NMAXPARTNERS]; /* partner particles for option PAIR_PARTICLES */ short int npartners; /* number of partner particles */ double partner_eqd[NMAXPARTNERS]; /* equilibrium distances between partners */ + double partner_eqa[NMAXPARTNERS]; /* equilibrium angle between partners */ int p0, p1; /* numbers of two first partners (for P_MOL_ANGLE color scheme) */ // short int mol_angle; /* for color scheme P_MOL_ANGLE */ int cluster; /* number of cluster */ + int cluster_color; /* color of cluster */ + int cluster_size; /* size of cluster */ int molecule; /* number of molecule */ short int tested, cactive; /* for cluster search */ short int coulomb; /* has value 1 if DNA-Coulomb interaction is attractive */ short int added; /* has value 1 if particle has been added */ short int reactive; /* has value 1 if particle can react */ short int paired; /* has value 1 if belongs to base-paired molecule */ + short int flip; /* keeps track of which particles in a cluster are flipped by PI */ int partner_molecule; /* number of partner molecule */ + int collision; /* number of collision */ } t_particle; +typedef struct +{ + short int active; /* has value 1 if cluster is active */ + short int thermostat; /* has value 1 if cluster is coupled to thermostat */ + short int selected; /* has value 1 if cluster is selected to be able to grow */ + double xg, yg; /* center of gravity */ + double vx, vy; /* velocity of center of gravity */ + double angle; /* orientation of cluster */ + double omega; /* angular velocity of cluster */ + double mass, mass_inv; /* mass of cluster and its inverse */ + double inertia_moment, inertia_moment_inv; /* moment of inertia */ + double fx, fy, torque; /* force and torque */ + double energy, emean; /* energy and averaged energy */ + double dirmean; /* time-averaged direction */ + int nparticles; /* number of particles in cluster */ + int particle[NMAXPARTINCLUSTER]; /* list of particles in cluster */ +// int angle_ref; /* reference particle for orientation */ +} t_cluster; + typedef struct { int number; /* total number of particles in cell */ @@ -433,6 +480,8 @@ typedef struct double pressure; /* pressure acting on segement */ double avrg_pressure; /* time-averaged pressure */ short int inactivate; /* set to 1 for segment to become inactive at time SEGMENT_DEACTIVATION_TIME */ + short int conveyor; /* set to 1 for segment to exert lateral force */ + double conveyor_speed; /* speed of conveyor belt */ } t_segment; typedef struct @@ -468,13 +517,24 @@ typedef struct typedef struct { int nparticles; /* number of particles */ - int particle[NPARTNERS+1]; /* list of particles */ + int particle[2*NPARTNERS+1]; /* list of particles */ int npartners; /* number of partner molecules */ int partner[NMAXPARTNERMOLECULES]; /* list of partner molecules */ int connection_type[NMAXPARTNERMOLECULES]; /* types of particles in connection */ short int added; /* has value 1 if molecule has been added */ } t_molecule; +typedef struct +{ + double x1, y1, x2, y2; /* positions of extremities */ + double width; /* width of belt */ + double speed; /* speed of conveyor belt */ + double position; /* position of belt (needed for display of rotating parts) */ + double length; /* distance between (x1,x2) and (y1,y2) */ + double angle; /* angle of (x1,x2) - (y1,y2) */ + double tx, ty; /* coordinates of tangent vector */ +} t_belt; + typedef struct { int nactive; /* number of active particles */ @@ -496,6 +556,6 @@ typedef struct -int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0; +int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0; FILE *lj_log; diff --git a/global_pdes.c b/global_pdes.c index 1286202..f94121d 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -103,6 +103,9 @@ #define D_RITCHEY_CHRETIEN_SPHERICAL 75 /* Ritchey-Chrétien telescope with spherical mirrors */ #define D_RITCHEY_CHRETIEN_HYPERBOLIC 751 /* Ritchey-Chrétien telescope with hyperbolic mirrors */ #define D_GRADIENT_INDEX_LENS 76 /* gradient index lens (only affects draw_billiard) */ +#define D_IMAGE 77 /* Taken from image file */ +#define D_MAGNETRON 78 /* simplified magnetron */ +#define D_MAGNETRON_CATHODE 781 /* simplified magnetron with central cathode */ /* for wave_sphere.c */ @@ -141,7 +144,9 @@ #define C_RINGS 20 /* obstacles arranged in concentric rings */ #define C_RINGS_T 201 /* obstacles arranged in concentric rings, triangular lattice */ -#define C_RINGS_SPIRAL 202 /* obstacles arranged on a "subflower" spiral, similar to C_GOLDEN_SPIRAL */ +#define C_RINGS_SPIRAL 202 /* obstacles arranged on a "sunflower" spiral, similar to C_GOLDEN_SPIRAL */ +#define C_RINGS_POISSONDISC 203 /* obstacles arranged in a Poisson disc pattern */ +#define C_RINGS_LOGSPIRAL 204 /* logarithmic spirals */ #define C_HEX_BOTTOM 101 /* hex/triangular lattice in lower half */ #define C_HEX_BOTTOM2 102 /* smaller hex/triangular lattice in lower half */ diff --git a/heat.c b/heat.c index 383b4ce..db20b5d 100644 --- a/heat.c +++ b/heat.c @@ -70,6 +70,7 @@ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ +#define PDISC_FACTOR 3.25 /* controls density of Poisson disc process (default: 3.25) */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 1.1 /* parameter controlling the dimensions of domain */ @@ -158,6 +159,7 @@ // #define SLOPE 0.1 /* sensitivity of color on wave amplitude */ #define SLOPE 0.2 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ +#define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #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 */ @@ -219,6 +221,7 @@ #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 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 */ #define MESSAGE_LDASH 14 /* length of dash for Morse code message */ diff --git a/lennardjones.c b/lennardjones.c index 9a512f2..4a9e907 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -36,8 +36,8 @@ #include #include -#define MOVIE 1 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define 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 while saving frames */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -58,10 +58,10 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.95 -#define INITXMAX 1.95 /* x interval for initial condition */ -#define INITYMIN -1.1 -#define INITYMAX 1.1 /* y interval for initial condition */ +#define INITXMIN -1.8 +#define INITXMAX 1.8 /* x interval for initial condition */ +#define INITYMIN 0.9 +#define INITYMAX 2.2 /* y interval for initial condition */ #define THERMOXMIN -1.25 #define THERMOXMAX 1.25 /* x interval for initial condition */ @@ -78,7 +78,7 @@ #define BCXMIN -2.0 #define BCXMAX 2.0 /* x interval for boundary condition */ #define BCYMIN -1.125 -#define BCYMAX 1.125 /* y interval for boundary condition */ +#define BCYMAX 2.4 /* y interval for boundary condition */ #define OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ @@ -91,14 +91,16 @@ #define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */ #define OBSTACLE_PATTERN 71 /* pattern of obstacles, see list in global_ljones.c */ -#define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */ -#define SEGMENT_PATTERN 29 /* pattern of repelling segments, see list in global_ljones.c */ +#define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */ +#define SEGMENT_PATTERN 182 /* 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 */ #define NOZZLE_SHAPE_B 6 /* shape of nozzle for second rocket, see list in global_ljones.c */ +#define BELT_SPEED1 20.0 /* speed of first conveyor belt */ +#define BELT_SPEED2 -10.0 /* speed of second conveyor belt */ -#define TWO_TYPES 1 /* set to 1 to have two types of particles */ +#define TWO_TYPES 0 /* set to 1 to have two types of particles */ #define TYPE_PROPORTION 0.5 /* proportion of particles of first type */ #define TWOTYPE_CONFIG 0 /* choice of types, see TTC_ list in global_ljones.c */ #define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ @@ -106,31 +108,31 @@ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ #define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ -#define INTERACTION 12 /* particle interaction, see list in global_ljones.c */ -#define INTERACTION_B 12 /* particle interaction for second type of particle, see list in global_ljones.c */ -#define SPIN_INTER_FREQUENCY 2.0 /* angular frequency of spin-spin interaction */ -#define SPIN_INTER_FREQUENCY_B 2.0 /* angular frequency of spin-spin interaction for second particle type */ -#define MOL_ANGLE_FACTOR 1.0 /* rotation angle for P_MOL_ANGLE color scheme */ +#define INTERACTION 171 /* 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 6.0 /* angular frequency of spin-spin interaction */ +#define SPIN_INTER_FREQUENCY_B 6.0 /* angular frequency of spin-spin interaction for second particle type */ +#define MOL_ANGLE_FACTOR 6.0 /* rotation angle for P_MOL_ANGLE color scheme */ #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 9.6 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 2.5 /* 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.0087 /* parameter controlling radius of particles */ -#define MU_B 0.012 /* parameter controlling radius of particles of second type */ -#define NPOLY 40 /* number of sides of polygon */ -#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define MU 0.025 /* parameter controlling radius of particles */ +#define MU_B 0.03 /* parameter controlling radius of particles of second type */ +#define NPOLY 6 /* number of sides of polygon */ +#define APOLY 0.075 /* 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 */ #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 40 /* number of grid point for grid of disks */ -#define NGRIDY 20 /* number of grid point for grid of disks */ +#define NGRIDX 36 /* number of grid point for grid of disks */ +#define NGRIDY 36 /* number of grid point for grid of disks */ #define EHRENFEST_RADIUS 0.9 /* radius of container for Ehrenfest urn configuration */ #define EHRENFEST_WIDTH 0.035 /* width of tube for Ehrenfest urn configuration */ #define TWO_CIRCLES_RADIUS_RATIO 0.8 /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */ @@ -146,11 +148,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2500 /* number of frames of movie */ -// #define NSTEPS 200 /* number of frames of movie */ -#define NVID 50 /* number of iterations between images displayed on screen */ +#define NSTEPS 2700 /* number of frames of movie */ +#define NVID 200 /* number of iterations between images displayed on screen */ #define NSEG 25 /* number of segments of boundary of circles */ -#define INITIAL_TIME 5 /* time after which to start saving frames */ +#define INITIAL_TIME 200 /* 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 */ @@ -161,45 +162,49 @@ #define SLEEP1 1 /* initial sleeping time */ #define SLEEP2 1 /* final sleeping time */ #define MID_FRAMES 100 /* number of still frames between parts of two-part movie */ -// #define END_FRAMES 250 /* number of still frames at end of movie */ #define END_FRAMES 100 /* number of still frames at end of movie */ /* Boundary conditions, see list in global_ljones.c */ -#define BOUNDARY_COND 3 +#define BOUNDARY_COND 1 /* Plot type, see list in global_ljones.c */ -#define PLOT 5 -#define PLOT_B 13 /* plot type for second movie */ +#define PLOT 11 +#define PLOT_B 14 /* plot type for second movie */ /* Background color depending on particle properties */ -#define COLOR_BACKGROUND 1 /* set to 1 to color background */ +#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 DRAW_BONDS 1 /* set to 1 to draw bonds between neighbours */ +#define DRAW_BONDS 0 /* set to 1 to draw bonds between neighbours */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ #define FILL_TRIANGLES 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 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 200 /* 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.995 /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */ +#define DRATIO 0.999 /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */ /* Color schemes */ #define COLOR_PALETTE 10 /* Color palette, see list in global_ljones.c */ #define COLOR_PALETTE_EKIN 10 /* Color palette for kinetic energy */ -#define COLOR_PALETTE_ANGLE 10 /* Color palette for angle representation */ -#define COLOR_PALETTE_DIRECTION 0 /* Color palette for direction representation */ +#define COLOR_PALETTE_ANGLE 0 /* Color palette for angle representation */ +#define COLOR_PALETTE_DIRECTION 17 /* Color palette for direction representation */ #define COLOR_PALETTE_INITIAL_POS 10 /* Color palette for initial position representation */ #define COLOR_PALETTE_DIFFNEIGH 10 /* Color palette for different neighbours representation */ #define COLOR_PALETTE_PRESSURE 11 /* Color palette for different neighbours representation */ #define COLOR_PALETTE_CHARGE 18 /* Color palette for charge representation */ -#define COLOR_PALETTE_CLUSTER 11 /* Color palette for cluster representation */ +#define COLOR_PALETTE_CLUSTER 14 /* Color palette for cluster representation */ +#define COLOR_PALETTE_CLUSTER_SIZE 13 /* Color palette for cluster size representation */ +#define COLOR_PALETTE_CLUSTER_SELECTED 11 /* Color palette for selected cluster representation */ +#define COLOR_HUE_CLUSTER_SELECTED 90.0 /* Color hue for selected cluster */ +#define COLOR_HUE_CLUSTER_NOT_SELECTED 220.0 /* Color hue for selected cluster */ #define BLACK 1 /* background */ @@ -234,12 +239,11 @@ #define ENERGY_HUE_MAX 50.0 /* color of saturated particle */ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ -// #define PARTICLE_EMAX 50000.0 /* energy of particle with hottest color */ #define PARTICLE_EMIN 10.0 /* energy of particle with coolest color */ -#define PARTICLE_EMAX 20000.0 /* energy of particle with hottest color */ -#define HUE_TYPE0 300.0 /* hue of particles of type 0 */ -#define HUE_TYPE1 00.0 /* hue of particles of type 1 */ -#define HUE_TYPE2 340.0 /* hue of particles of type 2 */ +#define PARTICLE_EMAX 100.0 /* energy of particle with hottest color */ +#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 */ #define HUE_TYPE3 260.0 /* hue of particles of type 3 */ #define HUE_TYPE4 200.0 /* hue of particles of type 4 */ #define HUE_TYPE5 60.0 /* hue of particles of type 5 */ @@ -249,39 +253,40 @@ #define RANDOM_RADIUS 0 /* set to 1 for random circle radius */ #define DT_PARTICLE 3.0e-6 /* time step for particle displacement */ -#define KREPEL 15.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 5.0 /* Lennard-Jones equilibrium distance */ -#define EQUILIBRIUM_DIST_B 5.0 /* Lennard-Jones equilibrium distance for second type of particle */ +#define KREPEL 50.0 /* constant in repelling force between particles */ +#define EQUILIBRIUM_DIST 2.3 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST_B 2.5 /* Lennard-Jones equilibrium distance for second type of particle */ #define REPEL_RADIUS 25.0 /* radius in which repelling force acts (in units of particle radius) */ #define DAMPING 500.0 /* damping coefficient of particles */ -#define INITIAL_DAMPING 5000.0 /* damping coefficient of particles during initial phase */ -#define DAMPING_ROT 1.0e6 /* damping coefficient for rotation of particles */ +#define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */ +#define DAMPING_ROT 1000.0 /* damping coefficient for rotation of particles */ #define PARTICLE_MASS 2.0 /* mass of particle of radius MU */ -#define PARTICLE_MASS_B 16.0 /* mass of particle of radius MU_B */ -#define PARTICLE_INERTIA_MOMENT 0.2 /* moment of inertia of particle */ -#define PARTICLE_INERTIA_MOMENT_B 0.02 /* moment of inertia of second type of particle */ +#define 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 50.0 /* initial velocity range */ -#define OMEGA_INITIAL 10.0 /* initial angular velocity range */ +#define OMEGA_INITIAL 50.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 */ +#define COULOMB_LJ_FACTOR 1.0 /* relative intensity of LJ interaction in I_COULOMB_LJ interaction (default: 0.01) */ #define V_INITIAL_TYPE 0 /* type of initial speed distribution (see VI_ in global_ljones.c) */ -#define THERMOSTAT 1 /* set to 1 to switch on thermostat */ +#define THERMOSTAT 0 /* set to 1 to switch on thermostat */ #define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ #define SIGMA 5.0 /* noise intensity in thermostat */ -#define BETA 0.00007 /* initial inverse temperature */ +#define BETA 0.004 /* initial inverse temperature */ #define MU_XI 0.005 /* friction constant in thermostat */ #define KSPRING_BOUNDARY 2.0e11 /* confining harmonic potential outside simulation region */ -#define KSPRING_OBSTACLE 2.0e11 /* harmonic potential of obstacles */ -#define NBH_DIST_FACTOR 5.0 /* radius in which to count neighbours */ -#define GRAVITY 0.0 /* gravity acting on all particles */ +#define KSPRING_OBSTACLE 1.0e9 /* harmonic potential of obstacles */ +#define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ +#define GRAVITY 4000.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 */ #define GRAVITY_SCHEDULE 1 /* type of gravity schedule, see list in global_ljones.c */ #define GRAVITY_FACTOR 10.0 /* factor by which to increase gravity */ -#define GRAVITY_INITIAL_TIME 200 /* time at start of simulation with constant gravity */ +#define GRAVITY_INITIAL_TIME 250 /* time at start of simulation with constant gravity */ #define GRAVITY_RESTORE_TIME 500 /* time at end of simulation with gravity restored to initial value */ #define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */ #define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */ @@ -291,9 +296,8 @@ #define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ #define BFIELD 2.666666667 /* value of magnetic field */ #define CHARGE 1.0 /* charge of particles of first type */ -#define CHARGE_B -1.5 /* charge of particles of second type */ +#define CHARGE_B 1.0 /* charge of particles of second type */ #define INCREASE_E 0 /* set to 1 to increase electric field */ -// #define EFIELD_FACTOR 2500000.0 /* factor by which 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 */ #define BFIELD_FACTOR 20000.0 /* factor by which to increase magnetic field */ @@ -301,25 +305,24 @@ #define OBSTACLE_CHARGE 3.0 /* charge of obstacles */ #define KCOULOMB_OBSTACLE 1000.0 /* Coulomb force constant for charged obstacles */ -#define ROTATION 0 /* set to 1 to include rotation of particles */ +#define ROTATION 1 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ #define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */ -#define KTORQUE -100.0 /* force constant in angular dynamics */ -#define KTORQUE_BOUNDARY 1.0e6 /* constant in torque from the boundary */ +#define KTORQUE 5.0e3 /* force constant in angular dynamics */ +#define KTORQUE_BOUNDARY 5.0e5 /* constant in torque from the boundary */ #define KTORQUE_B 10.0 /* force constant in angular dynamics */ -#define KTORQUE_DIFF -150.0 /* force constant in angular dynamics for different particles */ +#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 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 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 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 */ -#define INCREASE_BETA 1 /* set to 1 to increase BETA during simulation */ +#define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */ #define BETA_SCHEDULE 3 /* type of temperature schedule, see TS_* in global_ljones */ -#define BETA_FACTOR 10.0 /* factor by which to change BETA during simulation */ -// #define BETA_FACTOR 0.08 /* factor by which to change BETA during simulation */ +#define BETA_FACTOR 0.06 /* factor by which to change BETA during simulation */ #define TS_SLOPE 8.5 /* controls speed of change of BETA for TS_TANH schedule (default 1.0) */ #define N_TOSCILLATIONS 1.0 /* number of temperature oscillations in BETA schedule */ #define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */ @@ -380,7 +383,7 @@ #define SMOOTH_ROTATION 1 /* set to 1 to update segments at each time step (rather than at each movie frame) */ #define ROTATION_SCHEDULE 0 /* time-dependence of rotation angle, see ROT_* in global_ljones.c */ #define PERIOD_ROTATE_BOUNDARY 1000 /* period of rotating boundary */ -#define ROTATE_INITIAL_TIME 300 /* initial time without rotation */ +#define ROTATE_INITIAL_TIME 150 /* initial time without rotation */ #define ROTATE_FINAL_TIME 300 /* final time without rotation */ #define ROTATE_CHANGE_TIME 0.5 /* relative duration of acceleration/deceleration phases */ #define OMEGAMAX -2.0*PI /* maximal rotation speed */ @@ -388,7 +391,7 @@ #define MOVE_BOUNDARY 0 /* set to 1 to move repelling segments, due to force from particles */ #define SEGMENTS_MASS 40.0 /* mass of collection of segments */ #define DEACTIVATE_SEGMENT 1 /* set to 1 to deactivate last segment after a certain time */ -#define SEGMENT_DEACTIVATION_TIME 20 /* time at which to deactivate last segment */ +#define SEGMENT_DEACTIVATION_TIME 200 /* time at which to deactivate last segment */ #define RELEASE_ROCKET_AT_DEACTIVATION 0 /* set to 1 to limit segments velocity before segment release */ #define SEGMENTS_X0 1.5 /* initial position of segments */ #define SEGMENTS_Y0 0.0 /* initial position of segments */ @@ -405,6 +408,7 @@ #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 */ #define KSPRING_GROUPS 5.0e11 /* harmonic potential between segment groups */ +#define KSPRING_BELT 1.0e4 /* spring constant from belt */ #define GROUP_WIDTH 0.05 /* interaction width of groups */ #define GROUP_G_REPEL 0 /* set to 1 to add repulsion between centers of mass of groups */ #define GROUP_G_REPEL_RADIUS 1.2 /* radius within which centers of mass of groups repel each other */ @@ -418,20 +422,29 @@ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ #define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */ -#define RD_REACTION 256 /* type of reaction, see list in global_ljones.c */ -#define RD_TYPES 6 /* number of types in reaction-diffusion equation */ -#define RD_INITIAL_COND 10 /* initial condition of particles */ -#define REACTION_DIST 4.0 /* maximal distance for reaction to occur */ +#define RD_REACTION 262 /* type of reaction, see list in global_ljones.c */ +#define RD_TYPES 1 /* number of types in reaction-diffusion equation */ +#define RD_INITIAL_COND 0 /* initial condition of particles */ +#define REACTION_DIST 1.2 /* maximal distance for reaction to occur */ #define REACTION_PROB 1.0 /* probability controlling reaction term */ #define DISSOCIATION_PROB 0.0 /* probability controlling dissociation reaction */ -#define KILLING_PROB 0.0015 /* probability of enzymes being killed */ +#define KILLING_PROB 0.0015 /* probability of enzymes being killed */ +#define DELTAMAX 0.1 /* max orientation difference for pairing polygons */ #define CENTER_COLLIDED_PARTICLES 0 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */ #define EXOTHERMIC 0 /* set to 1 to make reaction exo/endothermic */ #define DELTA_EKIN 2000.0 /* change of kinetic energy in reaction */ -#define COLLISION_TIME 25 /* time during which collisions are shown */ -#define DELTAVMAX 1000.0 /* maximal deltav allowed for pairing molecules */ -#define AGREGMAX 11 /* maximal number of partners for CHEM_AGGREGATION reaction */ -#define AGREG_DECOUPLE 12 /* minimal number of partners to decouple from thermostat */ +#define COLLISION_TIME 25 /* time during which collisions are shown */ +#define COLLISION_RADIUS 2.0 /* radius of discs showing collisions, in units of MU */ +#define DELTAVMAX 200.0 /* maximal deltav allowed for pairing molecules */ +#define AGREGMAX 6 /* maximal number of partners for CHEM_AGGREGATION reaction */ +#define AGREG_DECOUPLE 12 /* minimal number of partners to decouple from thermostat */ +#define CLUSTER_PARTICLES 0 /* set to 1 for particles to form rigid clusters */ +#define CLUSTER_MAXSIZE 1000 /* max size of clusters */ +#define SMALL_CLUSTER_MAXSIZE 10 /* size limitation on smaller cluster */ +#define SMALL_NP_MAXSIZE 2 /* limitation on number of partners of particle in smaller cluster */ +#define NOTSELECTED_CLUSTER_MAXSIZE 0 /* limit on size of clusters that can merge with non-selected cluster */ +#define REPAIR_CLUSTERS 0 /* set to 1 to repair alignment in clusters */ +#define REPAIR_MIN_DIST 0.75 /* relative distance below which overlapping polygons are inactivated */ #define CHANGE_RADIUS 0 /* set to 1 to change particle radius during simulation */ #define MU_RATIO 0.666666667 /* ratio by which to increase radius */ @@ -458,35 +471,36 @@ #define PROP_MIN 0.1 /* min proportion of type 1 particles */ #define PROP_MAX 0.9 /* max proportion of type 1 particles */ -#define PAIR_PARTICLES 1 /* set to 1 to form particles pairs */ +#define PAIR_PARTICLES 0 /* set to 1 to form particles pairs */ #define RANDOMIZE_ANGLE 0 /* set to 1 for random orientation */ #define DEACIVATE_CLOSE_PAIRS 1 /* set to 1 to test for closeness to other particles */ #define PAIR_SAFETY_FACTOR 1.2 /* distance to deactivate divided by sum of radii */ #define THIRD_TYPE_PROPORTION 1.0 /* proportion of third type pairings, for certain pairing types */ -#define KSPRING_PAIRS 1.0e10 /* spring constant for pair interaction */ +#define KSPRING_PAIRS 5.0e10 /* spring constant for pair interaction */ #define KTORQUE_PAIRS 1.0e10 /* constant for angular coupling in pair interaction */ #define KTORQUE_PAIR_ANGLE 0.0 /* constant for coupling between orientation in pairs */ -#define NPARTNERS 8 /* number of partners of particles */ -#define NARMS 1 /* number of "arms" for certain paring types */ -#define PAIRING_TYPE 42 /* type of pairing, see POLY_ in global_ljones.c */ +#define NPARTNERS 4 /* number of partners of particles - for DNA, set NPARTNERS_DNA */ +#define NPARTNERS_DNA 8 /* number of partners of particles, case of DNA, should be at least 8 */ +#define NARMS 5 /* number of "arms" for certain paring types */ +#define PAIRING_TYPE 9 /* type of pairing, see POLY_ in global_ljones.c */ #define PARTNER_ANGLE 104.45 /* angle (in degrees) between ions for POLY_WATER case */ -#define PAIR_DRATIO 0.9 /* ratio between equilibrium distance and radius (default: 1.0) */ -#define MU_C 0.0125 /* radius of partner particle */ -#define PARTICLE_MASS_C 8.0 /* mass or partner particle */ -#define CHARGE_C -1.0 /* charge of partner particle */ -#define CLUSTER_COLOR_FACTOR 400 /* factor for initialization of cluster colors */ +#define PAIR_DRATIO 1.1 /* ratio between equilibrium distance and radius (default: 1.0) */ +#define MU_C 0.035 /* radius of partner particle */ +#define PARTICLE_MASS_C 2.0 /* mass or partner particle */ +#define CHARGE_C 1.0 /* charge of partner particle */ +#define CLUSTER_COLOR_FACTOR 40 /* factor for initialization of cluster colors */ #define ALTERNATE_POLY_CHARGE 1 /* set to 1 for alternating charges in molecule */ #define SECONDARY_PAIRING 0 /* set to 1 to pair with secondary partners, experimental */ #define DNA_RIGIDITY 0.5 /* controls rigidity for POLY_DNA_DOUBLE pairs, default = 1 */ #define PAIR_TYPEB_PARTICLES 1 /* set to 1 to pair particle of type 1 */ -#define NPARTNERS_B 2 /* number of partners of particles */ +#define NPARTNERS_B 5 /* number of partners of particles */ #define NARMS_B 1 /* number of "arms" for certain paring types */ -#define PAIRING_TYPE_B 2 /* type of pairing, see POLY_ in global_ljones.c */ -#define MU_D 0.008 /* radius of partner particle */ -#define PARTICLE_MASS_D 1.0 /* mass or partner particle */ -#define CHARGE_D 0.75 /* charge of partner particle */ +#define PAIRING_TYPE_B 81 /* type of pairing, see POLY_ in global_ljones.c */ +#define MU_D 0.035 /* radius of partner particle */ +#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 */ @@ -500,8 +514,8 @@ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 80 /* size of hashgrid in x direction */ -#define HASHY 40 /* size of hashgrid in y direction */ +#define HASHX 100 /* size of hashgrid in x direction */ +#define HASHY 88 /* 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 */ @@ -515,11 +529,12 @@ #define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO)) #define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE)) #define TWO_OBSTACLES ((SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)||(SEGMENT_PATTERN == S_TWO_ROCKETS)) -#define COMPUTE_EMEAN ((PLOT == P_EMEAN)||(PLOT_B == P_EMEAN)||(PLOT == P_LOG_EMEAN)||(PLOT_B == P_LOG_EMEAN)||(PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN)) +#define COMPUTE_EMEAN ((PLOT == P_EMEAN)||(PLOT_B == P_EMEAN)||(PLOT == P_LOG_EMEAN)||(PLOT_B == P_LOG_EMEAN)||(PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN)||(PLOT == P_EMEAN_DENSITY)||(PLOT_B == P_EMEAN_DENSITY)) #define COMPUTE_DIRMEAN ((PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN)) #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)))) +#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))) double xshift = 0.0; /* x shift of shown window */ double xspeed = 0.0; /* x speed of obstacle */ @@ -960,6 +975,8 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA if (initial_phase) damping = INITIAL_DAMPING; else damping = DAMPING; +// printf("Evolving particles\n"); + #pragma omp parallel for private(j,xi,totalenergy,a,move) for (j=0; j direction + PI) dmean -= DPI; particle[j].dirmean = DRATIO*dmean + (1.0-DRATIO)*direction; @@ -1092,6 +1111,135 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA return(totalenergy); } +double evolve_clusters(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], + t_hashgrid hashgrid[HASHX*HASHY], + double cqx[NMAXCIRCLES], double cqy[NMAXCIRCLES], double cqangle[NMAXCIRCLES], + double cpx[NMAXCIRCLES], double cpy[NMAXCIRCLES], double cpangle[NMAXCIRCLES], + double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase, int verbose) +{ + double a, totalenergy = 0.0, damping, direction, dmean, newx, newy, newangle, deltax, deltay, deltaangle; + static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; + int j, move, ncoup; + + if (initial_phase) damping = INITIAL_DAMPING; + else damping = DAMPING; + + #pragma omp parallel for private(j,xi,totalenergy,a,move) + for (j=0; j direction + PI) dmean -= DPI; + cluster[j].dirmean = DRATIO*dmean + (1.0-DRATIO)*direction; + if (cluster[j].dirmean < 0.0) cluster[j].dirmean += DPI; + else if (cluster[j].dirmean > DPI) cluster[j].dirmean -= DPI; + } + + if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(cluster[j].thermostat)) + cluster[j].energy += cpangle[j]*cpangle[j]*cluster[j].inertia_moment_inv; + + cqx[j] = cluster[j].xg + 0.5*DT_PARTICLE*cpx[j]*cluster[j].mass_inv; + cqy[j] = cluster[j].yg + 0.5*DT_PARTICLE*cpy[j]*cluster[j].mass_inv; + cqangle[j] = cluster[j].angle + 0.5*DT_PARTICLE*cpangle[j]*cluster[j].inertia_moment_inv; + + if ((THERMOSTAT_ON)&&(cluster[j].thermostat)) + { + cpx[j] *= exp(- 0.5*DT_PARTICLE*xi); + cpy[j] *= exp(- 0.5*DT_PARTICLE*xi); + } + if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(cluster[j].thermostat)) + cpangle[j] *= exp(- 0.5*DT_PARTICLE*xi); + } + + /* compute kinetic energy */ +// *nactive = 0; + ncoup = 1; + for (j=0; j 0) +// { +// compute_relative_positions(particle, hashgrid); +// update_hashgrid(particle, hashgrid, 0); /* REDUNDANT ? */ +// } + } + +// sleep(1); + + *ncoupled = ncoup; + return(totalenergy); +} void evolve_lid(double fboundary) { @@ -1394,16 +1542,18 @@ 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; + angle, theta, sum, alpha; double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds; - int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, + 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; + group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt; int *particle_numbers; static int imin, imax; static short int first = 1; t_particle *particle; + t_cluster *cluster; t_obstacle *obstacle; t_segment *segment; t_group_segments *segment_group; @@ -1413,6 +1563,7 @@ void animation() t_hashgrid *hashgrid; t_molecule *molecule; t_lj_parameters params; + t_belt *conveyor_belt; char message[100]; ratioc = 1.0 - ratio; @@ -1426,14 +1577,19 @@ void animation() params.gravity = GRAVITY; params.radius = MU; - particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ + particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ + + if (CLUSTER_PARTICLES) + cluster = (t_cluster *)malloc(NMAXCIRCLES*sizeof(t_cluster)); /* particle clusters */ + if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle)); /* circular obstacles */ if (ADD_FIXED_SEGMENTS) { segment = (t_segment *)malloc(NMAXSEGMENTS*sizeof(t_segment)); /* linear obstacles */ segment_group = (t_group_segments *)malloc(NMAXGROUPS*sizeof(t_group_segments)); + conveyor_belt = (t_belt *)malloc(NMAXBELTS*sizeof(t_belt)); } - + if (TRACER_PARTICLE) trajectory = (t_tracer *)malloc(TRAJECTORY_LENGTH*N_TRACER_PARTICLES*sizeof(t_tracer)); @@ -1450,6 +1606,16 @@ void animation() pangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); pressure = (double *)malloc(N_PRESSURES*sizeof(double)); + if (CLUSTER_PARTICLES) + { + cqx = (double *)malloc(NMAXCIRCLES*sizeof(double)); + cqy = (double *)malloc(NMAXCIRCLES*sizeof(double)); + cpx = (double *)malloc(NMAXCIRCLES*sizeof(double)); + cpy = (double *)malloc(NMAXCIRCLES*sizeof(double)); + cqangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); + cpangle = (double *)malloc(NMAXCIRCLES*sizeof(double)); + } + if (REACTION_DIFFUSION) { collisions = (t_collision *)malloc(2*NMAXCOLLISIONS*sizeof(t_collision)); @@ -1465,7 +1631,7 @@ void animation() lj_log = fopen("lj_logfile.txt", "w"); if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle); - if (ADD_FIXED_SEGMENTS) init_segment_config(segment); + if (ADD_FIXED_SEGMENTS) init_segment_config(segment, conveyor_belt); if ((MOVE_SEGMENT_GROUPS)&&(ADD_FIXED_SEGMENTS)) { @@ -1503,6 +1669,8 @@ void animation() params.nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment, molecule); printf("%i active particles\n", params.nactive); + + if (CLUSTER_PARTICLES) init_cluster_config(particle, cluster); // xi = 0.0; @@ -1510,7 +1678,7 @@ void animation() // { // printf("Particle %i at (%.3f, %.3f) of energy %.3f\n", i, particle[i].xc, particle[i].yc, particle[i].energy); // } - sleep(1); +// sleep(1); update_hashgrid(particle, hashgrid, 1); printf("Updated hashgrid\n"); @@ -1695,6 +1863,29 @@ void animation() } } } + /* TEST: average angles for clustered polygons */ +// if ((REACTION_DIFFUSION)&&(RD_REACTION == CHEM_POLYGON_AGGREGATION)) +// { +// np = particle[j].npartners; +// alpha = DPI/(double)NPOLY; +// if (np >= 2) +// { +// angle = particle[j].angle; +// while (angle > alpha) angle -= alpha; +// sum = angle; +// for (p=0; p alpha) angle -= alpha; +// sum += angle; +// } +// sum *= 1.0/(double)(np+1); +// particle[j].angle = sum; +// for (p=0; p= 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) @@ -1755,11 +1965,24 @@ void animation() } if ((MOVE_BOUNDARY)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segments(segment, i); + if (ADD_CONVEYOR_FORCE) for (belt = 0; belt < nbelts; belt++) + conveyor_belt[belt].position += conveyor_belt[belt].speed*DT_PARTICLE; + 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); } /* end of for (n=0; nINITIAL_TIME)&&(SAVE_TIME_SERIES)) { n_total_active = 0; @@ -1881,7 +2104,7 @@ void animation() /* case of reaction-diffusion equation */ if ((i > INITIAL_TIME)&&(REACTION_DIFFUSION)) { - ncollisions = update_types(particle, molecule, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1, &delta_energy); + ncollisions = update_types(particle, molecule, cluster, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1, &delta_energy); if (EXOTHERMIC) params.beta *= 1.0/(1.0 + delta_energy/totalenergy); params.nactive = 0; for (j=0; j ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD)) @@ -1956,8 +2179,8 @@ void animation() count_particle_number(particle, particle_numbers, i - INITIAL_TIME); draw_frame(i, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, - wall, pressure, pleft, pright, particle_numbers, 1, params, particle, - collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group); + wall, pressure, pleft, pright, particle_numbers, 1, params, particle, cluster, + collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt); if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers(); @@ -1988,8 +2211,8 @@ void animation() if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) { draw_frame(i, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, - wall, pressure, pleft, pright, particle_numbers, 0, params, particle, - collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group); + wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, + collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt); glutSwapBuffers(); save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter); counter++; @@ -2029,8 +2252,8 @@ void animation() { blank(); draw_frame(NSTEPS, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, - wall, pressure, pleft, pright, particle_numbers, 0, params, particle, - collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group); + wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, + collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt); } if (DOUBLE_MOVIE) for (i=0; i #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 */ @@ -48,13 +48,9 @@ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1150 /* window height */ -// #define NX 240 /* number of grid points on x axis */ -// #define NY 120 /* number of grid points on y axis */ -// #define NX 960 /* number of grid points on x axis */ -// #define NY 500 /* number of grid points on y axis */ -#define NX 780 /* number of grid points on x axis */ -#define NY 400 /* number of grid points on y axis */ -#define HRES 2 /* factor for high resolution plots */ +#define NX 960 /* number of grid points on x axis */ +#define NY 575 /* number of grid points on y axis */ +#define HRES 1 /* factor for high resolution plots */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ @@ -63,31 +59,32 @@ /* Choice of simulated equation */ -#define RDE_EQUATION 7 /* choice of reaction term, see list in global_3d.c */ +#define RDE_EQUATION 8 /* choice of reaction term, see list in global_3d.c */ #define NFIELDS 3 /* number of fields in reaction-diffusion equation */ #define NLAPLACIANS 0 /* number of fields for which to compute Laplacian */ #define SPHERE 1 /* set to 1 to simulate equation on sphere */ #define DPOLE 1 /* safety distance to poles */ #define DSMOOTH 1 /* size of neighbourhood of poles that are smoothed */ -#define SMOOTHPOLE 0.01 /* smoothing coefficient at poles */ -#define SMOOTHCOTPOLE 0.01 /* smoothing coefficient of cotangent at poles */ +#define SMOOTHPOLE 0.05 /* smoothing coefficient at poles */ +#define SMOOTHCOTPOLE 0.05 /* smoothing coefficient of cotangent at poles */ #define PHISHIFT 0.0 /* shift of phi in 2D plot (in degrees) */ #define SMOOTHBLOCKS 1 /* set to 1 to use blocks of points near the poles */ #define BLOCKDIST 64 /* distance to poles where points are blocked */ #define ZERO_MERIDIAN 190.0 /* choice of zero meridian (will be at left/right boundary of 2d plot) */ +#define POLE_NODRAW 10 /* distance around poles where wave is not drawn */ #define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */ #define ADD_MAGNETIC_FIELD 0 /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */ -#define ADD_FORCE_FIELD 1 /* set to 1 to add a foce field (for compressible Euler equation) */ +#define ADD_FORCE_FIELD 0 /* set to 1 to add a foce field (for compressible Euler equation) */ #define POTENTIAL 7 /* type of potential or vector potential, see list in global_3d.c */ #define FORCE_FIELD 6 /* type of force field, see list in global_3d.c */ #define ADD_CORIOLIS_FORCE 1 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */ -#define VARIABLE_DEPTH 0 /* set to 1 for variable depth in shallow water equation */ -#define SWATER_DEPTH 4 /* variable depth in shallow water equation */ +#define VARIABLE_DEPTH 1 /* set to 1 for variable depth in shallow water equation */ +#define SWATER_DEPTH 5 /* variable depth in shallow water equation */ #define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */ -#define ADAPT_STATE_TO_BC 1 /* to smoothly adapt initial state to obstacles */ +#define ADAPT_STATE_TO_BC 0 /* to smoothly adapt initial state to obstacles */ #define OBSTACLE_GEOMETRY 84 /* geometry of obstacles, as in B_DOMAIN */ #define BC_STIFFNESS 50.0 /* controls region of boundary condition control */ #define CHECK_INTEGRAL 1 /* set to 1 to check integral of first field */ @@ -105,6 +102,7 @@ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 300 /* number of points for Poisson C_RAND_POISSON arrangement */ +#define PDISC_FACTOR 3.25 /* controls density of Poisson disc process (default: 3.25) */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ @@ -120,6 +118,7 @@ #define NGRIDY 8 /* number of grid point for grid of disks */ #define REVERSE_TESLA_VALVE 1 /* set to 1 to orient Tesla valve in blocking configuration */ #define WALL_WIDTH 0.05 /* width of wall separating lenses */ +#define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -137,17 +136,11 @@ /* Physical parameters of wave equation */ -// #define DT 0.0000001 -#define DT 0.00000015 -// #define DT 0.0000002 -// #define DT 0.0000012 -// #define DT 0.000003 -// #define DT 0.0000022 -// #define DT 0.0000012 +#define DT 0.00000025 #define VISCOSITY 0.02 #define POISSON_STIFFNESS 1.0 /* stiffness of Poisson equation solver for incompressible Euler */ -#define DISSIPATION 0.0 +#define DISSIPATION 1.0e-8 #define RPSA 0.75 /* parameter in Rock-Paper-Scissors-type interaction */ #define RPSLZB 0.0 /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */ @@ -163,7 +156,7 @@ #define BZQ 0.0008 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */ #define B_FIELD 10.0 /* magnetic field */ -#define G_FIELD 0.03 /* gravity/Coriolis force */ +#define G_FIELD 0.002 /* gravity/Coriolis force */ #define BC_FIELD 1.0e-5 /* constant in repulsive field from obstacles */ #define AB_RADIUS 0.2 /* radius of region with magnetic field for Aharonov-Bohm effect */ #define K_EULER 50.0 /* constant in stream function integration of Euler equation */ @@ -172,8 +165,8 @@ #define SMOOTHEN_VORTICITY 0 /* set to 1 to smoothen vorticity field in Euler equation */ #define SMOOTHEN_VELOCITY 1 /* set to 1 to smoothen velocity field in Euler equation */ -#define SMOOTHEN_PERIOD 10 /* period between smoothenings */ -#define SMOOTH_FACTOR 0.05 /* factor by which to smoothen */ +#define SMOOTHEN_PERIOD 7 /* period between smoothenings */ +#define SMOOTH_FACTOR 0.15 /* factor by which to smoothen */ #define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */ #define OSCILLATING_SOURCE_PERIOD 1 /* period of oscillating source */ @@ -181,7 +174,7 @@ #define ADD_TRACERS 1 /* set to 1 to add tracer particles (for Euler equations) */ #define N_TRACERS 2000 /* number of tracer particles */ -#define TRACERS_STEP 0.005 /* step size in tracer evolution */ +#define TRACERS_STEP 0.1 /* step size in tracer evolution */ #define T_OUT 2.0 /* outside temperature */ #define T_IN 0.0 /* inside temperature */ @@ -215,8 +208,8 @@ #define LAMINAR_FLOW_YPERIOD 1.0 /* period of laminar flow in y direction */ #define PRESSURE_GRADIENT 0.2 /* amplitude of pressure gradient for Euler equation */ -#define SWATER_MIN_HEIGHT 0.5 /* min height of initial condition for shallow water equation */ -#define DEPTH_FACTOR 0.015 /* proportion of min height in variable depth */ +#define SWATER_MIN_HEIGHT 0.025 /* min height of initial condition for shallow water equation */ +#define DEPTH_FACTOR 0.075 /* proportion of min height in variable depth */ #define TANH_FACTOR 1.0 /* steepness of variable depth */ #define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */ @@ -232,9 +225,8 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1900 /* number of frames of movie */ -// #define NSTEPS 500 /* number of frames of movie */ -#define NVID 100 /* number of iterations between images displayed on screen */ +#define NSTEPS 950 /* number of frames of movie */ +#define NVID 85 /* number of iterations between images displayed on screen */ #define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */ #define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */ #define MAX_DT 0.024 /* maximal value of integration step */ @@ -256,24 +248,23 @@ #define PLOT_SPHERE 1 /* draws fields on a sphere */ #define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ -#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ -// #define ROTATE_ANGLE 90.0 /* total angle of rotation during simulation */ +#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ #define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ #define SHADE_2D 0 /* set to 1 to change luminosity according to normal vector */ -#define VIEWPOINT_TRAJ 0 /* type of viewpoint trajectory */ +#define VIEWPOINT_TRAJ 1 /* type of viewpoint trajectory */ #define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */ #define DRAW_PERIODICISED 0 /* set to 1 to repeat wave periodically in x and y directions */ /* Plot type - color scheme */ -#define CPLOT 62 -#define CPLOT_B 64 +#define CPLOT 70 +#define CPLOT_B 74 /* Plot type - height of 3D plot */ -#define ZPLOT 62 /* z coordinate in 3D plot */ -#define ZPLOT_B 64 /* z coordinate in second 3D plot */ +#define ZPLOT 70 /* z coordinate in 3D plot */ +#define ZPLOT_B 71 /* z coordinate in second 3D plot */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */ #define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */ @@ -283,7 +274,6 @@ #define DRAW_OUTSIDE_GRAY 0 /* experimental - draw outside of billiard in gray */ #define ADD_POTENTIAL_TO_Z 0 /* set to 1 to add the external potential to z-coordinate of plot */ #define ADD_POT_CONSTANT 0.35 /* constant added to wave height */ -#define ADD_HEIGHT_CONSTANT -0.5 /* constant added to wave height */ #define DRAW_DEPTH 0 /* set to 1 to draw water depth */ #define DEPTH_SCALE 0.75 /* vertical scaling of depth plot */ #define DEPTH_SHIFT -0.015 /* vertical shift of depth plot */ @@ -317,7 +307,7 @@ /* Color schemes, see list in global_pdes.c */ #define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 17 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 12 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ #define COLOR_OUT_R 1.0 /* color outside domain */ @@ -326,20 +316,22 @@ #define COLOR_SCHEME 3 /* choice of color scheme */ -#define COLOR_PHASE_SHIFT 0.25 /* phase shift of color scheme, in units of Pi */ +#define PHASE_SHIFT -0.25 /* phase shift of color scheme, in units of Pi (formerly COLOR_PHASE_SHIFT) */ #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 100.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSHIFT_AMPLITUDE 0.0 /* additional shift for wave amplitude */ +#define VSCALE_AMPLITUDE 15.0 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define CURL_SCALE 1.0 /* scaling factor for curl representation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ -#define SLOPE_SCHROD_LUM 100.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ +#define SLOPE_SCHROD_LUM 10.0 /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */ #define MIN_SCHROD_LUM 0.1 /* minimal luminosity in color scheme Z_ARGUMENT*/ #define VSCALE_PRESSURE 2.0 /* additional scaling factor for color scheme Z_EULER_PRESSURE */ #define PRESSURE_SHIFT 10.0 /* shift for color scheme Z_EULER_PRESSURE */ #define PRESSURE_LOG_SHIFT -2.5 /* shift for color scheme Z_EULER_PRESSURE */ -#define VSCALE_WATER_HEIGHT 0.4 /* vertical scaling of water height */ +#define VSCALE_WATER_HEIGHT 30.0 /* vertical scaling of water height */ +#define ADD_HEIGHT_CONSTANT -0.025 /* constant added to wave height */ #define SHADE_SCALE_2D 0.25 /* controls "depth" of 2D shading */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -359,10 +351,10 @@ #define VSCALE_DENSITY 30.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ #define VSCALE_VORTICITY 15.0 /* additional scaling factor for color scheme Z_EULERC_VORTICITY */ #define VORTICITY_SHIFT 0.0 /* vertical shift of vorticity */ -#define ZSCALE_SPEED 0.5 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ +#define ZSCALE_SPEED 300.0 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ #define ZSHIFT_SPEED 0.0 /* additional shift of z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ #define ZSCALE_NORMGRADIENT -0.0001 /* vertical scaling for Z_NORM_GRADIENT */ -#define VSCALE_SWATER 250.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ +#define VSCALE_SWATER 50.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ #define NXMAZE 7 /* width of maze */ #define NYMAZE 7 /* height of maze */ @@ -372,7 +364,7 @@ #define MAZE_WIDTH 0.04 /* 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 2.5 /* scale of color scheme bar */ #define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */ @@ -393,7 +385,6 @@ #define INITIAL_WAVELENGTH 0.1 /* wavelength of initial condition */ #define VSCALE_ENERGY 200.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 AVRG_E_FACTOR 0.99 /* controls time window size in P_AVERAGE_ENERGY scheme */ #define OSCILLATION_SCHEDULE 0 /* oscillation schedule, see list in global_pdes.c */ #define AMPLITUDE 0.8 /* amplitude of periodic excitation */ @@ -436,8 +427,8 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_3D representation */ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; -double light[3] = {-0.40824829, -0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {-8.0, -4.0, 4.0}; /* location of observer for REP_PROJ_3D representation */ +double light[3] = {0.40824829, -0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ +double observer[3] = {8.0, 4.0, 2.5}; /* location of observer for REP_PROJ_3D representation */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ /* constants for simulations on planets */ @@ -449,7 +440,7 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO #define DEM_SMOOTH_HEIGHT 2.0 /* relative height below which to smoothen */ #define DEM_MAXHEIGHT 9000.0 /* max height of DEM (estimated from Everest/Olympus Mons) */ #define DEM_MAXDEPTH -10000 /* max depth of DEM */ -#define PLANET_SEALEVEL 0.0 /* sea level for flooded planet */ +#define PLANET_SEALEVEL 2500.0 /* sea level for flooded planet */ #define VENUS_NODATA_FACTOR 0.5 /* altitude to assign to DEM points without data (fraction of mean altitude) */ #define Z_SCALING_FACTOR 0.8 /* overall scaling factor of z axis for REP_PROJ_3D representation */ @@ -458,27 +449,27 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO #define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */ #define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */ #define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */ -#define DRAW_ARROW 0 /* set to 1 to draw arrow above sphere */ - -#define RSCALE 0.01 /* scaling factor of radial component */ -#define RSHIFT -0.01 /* shift in radial component */ -#define RMAX 2.0 /* max value of radial component */ -#define RMIN 0.5 /* min value of radial component */ -// #define COS_VISIBLE -1.1 /* limit on cosine of normal to shown facets */ -#define COS_VISIBLE -0.3 /* limit on cosine of normal to shown facets */ +#define DRAW_ARROW 1 /* set to 1 to draw arrow above sphere */ +#define RSCALE 0.2 /* scaling factor of radial component */ +#define RSHIFT -0.01 /* shift in radial component */ +#define RMAX 2.0 /* max value of radial component */ +#define RMIN 0.5 /* min value of radial component */ +#define COS_VISIBLE -0.75 /* limit on cosine of normal to shown facets */ /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ #define VMAX 100.0 /* max value of wave amplitude */ #define TEST_GRADIENT 0 /* print norm squared of gradient */ + #define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT) /* to save computing time, to be improved */ #define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED))) #define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)||(PRINT_AVERAGE_SPEED)) #define COMPUTE_PRESSURE ((ZPLOT == Z_EULER_PRESSURE)||(CPLOT == Z_EULER_PRESSURE)||(ZPLOT_B == Z_EULER_PRESSURE)||(CPLOT_B == Z_EULER_PRESSURE)) #define ASYM_SPEED_COLOR (VMEAN_SPEED == 0.0) +#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE) int block_sizes[NY]; /* table of block sizes for blocking around poles */ int block_numbers[NY]; /* table of block numbers for blocking around poles */ @@ -710,84 +701,166 @@ void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc { int i, j; double dx, dy; - - if (FORCE_FIELD == GF_COMPUTE_FROM_BC) + + switch (FORCE_FIELD) { - dx = (XMAX - XMIN)/(double)NX; - dy = (YMAX - YMIN)/(double)NY; + case (GF_COMPUTE_FROM_BC): + { + dx = (XMAX - XMIN)/(double)NX; + dy = (YMAX - YMIN)/(double)NY; - #pragma omp parallel for private(i,j) - for (i=1; idepth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; break; } + case (SH_SPHERE_CUBE): + { + rmax = 0.0; + /* compute distance from cube to origin, given by 0.5/rmax */ + if (vabs(wsphere[i*NY+j].x) > rmax) rmax = vabs(wsphere[i*NY+j].x); + if (vabs(wsphere[i*NY+j].y) > rmax) rmax = vabs(wsphere[i*NY+j].y); + if (vabs(wsphere[i*NY+j].z) > rmax) rmax = vabs(wsphere[i*NY+j].z); + h = 1.0 - 0.5/rmax; +// printf("h = %.3lg\n", h); + if (h < 0.0) h = 0.0; + rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + break; + } + case (SH_SPHERE_OCTAHEDRON): + { + rmax = 0.0; + /* compute distance from octahedron to origin, given by 0.5/rmax */ + rmax = 1.0/(vabs(wsphere[i*NY+j].x) + vabs(wsphere[i*NY+j].y) + vabs(wsphere[i*NY+j].z)); + h = 1.0 - 0.5/rmax; + printf("h = %.3lg\n", h); + if (h < 0.0) h = 0.0; + rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + break; + } + case (SH_SPHERE_DODECAHEDRON): + { + rmax = 0.0; + x1 = wsphere[i*NY+j].x; + y1 = wsphere[i*NY+j].y; + z1 = wsphere[i*NY+j].z; + /* compute distance from dodecahedron */ + d = vabs(phi1*y1 - phi2*z1) - sq3; + if (d > rmax) rmax = d; + d = vabs(-phi2*x1 + phi1*z1) - sq3; + if (d > rmax) rmax = d; + d = vabs(phi1*x1 - phi2*y1) - sq3; + if (d > rmax) rmax = d; + d = vabs(phi1*y1 + phi2*z1) - sq3; + if (d > rmax) rmax = d; + d = vabs(phi2*x1 + phi1*z1) - sq3; + if (d > rmax) rmax = d; + d = vabs(phi1*x1 + phi2*y1) - sq3; + if (d > rmax) rmax = d; + h = rmax; + printf("h = %.3lg\n", h); + if (h < 0.0) h = 0.0; + rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + break; + } + case (SH_SPHERE_ICOSAHEDRON): + { + rmax = 0.0; + x1 = wsphere[i*NY+j].x; + y1 = wsphere[i*NY+j].y; + z1 = wsphere[i*NY+j].z; + /* compute distance from dodecahedron */ + d = vabs(phi2*(x1+y1+z1)) - rico; + if (d > rmax) rmax = d; + d = vabs(phi2*(-x1+y1+z1)) - rico; + if (d > rmax) rmax = d; + d = vabs(phi2*(x1-y1+z1)) - rico; + if (d > rmax) rmax = d; + d = vabs(phi2*(x1+y1-z1)) - rico; + if (d > rmax) rmax = d; + d = vabs(phi3*x1 + phi1*z1) - rico; + if (d > rmax) rmax = d; + d = vabs(phi3*z1 + phi1*y1) - rico; + if (d > rmax) rmax = d; + d = vabs(phi3*y1 + phi1*x1) - rico; + if (d > rmax) rmax = d; + d = vabs(phi3*x1 - phi1*z1) - rico; + if (d > rmax) rmax = d; + d = vabs(phi3*z1 - phi1*y1) - rico; + if (d > rmax) rmax = d; + d = vabs(phi3*y1 - phi1*x1) - rico; + if (d > rmax) rmax = d; + h = rmax; + printf("h = %.3lg\n", h); + if (h < 0.0) h = 0.0; + rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + break; + } } } -double initialize_water_depth(t_rde rde[NX*NY]) +double initialize_water_depth(t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) { int i, j, ncircles; double dx, dy, min, max, pscal, norm, vz = 0.01; @@ -856,7 +1011,7 @@ double initialize_water_depth(t_rde rde[NX*NY]) #pragma omp parallel for private(i,j) for (i=0; i DSMOOTH)&&(j < NY-DSMOOTH)) + { + phi_out[0][i*NY+j] = eta - intstep*(u*etax + v*etay + eta*(ux + vy)); + phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax); + phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etax); + } + else + { + phi_out[0][i*NY+j] = SWATER_MIN_HEIGHT; + phi_out[1][i*NY+j] = 0.0; + phi_out[2][i*NY+j] = 0.0; + } + } + else + { + phi_out[0][i*NY+j] = eta - intstep*(u*etax + v*etay + eta*(ux + vy)); + phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax); + phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etay); + } if (VISCOSITY > 0.0) { @@ -1513,6 +1689,12 @@ void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], t_ vy = phi[2][i*NY+j]; break; } + case (E_SHALLOW_WATER): + { + vx = phi[1][i*NY+j]; + vy = phi[2][i*NY+j]; + break; + } } // v = module2(vx, vy); @@ -1842,7 +2024,7 @@ void animation() if (VARIABLE_DEPTH) { // water_depth = (t_swater_depth *)malloc(NX*NY*sizeof(t_swater_depth)); - max_depth = initialize_water_depth(rde); + max_depth = initialize_water_depth(rde, wsphere); printf("Max depth = %.3lg\n", max_depth); } @@ -1883,7 +2065,8 @@ void animation() /* initialize field */ // init_random(0.5, 0.4, phi, xy_in); // init_random(0.5, 0.25, phi, xy_in, wsphere); -// init_random_smoothed(0.5, 0.25, phi, xy_in, wsphere); +// init_random_smoothed(1.0, 0.1, phi, xy_in, wsphere); + // init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in) // init_coherent_state(0.0, 0.0, 10.0, 0.0, 0.1, phi, xy_in); @@ -1914,29 +2097,12 @@ void animation() // init_laminar_flow(double amp, double xmodulation, double ymodulation, double xperiod, double yperiod, double yshift, double density_mod, double *phi[NFIELDS], short int xy_in[NX*NY]) - init_laminar_flow_earth(0.04, phi, xy_in); - - /* high pressure systems, northern hemisphere */ - init_vortex_state_sphere_mod(1, 0.5, 2.5*PID, PID + 0.6, 0.03, 0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, 0.6, 3.7*PID, PID + 0.7, 0.035, 0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, 0.2, 2.5*PID, PID + 1.3, 0.01, 0.01, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, 0.5, 1.0*PID, PID + 1.0, 0.03, 0.02, phi, xy_in, wsphere); -// init_vortex_state_sphere_mod(1, 0.1, 3.0*PID, 0.9*PI, 0.01, 0.01, phi, xy_in, wsphere); - /* low pressure systems, northern hemisphere */ - init_vortex_state_sphere_mod(1, -0.5, 3.0*PID, PID + 0.6, 0.03, -0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, -0.2, 2.3*PID, PID + 1.1, 0.01, -0.02, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, -0.5, 1.2*PID, PID + 0.4, 0.01, -0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, -0.8, 0.3*PID, PID + 0.6, 0.01, -0.04, phi, xy_in, wsphere); - /* high pressure systems, southern hemisphere */ - init_vortex_state_sphere_mod(1, -0.6, 2.4*PID, PID - 0.7, 0.03, 0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, -0.6, 1.6*PID, PID - 0.7, 0.02, 0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, -0.4, 0.5*PID, PID - 0.7, 0.02, 0.04, phi, xy_in, wsphere); - /* low pressure systems, southern hemisphere */ - init_vortex_state_sphere_mod(1, 0.5, 3.4*PID, PID - 0.7, 0.03, -0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, 0.2, 2.0*PID, PID - 0.1, 0.04, -0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, 0.3, 1.0*PID, PID - 0.6, 0.03, -0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, 0.4, 0.1*PID, PID - 0.5, 0.03, -0.04, phi, xy_in, wsphere); - init_vortex_state_sphere_mod(1, 0.2, 0.1*PID, 0.1, 0.1, -0.01, phi, xy_in, wsphere); +// init_laminar_flow_earth(0.01, phi, xy_in); +// +// add_random_onefield_smoothed(0, 0.0, 0.03, phi, xy_in, wsphere); +// add_random_onefield_smoothed(1, 0.0, 0.07, phi, xy_in, wsphere); +// add_random_onefield_smoothed(2, 0.0, 0.07, phi, xy_in, wsphere); + // init_pressure_gradient_flow(flow_speed_schedule(0), 1.0 + PRESSURE_GRADIENT, 1.0 - PRESSURE_GRADIENT, phi, xy_in, bc_field); @@ -1948,14 +2114,13 @@ void animation() // init_vortex_state_sphere(0, amp, phishift, PID + thetashift, 0.15, 0.3, phi, xy_in, wsphere); -// init_gaussian_wave(-1.0, 0.0, 0.005, 0.25, SWATER_MIN_HEIGHT, phi, xy_in); +// init_gaussian_wave(0.7, 0.0, 0.0225, 0.1, SWATER_MIN_HEIGHT, phi, xy_in); +// init_gaussian_wave_sphere(4.0, PID, 0.05, 0.05, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); + + // init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in); -// init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in); - -// init_linear_blob(0.0, -0.75, 0.025, 0.0, 0.001, 0.25, 0.5, SWATER_MIN_HEIGHT, phi, xy_in); + init_swater_laminar_flow_sphere(0, -0.75, 0.0075, 6, 1, 0.0025, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); -// init_elliptical_vibration(0.0, 0.0, 0.0075, LAMBDA, 1.0, 0.0003, 0.1, 1.0, SWATER_MIN_HEIGHT, phi, xy_in); - // add_gaussian_wave(-1.6, -0.5, 0.015, 0.25, SWATER_MIN_HEIGHT, phi, xy_in); if (SMOOTHBLOCKS) for (k=0; k NMAXSEGMENTS) + { + printf("NMAXSEGMENTS too small"); + return(0); + } + if (nbelts + 1 > NMAXBELTS) + { + printf("NMAXBELTS too small"); + return(0); + } -void init_segment_config(t_segment segment[NMAXSEGMENTS]) + segment[n].x1 = x1 - ty*width; + segment[n].y1 = y1 + tx*width; + + for (i=0; i<7; i++) + { + segment[n+1+i].x1 = x2 + width*sin(-angle + (double)i*PI/6.0); + segment[n+1+i].y1 = y2 + width*cos(-angle + (double)i*PI/6.0); + } + for (i=7; i<14; i++) + { + segment[n+1+i].x1 = x1 + width*sin(-angle + (double)(i-1)*PI/6.0); + segment[n+1+i].y1 = y1 + width*cos(-angle + (double)(i-1)*PI/6.0); + } + + for (i=0; i<13; i++) + { + segment[n+i].x2 = segment[n+i+1].x1; + segment[n+i].y2 = segment[n+i+1].y1; + } + segment[n+13].x2 = segment[n].x1; + segment[n+13].y2 = segment[n].y1; + + for (i=n; i n + 13) iplus = n; + angle1 = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID; + angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID; + if (angle2 < angle1) angle2 += DPI; + segment[i].angle1 = angle1; + segment[i].angle2 = angle2; + +// printf("i = %i, iplus = %i, iminus = %i, angle1 = %.0f, angle2 = %.0f\n", i, iplus, iminus, angle*360.0/DPI, angle2*360.0/DPI); + } + + nsegments += 14; + + conveyor_belt[nbelts].x1 = x1; + conveyor_belt[nbelts].x2 = x2; + conveyor_belt[nbelts].y1 = y1; + conveyor_belt[nbelts].y2 = y2; + conveyor_belt[nbelts].speed = speed; + conveyor_belt[nbelts].width = width; + conveyor_belt[nbelts].position = 0.0; + conveyor_belt[nbelts].length = length; + conveyor_belt[nbelts].angle = angle; + conveyor_belt[nbelts].tx = tx; + conveyor_belt[nbelts].ty = ty; + nbelts++; + + return(1); +} + +void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[NMAXBELTS]) /* initialise linear obstacle configuration */ { int i, j, cycle = 0, iminus, iplus, nsides, n, concave = 1; double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, yy, dy, ca, sa, padding; double lfoot[NTREES][2], rfoot[NTREES][2]; + /* set default to no conveyor */ + for (i=0; i 0.0)) fy += y - width; + if (y > 0.0) fy += y - width; + else fy += y + width; + } + else if ((x >= mu_i)&&(x < mu_i+width)) + { + r = module2(x-mu_i, y) + 1.0e-10; + if (r < width) + { + fx -= (width - r)*(x-mu_i)/r; + fy -= (width - r)*y/r; + } + } + else if ((x <= -mu_i)&&(x > -mu_i-width)) + { + r = module2(x+mu_i, y) + 1.0e-10; + if (r < width) + { + fx -= (width - r)*(x+mu_i)/r; + fy -= (width - r)*y/r; + } + } + + torque += x*fy - y*fx; + } + + if (charge) /* add Coulomb force between ends */ + { + cfact = CHARGE*KREPEL/KSPRING_OBSTACLE; + + r = module2(x-mu_i,y) + 1.0e-2; + r3 = r*r*r; + fx += cfact*(double)s*(x+mu_i)/r3; + fy += cfact*(double)s*y/r3; + + torque += mu_i*fx; + + r = module2(x+mu_i,y) + 1.0e-2; + r3 = r*r*r; + fx -= cfact*(double)s*(x+mu_i)/r3; + fy -= cfact*(double)s*y/r3; + + torque -= mu_i*fy; + } + } + + ffm[0] = cb*fx - sb*fy; + ffm[1] = sb*fx + cb*fy; + ffm[2] = torque; +} + +void polygon_interaction(t_particle part_i, t_particle part_k, double ca, double sa, double distance, double ffm[3], int charge) +/* compute interaction for I_POLYGON interaction */ +/* computes force and torque of particle k on particle i */ +{ + int s, k, s0, s1; + double x0, y0, x, y, f, fx, fy, cb, sb, cab, sab, r, r0, z2, phi, ckt, skt, d; + double torque, mu_i, mu_k, r3, width, gamma_beta, rot, cr, sr, xx[NPOLY], yy[NPOLY], dd[NPOLY], z, z0, z02, d1, dmin; + static double theta, twotheta, ctheta, stheta, cornerx[NPOLY+1], cornery[NPOLY+1], nx[NPOLY], ny[NPOLY]; + static int first = 1; + + if (first) + { + theta = PI/(double)NPOLY; + twotheta = 2.0*theta; + ctheta = cos(theta); + stheta = sin(theta); + + for (s=0; s width) f = width; + fx -= f*nx[k]; + fy -= f*ny[k]; + torque += x*fy - y*fx; + } + else for (s1=0; s1 width) f = width; + fx -= f*(x - mu_i*cornerx[s1])/d1; + fy -= f*(y - mu_i*cornery[s1])/d1; + } + torque += x*fy - y*fx; + } + } + } + + ffm[0] = cb*fx - sb*fy; + ffm[1] = sb*fx + cb*fy; + ffm[2] = torque; +} int dna_charge(int type1, int type2, int coulomb) /* returns -1 for matching nucleotide ends, 1 for non-matching, and 0 else */ @@ -3887,7 +4331,7 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque, /* compute repelling force and torque of particle #k on particle #i */ /* returns 1 if distance between particles is smaller than NBH_DIST_FACTOR*MU */ { - double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp, charge, f1; + double x1, y1, x2, y2, r, f, angle, aniso, fx, fy, ff[2], dist_scaled, spin_f, ck, sk, ck_rel, sk_rel, alpha, amp, charge, f1, ffm[3]; static double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN); int wwrapx, wwrapy, twrapx, twrapy; @@ -4022,7 +4466,7 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque, charge = particle[i].charge*particle[k].charge; f = -100.0*krepel*charge/(1.0e-12 + distance*distance); if (charge <= 0.0) - f += 0.01*krepel*lennard_jones_force(distance, particle[k]); + f += COULOMB_LJ_FACTOR*krepel*lennard_jones_force(distance, particle[k]); force[0] = f*ca; force[1] = f*sa; break; @@ -4096,6 +4540,60 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque, force[1] = f*sa; break; } + case (I_SEGMENT): + { + segment_interaction(particle[i], particle[k], ca, sa, distance, ffm, 0); + force[0] = KSPRING_OBSTACLE*ffm[0]; + force[1] = KSPRING_OBSTACLE*ffm[1]; + *torque = KTORQUE*ffm[2]; + + f = 0.1*krepel*lennard_jones_force(distance, particle[k]); + force[0] += f*ca; + force[1] += f*sa; + + break; + } + case (I_SEGMENT_CHARGED): + { + segment_interaction(particle[i], particle[k], ca, sa, distance, ffm, 1); + force[0] = KSPRING_OBSTACLE*ffm[0]; + force[1] = KSPRING_OBSTACLE*ffm[1]; + *torque = KTORQUE*ffm[2]; + +// printf("Force between particles %i and %i: (%.3lg, %.3lg)\n", i, k, force[0], force[1]); + +// f = 0.1*krepel*lennard_jones_force(distance, particle[k]); +// force[0] += f*ca; +// force[1] += f*sa; + + break; + } + case (I_POLYGON): + { + polygon_interaction(particle[i], particle[k], ca, sa, distance, ffm, 0); + force[0] = KSPRING_OBSTACLE*ffm[0]; + force[1] = KSPRING_OBSTACLE*ffm[1]; + *torque = KTORQUE*ffm[2]; + + f = 0.1*krepel*lennard_jones_force(distance, particle[k]); + force[0] += f*ca; + force[1] += f*sa; + + break; + } + case (I_POLYGON_ALIGN): + { + polygon_interaction(particle[i], particle[k], ca, sa, distance, ffm, 0); + force[0] = KSPRING_OBSTACLE*ffm[0]; + force[1] = KSPRING_OBSTACLE*ffm[1]; + *torque = KTORQUE*ffm[2]; + + f = 0.1*krepel*lennard_jones_force(distance, particle[k]); + force[0] += f*ca; + force[1] += f*sa; + + break; + } } if (ROTATION) @@ -4160,6 +4658,36 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque, break; } + case (I_SEGMENT): + { + /* Do nothing, torque already computed */ + break; + } + case (I_SEGMENT_CHARGED): + { + /* Do nothing, torque already computed */ + break; + } + case (I_POLYGON): + { + /* Do nothing, torque already computed */ + break; + } + case (I_POLYGON_ALIGN): + { + spin_f = particle[i].spin_freq; + if (twrapx||twrapy) + *torque += sin(spin_f*(-particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); + else + { + if (NPOLY%2 == 0) + *torque += sin(spin_f*(particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); + else + *torque -= sin(spin_f*(particle[k].angle - particle[i].angle))/(1.0e-8 + dist_scaled*dist_scaled); + + } + break; + } default: { spin_f = particle[i].spin_freq; @@ -4333,11 +4861,16 @@ void print_partners(int i, t_particle particle[NMAXCIRCLES]) // printf("(x, y) = (%.3lg, %.3lg), radius = %.3lg\n", particle[i].xc, particle[i].yc, particle[i].radius); printf("%i partners: ", particle[i].npartners); + fprintf(lj_log, "%i partners: ", particle[i].npartners); for (p=0; p nmolecules) nmolecules = m + 1; np = molecule[m].nparticles; +// printf("np = %i\n", np); if (np < NPARTNERS+1) { molecule[m].particle[np] = i; @@ -5138,14 +5885,20 @@ void init_molecule_data(t_particle particle[NMAXCIRCLES], t_molecule molecule[NM molecule[m].added = particle[i].added; } printf("Found %i molecules\n", nmolecules); + fprintf(lj_log, "Found %i molecules\n", nmolecules); /* for debugging */ for (m=0; m 0.0) + { + hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*ej/PARTICLE_EMAX; + if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN; + if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX; + } + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } case (P_LOG_EMEAN): { ej = particle.emean; @@ -5542,7 +6311,9 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl case (P_DIRECT_EMEAN): { hue = particle.dirmean + COLOR_HUESHIFT*PI; +// printf("dirmean = %.3lg\n", particle.dirmean); if (hue > DPI) hue -= DPI; + if (hue < 0.0) hue += DPI; hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*(hue)/DPI; ej = particle.emean; if (ej < 0.5*PARTICLE_EMAX) lum = 2.0*ej/PARTICLE_EMAX; @@ -5597,8 +6368,41 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl } case (P_CLUSTER): { - cluster = (double)particle.cluster/(double)ncircles; - hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*cluster; +// cluster = (double)(particle.cluster)/(double)(ncircles); + ccluster = (double)(particle.cluster_color)/(double)(ncircles); + ccluster -= (double)((int)ccluster); + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*ccluster; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_CLUSTER_SIZE): + { +// cluster = (double)(particle.cluster)/(double)(ncircles); + ccluster = 1.0 - 5.0/((double)particle.cluster_size + 4.0); + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*ccluster; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_CLUSTER_SELECTED): + { + cl = particle.cluster; + if (cluster[cl].selected) hue = COLOR_HUE_CLUSTER_SELECTED; + else hue = COLOR_HUE_CLUSTER_NOT_SELECTED; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } + case (P_COLLISION): + { + hue = (double)particle.collision; + if (hue > 0.0) hue = atan(0.25*(0.03*hue + 1.0))/PID; +// { +// hue += 10.0; +// hue *= 1.0/(40.0 + hue); +// } + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; *radius = particle.radius; *width = BOUNDARY_WIDTH; break; @@ -5715,6 +6519,20 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_CLUSTER); break; } + case (P_CLUSTER_SIZE): + { + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_CLUSTER_SIZE); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_CLUSTER_SIZE); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_CLUSTER_SIZE); + break; + } + case (P_CLUSTER_SELECTED): + { + hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_CLUSTER_SELECTED); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_CLUSTER_SELECTED); + hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_CLUSTER_SELECTED); + break; + } default: { hsl_to_rgb(hue, 0.9, 0.5, rgb); @@ -5821,7 +6639,7 @@ void draw_one_triangle(t_particle particle, int same_table[9*HASHMAX], int p, in } -void draw_triangles(t_particle particle[NMAXCIRCLES], int plot) +void draw_triangles(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], int plot) /* fill triangles between neighboring particles */ { int i, j, k, p, q, t0, tj, tmax, nsame = 0, same_table[9*HASHMAX], width; @@ -5854,7 +6672,7 @@ void draw_triangles(t_particle particle[NMAXCIRCLES], int plot) } else { - compute_particle_colors(particle[i], plot, rgb, rgbx, rgby, &radius, &width, particle); + compute_particle_colors(particle[i], cluster, plot, rgb, rgbx, rgby, &radius, &width, particle); } glColor3f(rgb[0], rgb[1], rgb[2]); @@ -6301,6 +7119,13 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK)) draw_colored_rhombus(xc1, yc1, radius, angle + APOLY*PID, rgb); + else if (particle.interaction == I_SEGMENT) + draw_colored_rotated_rect(xc1, yc1, particle.radius, 0.1*particle.radius, angle + APOLY*PID, rgb); + else if (particle.interaction == I_SEGMENT_CHARGED) + { + draw_colored_rotated_rect(xc1, yc1, particle.radius, 0.1*particle.radius, angle + APOLY*PID, rgb); + /* TODO: draw ends */ + } else if (cont) { if (nsides == NSEG) draw_colored_circle_precomp(xc1, yc1, radius, rgb); @@ -6357,6 +7182,13 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius, if ((particle.interaction == I_LJ_QUADRUPOLE)||(particle.interaction == I_LJ_DIPOLE)||(particle.interaction == I_VICSEK)||(particle.interaction == I_VICSEK_REPULSIVE)||(particle.interaction == I_VICSEK_SPEED)||(particle.interaction == I_VICSEK_SHARK)) draw_rhombus(xc1, yc1, radius, angle + APOLY*PID); + else if (particle.interaction == I_SEGMENT) + draw_rotated_rect(xc1, yc1, particle.radius, 0.1*particle.radius, angle + APOLY*PID); + else if (particle.interaction == I_SEGMENT_CHARGED) + { + draw_rotated_rect(xc1, yc1, particle.radius, 0.1*particle.radius, angle + APOLY*PID); + /* TODO: draw ends */ + } else if (cont) { if (nsides == NSEG) draw_circle_precomp(xc1, yc1, radius); @@ -6402,7 +7234,7 @@ void draw_collisions(t_collision *collisions, int ncollisions) } else y1 = y; - draw_colored_polygon(x1, y1, 5.0*MU, NSEG, 0.0, rgb); + draw_colored_polygon(x1, y1, COLLISION_RADIUS*MU, NSEG, 0.0, rgb); collisions[i].time--; } @@ -6624,9 +7456,9 @@ void color_background(t_particle particle[NMAXCIRCLES], int bg_color, t_hashgrid } -void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_collision *collisions, int ncollisions, int bg_color, t_hashgrid hashgrid[HASHX*HASHY], t_lj_parameters params) +void draw_particles(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], int plot, double beta, t_collision *collisions, int ncollisions, int bg_color, t_hashgrid hashgrid[HASHX*HASHY], t_lj_parameters params) { - int i, j, k, m, width, nnbg, nsides; + int i, j, k, m, width, nnbg, nsides, cl, p, q; 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]; @@ -6726,7 +7558,7 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c } /* fill triangles between particles */ - if (FILL_TRIANGLES) draw_triangles(particle, plot); + if (FILL_TRIANGLES) draw_triangles(particle, cluster, plot); /* draw collision discs */ if ((REACTION_DIFFUSION)&&(ncollisions > 0)) draw_collisions(collisions, ncollisions); @@ -6734,7 +7566,7 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c /* determine particle color and size */ for (j=0; j 0.0) + pos0 = conveyor_belt[i].position - bsegment*(double)((int)(conveyor_belt[i].position/bsegment)); + else + { + pos = -conveyor_belt[i].position; + pos0 = pos - bsegment*(double)((int)(pos/bsegment)); + pos0 = bsegment - pos0; + } + + while (pos0 < conveyor_belt[i].length) + { + x = x1 + tx*pos0; + y = y1 + ty*pos0; + draw_line(x - r*ty, y + r*tx, x - bwidth*ty, y + bwidth*tx); + pos0 += bsegment; + } + while (pos0 < conveyor_belt[i].length + PI*bwidth) + { + pos = (pos0 - conveyor_belt[i].length)/bwidth; + angle = PID + conveyor_belt[i].angle - pos; + ca = cos(angle); + sa = sin(angle); + draw_line(x2 + r*ca, y2 + r*sa, x2 + bwidth*ca, y2 + bwidth*sa); + pos0 += bsegment; + } + while (pos0 < 2.0*conveyor_belt[i].length + PI*bwidth) + { + pos = pos0 - conveyor_belt[i].length - PI*bwidth; + x = x2 - tx*pos; + y = y2 - ty*pos; + draw_line(x + r*ty, y - r*tx, x + bwidth*ty, y - bwidth*tx); + pos0 += bsegment; + } + while (pos0 < 2.0*conveyor_belt[i].length + DPI*bwidth) + { + pos = (pos0 - 2.0*conveyor_belt[i].length - PI*bwidth)/bwidth; + angle = -PID + conveyor_belt[i].angle - pos; + ca = cos(angle); + sa = sin(angle); + draw_line(x1 + r*ca, y1 + r*sa, x1 + bwidth*ca, y1 + bwidth*sa); + pos0 += bsegment; + } + + } } switch (BOUNDARY_COND) { @@ -7840,7 +8786,7 @@ void print_particles_speeds(t_particle particle[NMAXCIRCLES]) 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 i, k; - 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, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin, charge, d2; + 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, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin, charge, d2, speed; /* compute force from fixed circular obstacles */ if (ADD_FIXED_OBSTACLES) for (i=0; i 0.5*(BCYMAX - BCYMIN)) dy -= (BCYMAX-BCYMIN); else if (dy < -0.5*(BCYMAX - BCYMIN)) dy += (BCYMAX-BCYMIN); - r = module2(dx, dy); - if (r < 1.0e-10) r = 1.0e-10; - if (r > 2.0*eq_distance) r = 2.0*eq_distance; -// if (r > 1.5*eq_distance) r = 1.5*eq_distance; - ca = dx/r; - sa = dy/r; - /* TODO: adjust max distance */ - if (r < 1.5*eq_distance) force = KSPRING_PAIRS*(r - eq_distance); - else + + if ((PAIR_PARTICLES)&&(PAIRING_TYPE == POLY_SEG_POLYGON)) { + rmax = 1.0*particle[j].radius; + for (p=-1; p<=1; p+=2) + for (q=-1; q<=1; q+=2) + { + x1 = particle[n].xc + (double)p*particle[n].radius*cos(particle[n].angle); + y1 = particle[n].yc + (double)p*particle[n].radius*sin(particle[n].angle); + + x2 = particle[j].xc + (double)q*particle[n].radius*cos(particle[j].angle); + y2 = particle[j].yc + (double)q*particle[n].radius*sin(particle[j].angle); + + r = module2(x2-x1, y2-y1); + if ((r>0.0)&&(r < rmax)) + { +// f[0] += KSPRING_PAIRS*(x1-x2)/r; +// f[1] += KSPRING_PAIRS*(y1-y2)/r; + f[0] = 1.0e8*(x1-x2)/r; + f[1] = 1.0e8*(y1-y2)/r; + } + } + } +// else + { + r = module2(dx, dy); + if (r < 1.0e-10) r = 1.0e-10; + if (r > 2.0*eq_distance) r = 2.0*eq_distance; +// if (r > 1.5*eq_distance) r = 1.5*eq_distance; + ca = dx/r; + sa = dy/r; + + /* TODO: adjust max distance */ + if (r < 1.5*eq_distance) force = KSPRING_PAIRS*(r - eq_distance); + else + { // printf("Dissociating partners %i and %i because max distance exceeded\n", j, n); - force = 0.0; + force = 0.0; // dissociate_particles_findp(j, n, particle); // dissociate_molecule(j, n, particle); - } + } - f[0] = force*ca; - f[1] = force*sa; + f[0] = force*ca; + f[1] = force*sa; + + /* TEST */ + f[0] -= particle[j].vx*DAMPING; + f[1] -= particle[j].vy*DAMPING; + } if (ROTATION) { - sangle = sin(SPIN_INTER_FREQUENCY*(particle[n].angle - particle[j].angle)); - if (sangle > 0.0) *torque = KTORQUE_PAIRS*(1.0 + sangle); - else *torque = KTORQUE_PAIRS*(-1.0 + sangle); - - if (COMPUTE_PAIR_TORQUE) + *torque = 0.0; + if ((REACTION_DIFFUSION)&&(RD_REACTION == CHEM_POLYGON_AGGREGATION)) { - torque2 = KTORQUE_PAIR_ANGLE*sangle; - f[0] -= torque2*sa; - f[1] += torque2*ca; +// r = module2(dx, dy); +// if (r < 2.0*MU) +// { +// if (NPOLY%2 == 0) +// sangle = sin((double)NPOLY*(particle[n].angle - particle[j].angle)); +// else sangle = sin((double)NPOLY*(particle[n].angle - particle[j].angle) - PI); +// sangle = sin((double)NPOLY*(particle[n].angle - particle[j].angle - eq_angle)); +// *torque += KTORQUE_PAIRS*sangle; + + /* stabilize relative angle */ + q = particle[j].partner[0]; + if ((particle[j].npartners == 1)&&(particle[q].npartners == 1)) + { + sangle = sin((double)NPOLY*(particle[n].angle - particle[j].angle - eq_angle)); + *torque += KTORQUE_PAIRS*sangle; + phi = argument(dx, dy); + sangle = sin((double)NPOLY*(phi - particle[j].angle) - PI); + *torque += KTORQUE_PAIRS*sangle; + } + + /* TEST: add some damping to rotation */ + *torque -= particle[j].omega*DAMPING_ROT; + +// sangle = sin((double)NPOLY*(phi - particle[j].angle) - PI); +// *torque += KTORQUE_PAIRS*sangle; +// } + } + else if ((PAIR_PARTICLES)&&(PAIRING_TYPE == POLY_SEG_POLYGON)) + { + /* TODO */ + sangle = sin(particle[n].angle - particle[j].angle - eq_angle); +// if (sangle > 0.0) *torque = KTORQUE_PAIRS*(1.0 + sangle); +// else *torque = KTORQUE_PAIRS*(-1.0 + sangle); + *torque = KTORQUE_PAIRS*sangle; + + +// *torque = 0.0; + +// if (COMPUTE_PAIR_TORQUE) +// { +// torque2 = KTORQUE_PAIR_ANGLE*sangle; +// f[0] -= torque2*sa; +// f[1] += torque2*ca; +// } + } + else + { + sangle = sin(SPIN_INTER_FREQUENCY*(particle[n].angle - particle[j].angle)); + if (sangle > 0.0) *torque = KTORQUE_PAIRS*(1.0 + sangle); + else *torque = KTORQUE_PAIRS*(-1.0 + sangle); + + if (COMPUTE_PAIR_TORQUE) + { + torque2 = KTORQUE_PAIR_ANGLE*sangle; + f[0] -= torque2*sa; + f[1] += torque2*ca; + } } } } @@ -8521,7 +9562,7 @@ void compute_partner_force(int j, int n, double eq_distance, double f[2], double void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY]) /* compute force from other particles on particle j */ { - int i0, j0, m0, k, l, m, q, close, n, test; + int i0, j0, m0, k, l, m, q, close, n, test, different_cluster = 1; double fx = 0.0, fy = 0.0, force[2], torque = 0.0, torque_ij, x, y; particle[j].neighb = 0; @@ -8530,11 +9571,15 @@ void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLE for (k=0; k= cluster[j].nparticles) + { + il = i; + is = j; + pl = pi; + ps = pj; + } + else + { + il = j; + is = i; + pl = pj; + ps = pi; + } + ns = cluster[is].nparticles; + nl = cluster[il].nparticles; + + if (ns + nl >= NMAXPARTINCLUSTER) + { + printf("Warning: NMAXPARTINCLUSTER is too small\n"); + printf("Try increasing to %i\n", cluster[i].nparticles + cluster[j].nparticles); + return(0); + } + + printf("Merging clusters %i and %i, having %i and %i particles\n", il, is, nl, ns); + fprintf(lj_log, "Merging clusters %i and %i, having %i and %i particles\n", il, is, nl, ns); + + fprintf(lj_log, "Large cluster %i:\n", il); + for (p=0; p= alpha2) particle[q].angle -= alpha2; + while (particle[q].angle < 0.0) particle[q].angle += alpha2; + } + q = cluster[il].particle[0]; + newangle = particle[q].angle; + cluster[il].angle = newangle; +// if (NPOLY%2==1) newangle += alpha; + for (p=nl; p maxpartners) return(ncollisions); + + m1 = 1.0/particle[i].mass_inv; + for (p=0; p maxpartners) return(0); + + if ((particle[k].active)&&(particle[k].type == type2)&&(nq <= maxpartners)) + { + distance = module2(particle[i].deltax[p], particle[i].deltay[p]); + rel_angle = argument(particle[i].deltax[p], particle[i].deltay[p]); + poly_condition = (distance < REACTION_DIST*(ri+rk)); + + rorient = particle[k].angle - particle[i].angle; + if (NPOLY%2 == 1) rorient -= alpha; + norient = (int)(rorient/alpha2 + 0.5); + delta = rorient - norient*alpha2; + poly_condition *= (vabs(delta) < DELTAMAX); + + if ((poly_condition)&&((double)rand()/RAND_MAX < reaction_prob)) + { + /* find nearest orientation facing particle k */ + rorient_new = norient*alpha2; + nangle = (int)((particle[i].angle - rel_angle)/alpha2); + anglei_new = rel_angle + (double)(2*nangle+1)*alpha; + + m2 = 1.0/particle[k].mass_inv; + for (q=0; q (%.3lg, %.3lg)\n", particle[q].xc, particle[q].yc, newx, newy); + dx = newx - particle[q].xc; + dy = newy - particle[q].yc; + r = module2(dx, dy); + +// printf("Particle %i distance to target = %.3lg\n", k, r); +// if ((t == nsteps-1)&&(q == 470)&&(k == 472)) fprintf(lj_log, "Particle %i relative distance to target %i = %.5lg\n", q, k, r/eqdist); + + fx += r*dx; + fy += r*dy; + } + else + { + eqdist2 = 1.0*(particle[k].radius + particle[q].radius)*sa; + if (dist < eqdist2) + { + fx -= 0.1*(eqdist2 - dist)*dx/dist; + fy -= 0.1*(eqdist2 - dist)*dy/dist; + } + } + } + + if ((kill_overlays)&&(dist < REPAIR_MIN_DIST*eqdist)) + { + nq = particle[q].npartners; + np = particle[k].npartners; + if (nq == 1) particle[q].active = 0; + if (np == 1) particle[k].active = 0; +// if (np > nq) particle[q].active = 0; +// else particle[k].active = 0; + + } + } + } + particle[q].xc += kspring*fx*DT_PARTICLE; + particle[q].yc += kspring*fy*DT_PARTICLE; + } + } +} + + +int chem_multi_glue_polygon2(int i, int maxpartners, t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], t_collision *collisions, int ncollisions, double reaction_prob, int mergeclusters, int singlecluster) +/* simplified version of chem_multi_glue_molecule without triangle and similar */ +/* conditions, but taking polygon geometry into account */ +/* version 2, for CHEM_POLYGON_CLUSTER interaction */ +{ + int k, p, q, p1, q1, p2, np, nq, n, m, closeby = 0, reaction = 0, jp[NMAXPARTNERS], jq[NMAXPARTNERS], jtot[NMAXPARTNERS], r, kk, different, np2, moli, molk, mp, newpartner, pp, type1, poly_condition, nangle, norient, cl1, cl2, cp, cq, newcluster, nmin, clmin, clmax, newreaction; + double distance, angle, move_factor, xnew, ynew, vxnew, vynew, m1, m2, mr1, mr2, deltav, x, y, ri, rk, alpha, alpha2, ca, rel_angle, rorient, delta, rorient_new, anglei_new, deltax, deltay, x1, y1, aratio, delta2, eqdistance, dx, dy, dalpha1, dalpha2, ssize1, ssize2; + static short int first = 1; + static int scluster; + + if ((singlecluster)&&(first)) + { + scluster = 0; + while (cluster[scluster].active == 0) scluster = rand()%ncircles; + cluster[scluster].selected = 1; + particle[scluster].collision = 1; + printf("Selected cluster %i as single cluster\n", scluster); + first = 0; + } + + type1 = particle[i].type; + np = particle[i].npartners; + if (np > maxpartners) return(ncollisions); + + cl1 = particle[i].cluster; + cp = cluster[cl1].nparticles; + m1 = cluster[cl1].mass; + + alpha = PI/(double)NPOLY; + alpha2 = 2.0*alpha; + ca = cos(alpha); + ri = particle[i].radius*ca; + rk = particle[k].radius*ca; + + p = 0; + while ((p maxpartners) return(ncollisions); + + if ((!reaction)&&(particle[k].active)&&(np+nq+1 < maxpartners)&&(cl2!=cl1)&&(nmin <= SMALL_NP_MAXSIZE)&&(clmin <= SMALL_CLUSTER_MAXSIZE)) + { + cq = cluster[cl2].nparticles; +// poly_condition = (cl1!=cl2); + + /* test distance */ + distance = module2(particle[i].deltax[p], particle[i].deltay[p]); + rel_angle = argument(particle[i].deltax[p], particle[i].deltay[p]); + poly_condition = (distance < REACTION_DIST*(ri+rk)); + + /* test relative angle */ + rorient = particle[k].angle - particle[i].angle; + if (NPOLY%2 == 1) rorient -= alpha; + + while (rorient > alpha) rorient -= alpha2; + while (rorient < -alpha) rorient += alpha2; + poly_condition *= (vabs(rorient) < DELTAMAX); + +// norient = (int)(rorient/alpha2 + 0.5 + 2.0*(double)NPOLY); +// delta = rorient - (double)(norient-2*NPOLY)*alpha2; +// poly_condition *= (vabs(delta) < DELTAMAX); + + /* test orientation */ + aratio = 0.5*vabs(rel_angle)/alpha; + delta2 = aratio - (double)((int)aratio); + poly_condition *= (vabs(delta2 - 0.5) < 0.5*DELTAMAX*alpha); + + /* option singlecluster: only allow merger if one of the clusters is selected cluster */ + if (singlecluster) poly_condition *= ((cluster[cl1].selected)||(cluster[cl2].selected)||(clmax <= NOTSELECTED_CLUSTER_MAXSIZE)); + + newreaction = 0; + + if ((poly_condition)&&((double)rand()/RAND_MAX < reaction_prob)) + { + /* find nearest orientation facing particle k */ + rorient_new = norient*alpha2; + nangle = (int)((particle[i].angle - rel_angle)/alpha2); + anglei_new = rel_angle + (double)(2*nangle+1)*alpha; + + m2 = cluster[cl2].mass; + + mr1 = m1/(m1 + m2); + mr2 = 1.0 - mr1; + + if ((np == 0)&&(nq == 0)) + { + reaction = 1; + newreaction = 1; + printf("Merging molecule %i with particle %i\n", i, k); + + eqdistance = ri + rk; + particle[i].angle = anglei_new; + particle[k].angle = particle[i].angle + rorient_new; + if (NPOLY%2 == 1) particle[k].angle += alpha; + + dx = 0.5*(distance - eqdistance)*cos(rel_angle); + dy = 0.5*(distance - eqdistance)*sin(rel_angle); + + translate_cluster(cl1, cluster, particle, dx, dy); + translate_cluster(cl2, cluster, particle, -dx, -dy); + + particle[i].npartners++; + particle[i].partner[0] = k; + + particle[k].npartners++; + particle[k].partner[0] = i; + + /* equalize speeds */ + vxnew = mr1*particle[i].vx + mr2*particle[k].vx; + vynew = mr1*particle[i].vy + mr2*particle[k].vy; + + particle[i].vx = vxnew; + particle[i].vy = vynew; + particle[k].vx = vxnew; + particle[k].vy = vynew; + + if (mergeclusters) + newcluster = merge_clusters(particle[i].cluster, particle[k].cluster, i, k, particle, cluster, 0, (particle[i].flip == particle[k].flip)); + } +// else if ((np+nq+2 < NMAXPARTNERS)&&(vabs(delta) < 0.5*DELTAMAX)) + else if ((np+nq+2 < NMAXPARTNERS)&&(cp+cq <= CLUSTER_MAXSIZE)) + { + deltav = module2(particle[i].vx - particle[k].vx, particle[i].vy - particle[k].vy); +// printf("Delta v is %.3lg\n", deltav); + + if (deltav < DELTAVMAX) + { + reaction = 1; + newreaction = 1; + printf("Pairing clusters containing particles %i and %i\n", i, k); + printf("np = %i, nq = %i\n", np, nq); + fprintf(lj_log, "Pairing clusters containing particles %i and %i\n", i, k); + fprintf(lj_log, "np = %i, nq = %i\n", np, nq); + + eqdistance = ri + rk; + + dx = 0.5*(distance - eqdistance)*cos(rel_angle); + dy = 0.5*(distance - eqdistance)*sin(rel_angle); + dalpha1 = anglei_new - particle[i].angle; + dalpha2 = anglei_new + rorient_new - particle[k].angle; + if (NPOLY%2 == 1) dalpha2 += alpha; + /* round to closest multiple of alpha2 */ + while (dalpha1 > 0.5*DELTAMAX) dalpha1 -= alpha2; + while (dalpha1 < -0.5*DELTAMAX) dalpha1 += alpha2; + while (dalpha2 > 0.5*DELTAMAX) dalpha2 -= alpha2; + while (dalpha2 < -0.5*DELTAMAX) dalpha2 += alpha2; + + translate_cluster(cl1, cluster, particle, dx, dy); + translate_cluster(cl2, cluster, particle, -dx, -dy); + + ssize1 = 0.02*sqrt((double)cluster[cl1].nparticles); + ssize2 = 0.02*sqrt((double)cluster[cl2].nparticles); + + if (vabs(dalpha1) + ssize1 < 0.5*DELTAMAX) + rotate_cluster_around_particle(cl1, i, cluster, particle, dalpha1); + else printf("Did not rotate cluster of size %.3lg\n", ssize1); + if (vabs(dalpha2) + ssize2 < 0.5*DELTAMAX) + rotate_cluster_around_particle(cl2, k, cluster, particle, dalpha2); + else printf("Did not rotate cluster of size %.3lg\n", ssize2); + + particle[i].npartners++; + particle[i].partner[np] = k; + + particle[k].npartners++; + particle[k].partner[nq] = i; + + /* equalize speeds */ + vxnew = mr1*particle[i].vx + mr2*particle[k].vx; + vynew = mr1*particle[i].vy + mr2*particle[k].vy; + + for (p1=0; p1 SMALL_NP_MAXSIZE)||(clmin > SMALL_CLUSTER_MAXSIZE)) +// { +// printf("No merger - Parameters: nmin = %i, clmin = %i\n", nmin, clmin); +// } + p++; + } + return(ncollisions); +} + + int chem_split_molecule(int i, t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_collision *collisions, int ncollisions) /* split molecule containing particle i from all other molecules */ { @@ -11466,10 +13358,11 @@ int check_dna_pairing(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAX } -int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_collision *collisions, int ncollisions, int *particle_numbers, int time, double *delta_e) +int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], t_collision *collisions, int ncollisions, int *particle_numbers, int time, double *delta_e) /* update the types in case of reaction-diffusion equation */ { int i, j, k, n, n3, n4, p, type, atype, btype, oldncollisions, delta_n; + short int reacted; double distance, rnd, p1, p2, reac_dist; static double inv_masses[RD_TYPES+1], radii[RD_TYPES+1]; static int first = 1, repair_counter = 0; @@ -12226,6 +14119,81 @@ int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCL // sleep(3); + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_POLYGON_AGGREGATION): + { + for (i=0; i= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + + /* update cluster color scheme */ + if ((PLOT == P_CLUSTER)||(PLOT_B == P_CLUSTER)) + update_cluster_color(particle); + + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + } + case (CHEM_POLYGON_CLUSTER): + { + for (i=0; i= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + + /* repair clusters */ + if (!REPAIR_CLUSTERS) for (i=0; i= 2)) + repair_cluster(i, particle, cluster, 1000, 0); + + printf("%i collisions\n", ncollisions); + delta_n = ncollisions - oldncollisions; + printf("delta_n = %i\n", delta_n); + return(ncollisions); + }case (CHEM_POLYGON_ONECLUSTER): + { + i = 0; + reacted = 0; + while ((i oldncollisions) reacted = 1; + } + + /* decouple particles with several partners from thermostat */ + if (particle[i].npartners >= AGREG_DECOUPLE) particle[i].thermostat = 0; + } + i++; + } + + /* repair clusters */ + if (!REPAIR_CLUSTERS) for (i=0; i= 2)) + repair_cluster(i, particle, cluster, 1000, 1); + printf("%i collisions\n", ncollisions); delta_n = ncollisions - oldncollisions; printf("delta_n = %i\n", delta_n); @@ -12615,19 +14583,110 @@ void reset_energy(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], doub } -void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_position, int traj_length, +int init_cluster_config(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES]) +/* initialize the clusters for option CLUSTER_PARTICLES */ +/* returns number of active clusters */ +{ + int i, j, tmp, nclusters = 0; + + for (i=0; i cluster %i -> particle %i\n", i, particle[i].cluster, cluster[particle[i].cluster].particle[0]); +// } + + return(nclusters); +} + + + +void compute_cluster_force(t_cluster cluster[NMAXCIRCLES], t_particle particle[NMAXCIRCLES]) +/* compute force and torque on clusters from force and torque on their particles */ +{ + int i, p, j; + double fx, fy, torque, energy; + + for (i=0; i= DPI) angle -= DPI; rde[i*NY+j].field_arg = angle; @@ -2832,7 +3130,7 @@ void compute_field_argument(double *phi[NFIELDS], t_rde rde[NX*NY]) for (i=0; i= DPI) arg -= DPI; rde[i*NY+j].field_arg = arg; @@ -2953,6 +3251,7 @@ void adjust_height(double *phi[NFIELDS], t_rde rde[NX*NY]) { value = VSCALE_WATER_HEIGHT*(phi[0][i*NY+j] + ADD_HEIGHT_CONSTANT); rde[i*NY+j].height = value; +// printf("[adjust_height] value = %.3lg\n", value); } } @@ -2961,6 +3260,14 @@ void compute_direction(double *phi[NFIELDS], t_rde rde[NX*NY]) { int i, j; double value; +// static double phaseshift; +// static int first = 1; +// +// if (first) +// { +// phaseshift = PHASE_SHIFT*PID/90.0; +// first = 0; +// } #pragma omp parallel for private(i,j,value) for (i=0; i= DPI) angle -= DPI; // phi_arg[i*NY+j] = 360.0*angle/DPI; @@ -4852,7 +5159,7 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3] case (Z_POLAR): { // hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette); - value += COLOR_PHASE_SHIFT*PI; + value += PHASE_SHIFT*PI; if (value > DPI) value -= DPI; color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb); break; @@ -5011,7 +5318,7 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3] case (Z_SWATER_DIRECTION_SPEED): { value *= 360.0/DPI; - value += 180.0*COLOR_PHASE_SHIFT; + value += 180.0*PHASE_SHIFT; if (value > 360.0) value -= 360.0; hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette); break; @@ -5163,6 +5470,7 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, { int i, j; +// printf("[compute_rde_fields] zplot = %i\n", zplot); switch (RDE_EQUATION) { case (E_RPS): { @@ -5233,11 +5541,13 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, } case (E_SHALLOW_WATER): { +// printf("[compute_rde_fields]\n"); if (zplot == Z_SWATER_HEIGHT) adjust_height(phi, rde); if ((zplot == Z_SWATER_SPEED)||(cplot == Z_SWATER_SPEED)||(zplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)) compute_speed(phi, rde, wsphere); if ((zplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)) compute_direction(phi, rde); + break; } default : break; } @@ -5605,6 +5915,7 @@ void init_cfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int cplot, t_ { #pragma omp parallel for private(i,j) for (i=0; i= NX-1) imin = NX-2; if (imax >= NX-1) imax = NX-2; - jmax = NY-3; + jmax = NY-POLE_NODRAW; + jmin = POLE_NODRAW; // jmax = NY - 50; jmid = (int)((double)NY*(observer_latitude + PID)/PI); imid = (imin + imax)/2; @@ -5910,12 +6225,12 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX { for (i=imax; i>imid; i--) { - for (j=0; j<=jmax; j++) + for (j=jmin; j<=jmax; j++) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imid; i>imin; i--) { - for (j=0; j<=jmid; j++) + for (j=jmin; j<=jmid; j++) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=jmax; j>=jmid; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); @@ -5923,16 +6238,16 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX for (i=imax+1; i=jmid; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); @@ -5950,13 +6265,13 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX { for (i=imax; i=jmid; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); @@ -5964,16 +6279,16 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX for (i=imax-1; i>=0; i--) { - for (j=0; j<=jmax; j++) + for (j=jmin; j<=jmax; j++) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } - for (j=0; j<=jmax; j++) + for (j=jmin; j<=jmax; j++) draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (i=NX-2; i>=imin; i--) { - for (j=0; j<=jmid; j++) + for (j=jmin; j<=jmid; j++) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (j=jmax; j>=jmid; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); @@ -6006,7 +6321,7 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX { for (i=imax; i>imid; i--) { - for (j=jmax; j>=0; j--) + for (j=jmax; j>=jmin; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } @@ -6014,59 +6329,102 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX { for (j=jmax; j>=jmid; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); - for (j=0; j<=jmid; j++) + for (j=jmin; j<=jmid; j++) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imax+1; i=0; j--) + for (j=jmax; j>=jmin; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } - for (j=jmax; j>=0; j--) + for (j=jmax; j>=jmin; j--) draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); for (i=0; i<=imin; i++) { for (j=jmax; j>=jmid; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); - for (j=0; j<=jmid; j++) + for (j=jmin; j<=jmid; j++) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } + + /* TEST */ + for (j=jmin; j<=jmid; j++) + for (i=-1; i<2; i++) if ((imin+i >= 0)&&(imin+i+1 < NX)) + draw_wave_sphere_rde_ij(imin+i, imin+i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + for (j=jmax; j>=jmid; j--) + for (i=-1; i<2; i++) if ((imin+i >= 0)&&(imin+i+1 < NX)) + draw_wave_sphere_rde_ij(imin+i, imin+i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + +// for (i=0; i<2; i++) if (imin >= i) +// { +// for (j=2*jmax/3; j>=jmid; j--) +// draw_wave_sphere_rde_ij(imin-i, imin-i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); +// for (j=jmin; j<=jmid; j++) +// draw_wave_sphere_rde_ij(imin-i, imin-i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); +// } + } else { for (i=imax; i=0; j--) + for (j=jmax; j>=jmin; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } - for (i=imid; i=jmid; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); - for (j=0; j<=jmid; j++) + for (j=jmin; j<=jmid; j++) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } for (i=imax-1; i>=0; i--) { - for (j=jmax; j>=0; j--) + for (j=jmax; j>=jmin; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } - for (j=jmax; j>=0; j--) + for (j=jmax; j>=jmin; j--) draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); - for (i=NX-2; i>=imin-1; i--) + for (i=NX-2; i>=imin; i--) { for (j=jmax; j>=jmid; j--) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); - for (j=0; j<=jmid; j++) + for (j=jmin; j<=jmid; j++) draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); } + + for (i=0; i<2; i++) if (imin >= i) + { + for (j=2*jmax/3; j>=jmid; j--) + draw_wave_sphere_rde_ij(imin-i, imin-i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + for (j=jmin; j<=jmid; j++) + draw_wave_sphere_rde_ij(imin-i, imin-i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + } + + /* TEST */ + for (j=jmin; j<=jmid; j++) + for (i=-1; i<2; i++) if ((imin+i >= 0)&&(imin+i+1 < NX)) + draw_wave_sphere_rde_ij(imin+i, imin+i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + for (j=jmax; j>=jmid; j--) + for (i=-1; i<2; i++) if ((imin+i >= 0)&&(imin+i+1 < NX)) + draw_wave_sphere_rde_ij(imin+i, imin+i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); + +// if (imin >= 1) +// { +// /* experimental */ +// for (j=2*jmax/3; j>=jmid; j--) +// draw_wave_sphere_rde_ij(imin-1, imin, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); +// for (j=0; j<=jmid; j++) +// draw_wave_sphere_rde_ij(imin-1, imin, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value); +// } + } /* South pole */ @@ -6329,9 +6687,9 @@ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rd if (refresh) { -// printf("Computing fields\n"); + printf("Computing fields\n"); compute_rde_fields(phi, xy_in, zplot, cplot, rde, wsphere); -// printf("Computed fields\n"); + printf("Computed fields\n"); if ((PLOT_3D)&&(SHADE_3D)) { if (SPHERE) diff --git a/sub_wave.c b/sub_wave.c index e7957f3..d8a65c7 100644 --- a/sub_wave.c +++ b/sub_wave.c @@ -685,12 +685,74 @@ double michelson_schedule(int i) } +int generate_poisson_discs(t_circle circles[NMAXCIRCLES], double xmin, double xmax, double ymin, double ymax, double dpoisson) +/* generate a Poisson disc sample in a given rectangle */ +{ + int i, j, k, n_p_active, ncandidates=5000, naccepted; + double r, phi, x, y; + short int active_poisson[NMAXCIRCLES], far; + + printf("Generating Poisson disc sample\n"); + /* generate first circle */ + circles[0].xc = (xmax - xmin)*(double)rand()/RAND_MAX + xmin; + circles[0].yc = (ymax - ymin)*(double)rand()/RAND_MAX + ymin; + active_poisson[0] = 1; + n_p_active = 1; + ncircles = 1; + + while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES)) + { + /* randomly select an active circle */ + i = rand()%(ncircles); + while (!active_poisson[i]) i = rand()%(ncircles); + /* 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 circle at (%.3f,%.3f) accepted\n", x, y); + circles[ncircles].xc = x; + circles[ncircles].yc = y; + circles[ncircles].radius = MU; + circles[ncircles].active = 1; + active_poisson[ncircles] = 1; + ncircles++; + n_p_active++; + naccepted++; + } + } + if (naccepted == 0) /* inactivate circle i */ + { + active_poisson[i] = 0; + n_p_active--; + } + printf("%i active circles\n", n_p_active); + } + + printf("Generated %i circles\n", ncircles); + return(ncircles); +} + + int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern) /* initialise the arrays circlex, circley, circlerad and circleactive */ /* for billiard shape D_CIRCLES */ { int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted; - double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = 3.25*MU, xx[4], yy[4], dr, dphi; + double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = PDISC_FACTOR*MU, xx[4], yy[4], dr, dphi; short int active_poisson[NMAXCIRCLES], far; switch (circle_pattern) { @@ -737,7 +799,8 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern { n = NGRIDY*i + j; circles[n].xc = ((double)(i-NGRIDX/2) + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; - circles[n].yc = YMIN + 0.5 + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; + circles[n].yc = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; +// circles[n].yc = YMIN + 0.5 + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy; circles[n].radius = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5)); circles[n].active = 1; } @@ -955,7 +1018,7 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern r0 = 2.0*MU; r = r0 + MU; - for (i=0; i<1000; i++) + for (i=0; i<2000; i++) { x = r*cos(phi); y = r*sin(phi); @@ -1003,13 +1066,13 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern { ncircles = NGRIDX*NGRIDY; dphi = DPI/((double)NGRIDX); - dr = 0.5*LAMBDA/(double)NGRIDY; + dr = (1.0 - RADIUS_FACTOR)*LAMBDA/(double)NGRIDY; for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; phi = (double)i*dphi; - r = 0.5*LAMBDA + (double)j*dr; + r = RADIUS_FACTOR*LAMBDA + (double)j*dr; circles[n].xc = r*cos(phi); circles[n].yc = r*sin(phi); circles[n].radius = MU; @@ -1023,14 +1086,14 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern { ncircles = NGRIDX*NGRIDY; dphi = DPI/((double)NGRIDX); - dr = 0.5*LAMBDA/(double)NGRIDY; + dr = (1.0 - RADIUS_FACTOR)*LAMBDA/(double)NGRIDY; for (i = 0; i < NGRIDX; i++) for (j = 0; j < NGRIDY; j++) { n = NGRIDY*i + j; phi = (double)i*dphi; phi += 0.5*(double)j*dphi; - r = 0.5*LAMBDA + (double)j*dr; + r = RADIUS_FACTOR*LAMBDA + (double)j*dr; circles[n].xc = r*cos(phi); circles[n].yc = r*sin(phi); circles[n].radius = MU; @@ -1048,10 +1111,11 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern gamma = (sqrt(5.0) - 1.0)*PI; /* golden mean times 2Pi */ phi = 0.0; - r0 = 0.5*LAMBDA; +// r0 = 0.5*LAMBDA; + r0 = RADIUS_FACTOR*LAMBDA; r = r0 + MU; - for (i=0; i<1000; i++) + for (i=0; i<(int)(1000.0/RADIUS_FACTOR); i++) { x = r*cos(phi); y = r*sin(phi); @@ -1075,6 +1139,40 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern } break; } + case (C_RINGS_POISSONDISC): + { + ncircles = generate_poisson_discs(circles, YMIN, YMAX, YMIN, YMAX, PDISC_FACTOR*MU); + for (i=0; i RADIUS_FACTOR*LAMBDA)) circles[i].active = 1; + else circles[i].active = 0; + } + break; + } + case (C_RINGS_LOGSPIRAL): + { + ncircles = NGRIDX*NGRIDY; + dr = (1.0 - RADIUS_FACTOR)*LAMBDA/(double)NGRIDY; + dphi = DPI/((double)NGRIDX); + for (i = 0; i < NGRIDX; i++) + for (j = 0; j < NGRIDY; j++) + { + n = NGRIDY*i + j; + r = RADIUS_FACTOR*LAMBDA + (double)j*dr; + phi = (double)i*dphi; + phi += log(r/(RADIUS_FACTOR*LAMBDA)); + circles[n].xc = r*cos(phi); + circles[n].yc = r*sin(phi); + circles[n].radius = MU; + /* activate only circles that intersect the domain */ + if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1; + else circles[n].active = 0; + } + break; + } case (C_HEX_BOTTOM): { ncircles = NGRIDX*(NGRIDY+1); @@ -2545,6 +2643,60 @@ int rc_hyp(double x, double y) } } + +int init_xyin_from_image(short int * xy_in[NX]) +/* initialize table xy_in from an image file */ +{ + FILE *image_file; + int nx, ny, maxrgb, i, j, k, ii, jj, nmaxpixels = 2000000, rgbtot, scan, rgbval; + int *rgb_values; + double scalex, scaley; + + image_file = fopen("PHOTONSband.ppm", "r"); + scan = fscanf(image_file,"%i %i\n", &nx, &ny); + scan = fscanf(image_file,"%i\n", &maxrgb); + + rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int)); + + scalex = (double)nx/(double)NX; + scaley = (double)ny/(double)NY; + + if (nx*ny > nmaxpixels) + { + printf("Image too large, increase nmaxpixels in init_xyin_from_image()\n"); + 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 > 3*maxrgb/2) xy_in[i][j] = 1; + else xy_in[i][j] = 0; + } + + fclose(image_file); + free(rgb_values); + return(1); +} + + int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_circle *circles) /* returns 1 if (x,y) represents a point in the billiard */ { @@ -3547,6 +3699,52 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ /* use IOR_GRADIENT_INDEX_LENS for index of refraction control */ return(1); } + case (D_IMAGE): + { + /* Not needed, uses option XYIN_INITIALISED */ + return(1); + } + case (D_MAGNETRON): + { + r = module2(x,y); + if (r > LAMBDA) return(1); + if (r < 0.33*LAMBDA) return(1); + if ((vabs(y) < WALL_WIDTH)&&(x > -0.66*LAMBDA)) return(1); + if ((vabs(x) < WALL_WIDTH)&&(vabs(y) < 0.66*LAMBDA)) return(1); + if ((vabs(x-y) < WALL_WIDTH*sqrt(2.0))&&(r < 0.66*LAMBDA)) return(1); + if ((vabs(x+y) < WALL_WIDTH*sqrt(2.0))&&(r < 0.66*LAMBDA)) return(1); + for (k=0; k<8; k++) + { + x1 = 0.66*cos((double)k*0.5*PID); + y1 = 0.66*sin((double)k*0.5*PID); + r = module2(x - x1, y - y1); + if (r < MU) return(1); + } + + + return(0); + } + case (D_MAGNETRON_CATHODE): + { + r = module2(x,y); + if (r < WALL_WIDTH) return(0); + if (r > LAMBDA) return(1); + if (r < 0.33*LAMBDA) return(1); + if ((vabs(y) < WALL_WIDTH)&&(x > -0.66*LAMBDA)) return(1); + if ((vabs(x) < WALL_WIDTH)&&(vabs(y) < 0.66*LAMBDA)) return(1); + if ((vabs(x-y) < WALL_WIDTH*sqrt(2.0))&&(r < 0.66*LAMBDA)) return(1); + if ((vabs(x+y) < WALL_WIDTH*sqrt(2.0))&&(r < 0.66*LAMBDA)) return(1); + for (k=0; k<8; k++) + { + x1 = 0.66*cos((double)k*0.5*PID); + y1 = 0.66*sin((double)k*0.5*PID); + r = module2(x - x1, y - y1); + if (r < MU) return(1); + } + + + return(0); + } case (D_MENGER): { x1 = 0.5*(x+1.0); @@ -5630,6 +5828,21 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound draw_line(WAVE_PROFILE_X, YMIN, WAVE_PROFILE_X, YMAX); break; } + case (D_IMAGE): + { + /* do nothing */ + break; + } + case (D_MAGNETRON): + { + /* TODO */ + break; + } + case (D_MAGNETRON_CATHODE): + { + /* TODO */ + break; + } case (D_MENGER): { glLineWidth(3); @@ -6337,7 +6550,7 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2, void draw_circular_color_scheme_palette_fade(double x1, double y1, double radius, int plot, double min, double max, int palette, int fade, double fade_value) { int j, k; - double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2], zscale = 0.9; + double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2], zscale = 0.85; // glBegin(GL_TRIANGLE_FAN); xy_to_pos(x1, y1, xy); @@ -6465,6 +6678,13 @@ void draw_circular_color_scheme_palette_fade(double x1, double y1, double radius color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb); break; } + case (Z_SWATER_DIRECTION_SPEED): + { + value = dy_phase*(double)(j) - 0.5*PHASE_SHIFT + 0.5; + if (value > 1.0) value -= 1.0; + color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); + break; + } default: { value = min + 1.0*dy*(double)(j); diff --git a/sub_wave_3d_rde.c b/sub_wave_3d_rde.c index b8e35fd..c5c2b93 100644 --- a/sub_wave_3d_rde.c +++ b/sub_wave_3d_rde.c @@ -604,6 +604,145 @@ void init_earth_map_rde(t_wave_sphere *wsphere, int res) } } + +void init_planet_map_rde(t_wave_sphere wsphere[NX*NY], int planet, int res) +/* init file from planetary map */ +{ + int i, j, ii, jj, k, nx, ny, maxrgb, nmaxpixels, scan, rgbval, diff, sshift, nshift, ishift, dem_number; + int *rgb_values; + double cratio, rx, ry, cy, vshift; + FILE *image_file; + + switch (planet){ + case (D_SPHERE_MARS): + { + printf("Reading Mars map\n"); + nmaxpixels = 8388608; + image_file = fopen("Mars_Viking_ClrMosaic_global_925m_scaled.ppm", "r"); + dem_number = DEM_MARS; + break; + } + case (D_SPHERE_MOON): + { + printf("Reading Moon map\n"); + nmaxpixels = 2048*1024; + image_file = fopen("Moon_photo_map.ppm", "r"); + dem_number = DEM_MOON; + break; + } + case (D_SPHERE_VENUS): + { + printf("Reading Venus map\n"); + nmaxpixels = 1440*720; + image_file = fopen("Venus_map_NASA_JPL_Magellan-Venera-Pioneer.ppm", "r"); + dem_number = DEM_VENUS; + break; + } + case (D_SPHERE_MERCURY): + { + printf("Reading Mercury map\n"); + nmaxpixels = 2304*1152; + image_file = fopen("Mercury_color_photo.ppm", "r"); + dem_number = DEM_MERCURY; + break; + } + } + + scan = fscanf(image_file,"%i %i\n", &nx, &ny); + scan = fscanf(image_file,"%i\n", &maxrgb); + printf("nx*ny = %i\n", nx*ny); + rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int)); + + if (nx*ny > nmaxpixels) + { + printf("Image too large, increase nmaxpixels in init_planet_map()\n"); + exit(0); + } + + /* shift due to min/max latitudes of image */ + sshift = 0 + DPOLE; + nshift = 0 + DPOLE; + + /* choice of zero meridian */ + ishift = (int)(nx*ZERO_MERIDIAN/360.0); + + /* read rgb values */ + for (j=0; j nx-1) ii -= nx; +// jj = (int)(-ry*(double)j + cy); +// jj = (int)(ry*(double)(NY-nshift - j)) + sshift; + jj = (int)(ry*(double)(res*NY-nshift - j)); + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + wsphere[i*res*NY+j].r = (double)rgb_values[3*(jj*nx+ii)]*cratio; + wsphere[i*res*NY+j].g = (double)rgb_values[3*(jj*nx+ii)+1]*cratio; + wsphere[i*res*NY+j].b = (double)rgb_values[3*(jj*nx+ii)+2]*cratio; + +// printf("RGB at (%i, %i) = (%.3lg, %3.lg, %.3lg)\n", i, j, wsphere[i*NY+j].r, wsphere[i*NY+j].g, wsphere[i*NY+j].b); + + /* decide which points are in the Sea */ + wsphere[i*NY+j].indomain = 1; + wsphere[i*NY+j].draw_wave = 1; + } + + + /* smooth colors in case of high resolution */ + if ((res*NX > nx)||(res*NY > ny)) + for (i=1; i 1.0) value -= 1.0; + if (value < 0.0) value += 1.0; color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb); break; } diff --git a/wave_3d.c b/wave_3d.c index 7a8ba4c..2ba99ca 100644 --- a/wave_3d.c +++ b/wave_3d.c @@ -43,7 +43,7 @@ #include #include -#define MOVIE 0 /* set to 1 to generate movie */ +#define MOVIE 1 /* set to 1 to generate movie */ #define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 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 */ @@ -83,6 +83,7 @@ #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 RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ #define LAMBDA 1.0 /* parameter controlling the dimensions of domain */ @@ -97,6 +98,7 @@ #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 RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -227,6 +229,7 @@ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define VSCALE_AMPLITUDE 2.0 /* 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 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 */ @@ -287,6 +290,7 @@ #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 XYIN_INITIALISED (B_DOMAIN == D_IMAGE) /* end of constants only used by sub_wave and sub_maze */ diff --git a/wave_billiard.c b/wave_billiard.c index 475ee82..f31107f 100644 --- a/wave_billiard.c +++ b/wave_billiard.c @@ -57,8 +57,6 @@ #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 3840 /* number of grid points on x axis */ #define NY 2300 /* number of grid points on y axis */ @@ -74,20 +72,21 @@ /* Choice of the billiard table */ -#define B_DOMAIN 20 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 6 /* choice of domain shape, see list in global_pdes.c */ -#define CIRCLE_PATTERN 1 /* pattern of circles or polygons, see list in global_pdes.c */ +#define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */ #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.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ -#define NPOISSON 1000 /* number of points for Poisson C_RAND_POISSON arrangement */ +#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 RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ -#define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ -#define MU 0.012 /* parameter controlling the dimensions of domain */ +#define LAMBDA 0.3 /* parameter controlling the dimensions of domain */ +#define MU 1.0 /* parameter controlling the dimensions of domain */ #define MU_B 1.0 /* parameter controlling the dimensions of domain */ #define NPOLY 6 /* number of sides of polygon */ #define APOLY -0.666666666666 /* angle by which to turn polygon, in units of Pi/2 */ @@ -96,9 +95,10 @@ #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 24 /* number of grid point for grid of disks */ -#define NGRIDY 40 /* number of grid point for grid of disks */ -#define WALL_WIDTH 0.1 /* width of wall separating lenses */ +#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 RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -116,7 +116,7 @@ /* Physical parameters of wave equation */ -#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ #define OSCILLATE_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ #define OSCILLATION_SCHEDULE 61 /* oscillation schedule, see list in global_pdes.c */ @@ -128,8 +128,10 @@ #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.2 /* Courant number */ -#define COURANTB 0.0 /* Courant number in medium B */ +#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 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 */ @@ -144,16 +146,16 @@ /* 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 8 /* period of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 3 /* period of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ -#define N_SOURCES 2 /* number of sources, for option draw_sources */ +#define N_SOURCES 1 /* number of sources, for option draw_sources */ #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 */ #define N_WAVE_PACKETS 5 /* number of wave packets */ #define WAVE_PACKET_RADIUS 50 /* radius of wave packets */ -#define USE_INPUT_TIMESERIES 1 /* set to 1 to use a time series (Morse code) as input * / +#define USE_INPUT_TIMESERIES 0 /* set to 1 to use a time series (Morse code) as input * / /* Boundary conditions, see list in global_pdes.c */ @@ -161,9 +163,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 1600 /* 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 NSTEPS 1000 /* number of frames of movie */ +// #define NSTEPS 2000 /* number of frames of movie */ +#define NVID 10 /* number of iterations between images displayed on screen */ +#define NSEG 1000 /* number of segments of boundary */ #define INITIAL_TIME 0 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PRINT_SPEED 0 /* print speed of moving source */ @@ -179,7 +182,7 @@ /* Parameters of initial condition */ -#define INITIAL_AMP 1.5 /* amplitude of initial condition */ +#define INITIAL_AMP 0.75 /* amplitude of initial condition */ #define INITIAL_VARIANCE 0.00001 /* variance of initial condition */ #define INITIAL_WAVELENGTH 0.025 /* wavelength of initial condition */ @@ -192,7 +195,7 @@ /* Color schemes */ #define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 13 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -202,15 +205,19 @@ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ +// #define COLOR_PHASE_SHIFT 0.0 /* phase shift of color scheme, in units of Pi */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define E_SCALE 55.0 /* scaling factor for energy representation */ -#define LOG_SCALE 1.0 /* scaling factor for energy log representation */ -#define LOG_SHIFT 0.5 /* shift of colors on log scale */ -#define FLUX_SCALE 5.0e3 /* scaling factor for energy flux represtnation */ -#define AVRG_E_FACTOR 0.95 /* controls time window size in P_AVERAGE_ENERGY scheme */ +#define VSHIFT_AMPLITUDE 0.0 /* additional shift for wave amplitude */ +#define VSCALE_AMPLITUDE 0.5 /* additional scaling factor for wave amplitude */ +#define E_SCALE 50.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.75 /* 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 0 /* set to 1 to fade color inside obstacles */ +#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 */ +// #define SHADE_SCALE_2D 0.25 /* lower value increases sensitivity of shading */ #define SHADE_SCALE_2D 0.05 /* lower value increases sensitivity of shading */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -220,24 +227,24 @@ #define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ -#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 2.0 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 0.8 /* scale of color scheme bar for 2nd part */ +#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 1.4 /* scale of color scheme bar */ +#define COLORBAR_RANGE_B 1.2 /* 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 */ -#define DRAW_WAVE_PROFILE 1 /* set to 1 to draw a profile of the wave */ +#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define HORIZONTAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */ #define VERTICAL_WAVE_PROFILE 1 /* set to 1 to draw wave profile vertically */ #define WAVE_PROFILE_X 1.9 /* value of x to sample wave profile */ -#define WAVE_PROFILE_Y 0.075 /* value of y to sample wave profile */ +#define WAVE_PROFILE_Y 0.12 /* value of y to sample wave profile */ #define PROFILE_AT_BOTTOM 1 /* draw wave profile at bottom instead of top */ -#define AVERAGE_WAVE_PROFILE 1 /* set to 1 to draw time-average of wave profile squared*/ +#define AVERAGE_WAVE_PROFILE 0 /* set to 1 to draw time-average of wave profile squared*/ #define DRAW_WAVE_TIMESERIES 0 /* set to 1 to draw a time series of the wave, 2 to also draw it at the top */ #define TIMESERIES_NVALUES 400 /* number of values plotted in time series */ #define SAVE_TIME_SERIES 0 /* set to 1 to save wave time series at a point */ -#define DRAW_WAVE_SOURCE 1 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y), set to 2 to also draw focus */ +#define DRAW_WAVE_SOURCE 0 /* set to 1 to draw source of wave at (wave_source_x, wave_source_y), set to 2 to also draw focus */ #define MESSAGE_LDASH 14 /* length of dash for Morse code message */ #define MESSAGE_LDOT 8 /* length of dot for Morse code message */ @@ -265,6 +272,7 @@ #define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX) #define REFRESH_IOR ((IOR == IOR_PERIODIC_WELLS_ROTATING)||(IOR == IOR_PERIODIC_WELLS_ROTATING_LARGE)) +#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE) double light[2] = {0.40824829, 0.816496581}; /* location of light source for SHADE_2D option*/ @@ -588,10 +596,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 = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle = 0.0, x1, sign1 = 1.0, ior_angle = 0.0, omega, phase_shift, vshift, dsource, finv; + 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 *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]; + 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; // static int image_counter = 0; int image_counter = 0; long int wave_value; @@ -628,6 +636,8 @@ void animation() init_wave_packets(packet, WAVE_PACKET_RADIUS); } + for (i=0; i vmax[pnb]) vmax[pnb] = average_vals[i]; + if (average_vals[pnb*NX+i] > vmax[pnb]) vmax[pnb] = average_vals[pnb*NX+i]; } - else /*if (!fade)*/ + else if (!fade) { for (i=imin; i 0.0)) - { - if (vmax[pnb] > -vmin[pnb]) vmin[pnb] = -vmax[pnb]; - else if (vmax[pnb] < -vmin[pnb]) vmax[pnb] = -vmin[pnb]; - } +// if ((vmin[pnb] < 0.0)&&(vmax[pnb] > 0.0)) +// { +// if (vmax[pnb] > -vmin[pnb]) vmin[pnb] = -vmax[pnb]; +// else if (vmax[pnb] < -vmin[pnb]) vmax[pnb] = -vmin[pnb]; +// } // printf("vmin = %.3lg, vmax = %.3lg\n", vmin, vmax); deltav = vmax[pnb]-vmin[pnb]; if (deltav <= 0.0) deltav = 0.01; a = deltaj/deltav; b = ymin - a*vmin[pnb]; - + erase_area_ij(imin, jmin, imax, jmax); if (fade) glColor3f(fade_value, fade_value, fade_value); @@ -1232,7 +1248,7 @@ void draw_wave_profile_horizontal(double *values, int size, int average, int pnb for (i=imin; i ymin) glVertex2d((double)i, y); } @@ -1250,6 +1266,8 @@ void draw_wave_profile_horizontal(double *values, int size, int average, int pnb { vmax[pnb] *= 0.99; vmin[pnb] *= 0.99; +// vmax[pnb] *= 0.95; +// vmin[pnb] *= 0.95; } } @@ -1260,7 +1278,7 @@ void draw_wave_profile_vertical(double *values, int size, int average, int pnb, int j, k, ij[2]; double deltav, a, b, x; static int imin, imax, ival, jmin, jmax, imid, d, first = 1; - static double vmin[2], vmax[2], deltai, xmin, average_vals[NY]; + static double vmin[2], vmax[2], deltai, xmin, average_vals[2*NY]; if (first) { @@ -1276,7 +1294,7 @@ void draw_wave_profile_vertical(double *values, int size, int average, int pnb, d = (imax - imin)/10 + 1; deltai = (double)(imax - imin - 2*d); xmin = (double)(imin + d); - for (j=0; j vmax[pnb]) vmax[pnb] = average_vals[j]; + if (average_vals[pnb*NY+j] > vmax[pnb]) vmax[pnb] = average_vals[pnb*NY+j]; } @@ -1345,7 +1363,7 @@ void draw_wave_profile_vertical(double *values, int size, int average, int pnb, glBegin(GL_LINE_STRIP); for (j=jmin; j= xmin) glVertex2d(x, (double)j); @@ -1722,25 +1740,33 @@ void init_fade_table(double *tcc_table[NX], double *fade_table) for (i=0; i