From b0dbe991d1f747eb9862409e3bce9c94d8cfa9c2 Mon Sep 17 00:00:00 2001 From: Nils Berglund <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Tue, 30 Dec 2025 18:32:47 +0100 Subject: [PATCH] Add files via upload Bugs corrected in lennardjones.c New domains for wave_billiard.c --- global_ljones.c | 37 ++- global_pdes.c | 12 +- lennardjones.c | 303 ++++++++++++------- sub_lj.c | 777 ++++++++++++++++++++++++++++++++++++++++------- sub_lj_sphere.c | 340 +++++++++++++++++++-- sub_wave.c | 790 +++++++++++++++++++++++++++++++++++++++++++++++- sub_wave_3d.c | 29 +- wave_billiard.c | 81 ++--- 8 files changed, 2068 insertions(+), 301 deletions(-) diff --git a/global_ljones.c b/global_ljones.c index 9225759..886c8f6 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -26,6 +26,7 @@ #define NMAXSHOVELS 50 /* max number of shovels */ #define NMAX_TRIANGLES_PER_OBSTACLE 10 /* max number of triangles per obstacle */ #define NMAX_TRIANGLES_PER_FACET 4 /* max number of triangles per facet between obstacles */ +#define NMAX_ABSORBERS 10 /* max number of absorbing discs */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -36,6 +37,7 @@ #define C_CLOAK_A 6 /* first optimized invisibility cloak */ #define C_LASER 7 /* laser fight in a room of mirrors */ #define C_POISSON_DISC 8 /* Poisson disc sampling */ +#define C_POISSON_DISC_SPHERE 81 /* Poisson disc sampling in spherical geometry */ #define C_GOLDEN_MEAN 10 /* pattern based on vertical shifts by golden mean */ #define C_GOLDEN_SPIRAL 11 /* spiral pattern based on golden mean */ @@ -43,6 +45,8 @@ #define C_POOL_TABLE 20 /* pool table initial position */ +#define C_SQUARE_SPHERE 30 /* "square" grid of circles on a sphere */ + #define C_ONE 97 /* one single circle, as for Sinai */ #define C_TWO 98 /* two concentric circles of different type */ #define C_NOTHING 99 /* no circle at all, for comparisons */ @@ -83,6 +87,7 @@ /* pattern of additional repelling segments */ #define S_RECTANGLE 0 /* segments forming a rectangle */ #define S_CYLINDRICAL 100 /* two lines at top and bottom */ +#define S_BOTTOM 110 /* one line at bottom */ #define S_CUP 1 /* segments forming a cup (for increasing gravity) */ #define S_HOURGLASS 2 /* segments forming an hour glass */ #define S_PENTA 3 /* segments forming a pentagon with 3 angles of 120° and 2 right angles */ @@ -252,6 +257,7 @@ #define CHEM_H2O_H_OH 20 /* H2O <-> H+ + OH- */ #define CHEM_2H2O_H3O_OH 21 /* 2 H2O <-> H3O+ + OH- */ #define CHEM_AGGREGATION 22 /* agregation of molecules coming close */ +#define CHEM_AGGREGATION_MAX 221 /* agregation of molecules coming close, with max total partner number */ #define CHEM_AGGREGATION_CHARGE 23 /* agregation of charged molecules coming close */ #define CHEM_AGGREGATION_CHARGE_NOTRIANGLE 231 /* agregation of charged molecules coming close, no multiple pairings */ #define CHEM_AGGREGATION_NNEIGH 24 /* agregation of molecules with limitation on neighbours */ @@ -409,6 +415,14 @@ #define VP_ORBIT2 11 /* rotate in a plane specified by max latitude */ #define VP_POLAR 2 /* polar orbit */ +/* Absorber pattern */ + +#define AP_ONE 0 /* one single absorber */ +#define AP_TWO 1 /* two absorbers, symmetric wrp x-axis */ +#define AP_CUBE 2 /* absorbers on the vertices of a cube */ +#define AP_DODECA 3 /* absorbers on the vertices of a dodecahedron */ +#define AP_ICOSA 4 /* absorbers on the vertices of an icosahedron */ + /* Color schemes */ #define C_LUM 0 /* color scheme modifies luminosity (with slow drift of hue) */ @@ -571,6 +585,12 @@ typedef struct short int chessboard; /* has value 1 on a chessboard, for some arrangements */ } t_obstacle; +typedef struct +{ + double xc, yc; /* center of absorber */ + double radius; /* radius of absorber */ +} t_absorber; + typedef struct { int i[3]; /* indices of obstacles forming triangle */ @@ -678,6 +698,7 @@ typedef struct typedef struct { int nactive; /* number of active particles */ + int nabsorbed; /* number of absorbed particles */ double beta; /* inverse temperature */ double mean_energy; /* mean energy */ double krepel; /* force constant */ @@ -704,24 +725,12 @@ typedef struct double reg_cottheta; /* regularized cotangent of theta */ double x, y, z; /* x, y, z coordinates of point on sphere */ double radius; /* radius with wave height */ -// double radius_dem; /* radius with digital elevation model */ double r, g, b; /* RGB values for image */ -// short int indomain; /* has value 1 if lattice point is in domain */ -// short int draw_wave; /* has value 1 if wave instead of DEM is drawn */ -// short int evolve_wave; /* has value 1 where there is wave evolution */ -// double x2d, y2d; /* x and y coordinates for 2D representation */ -// double altitude; /* altitude in case of Earth with digital elevation model */ double cos_angle; /* cosine of light angle */ double cos_angle_sphere; /* cosine of light angle for perfect sphere */ -// double force; /* external forcing */ -// double phigrid, thetagrid; /* phi, theta angles for alt simulation grid on sphere */ -// short int nneighb; /* number of neighbours, for Kuramoto model on sphere */ -// int neighbor[NMAX_SPHERE_NEIGHB]; /* list of neighbours */ -// int convert_grid; /* convert field from simulation grid to longitude-latitude */ -// short int edge; /* has value 1 on edges of cubic simulation grid */ -// double cos_grid, sin_grid; /* cosine and sine of grid angle */ + short int locked; /* has value 1 if color of pixel is fixed */ } t_lj_sphere; -int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0, n_otriangles = 0, n_ofacets = 0; +int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0, n_otriangles = 0, n_ofacets = 0, nabsorbers = 0; FILE *lj_log; diff --git a/global_pdes.c b/global_pdes.c index be6d69b..65b956d 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -41,6 +41,11 @@ #define D_TWO_PARABOLAS 19 /* two facing parabolic antennas */ #define D_TWO_PARABOLAS_ASYM 191 /* confocal facing parabolas of different size */ #define D_TWO_PARABOLAS_ASYM_GUIDE 192 /* confocal facing parabolas of different size with a wave guide */ +#define D_TWO_PARABOLAS_SEMITRANS 193 /* 192 with a semitransparent part */ +#define D_TWO_PARABOLAS_LENS 194 /* 193 with a lens */ +#define D_TWO_PARABOLAS_CLOSED_SEMITRANS 195 /* cavity with two parabolic sides and semitransparent part */ +#define D_TWO_PARABOLAS_CLOSED_LENS 1951 /* 105 with a lens as window */ +#define D_TWO_PARABOLAS_CLOSED_LENS_R 1952 /* 105 with a lens of adaptable radius as window */ #define D_CIRCLES 20 /* several circles */ #define D_CIRCLES_IN_RECT 201 /* several circles in a rectangle */ @@ -57,6 +62,7 @@ #define D_POLYGONS 40 /* several polygons */ #define D_VONKOCH 41 /* von Koch snowflake fractal */ #define D_STAR 42 /* star shape */ +#define D_STAR_CHANNEL 421 /* star shape with an attached channel */ #define D_FRESNEL 43 /* Fresnel lens */ #define D_NOISEPANEL 44 /* zigzag noise insulating panel */ #define D_NOISEPANEL_RECT 441 /* comparison between zigzag noise insulating panel and flat walls */ @@ -102,6 +108,7 @@ #define D_TESLA_FOUR 72 /* four Tesla valves */ #define D_TREE 73 /* Christmas tree, to use with IOR_TREE */ +#define D_TREE_BEAM 731 /* Christmas tree with an attached wave guide */ #define D_MICHELSON 74 /* Michelson interferometer, to use with IOR_MICHELSON */ #define D_MICHELSON_MOVING 741 /* moving Michelson interferometer, to use with IOR_MICHELSON */ #define D_RITCHEY_CHRETIEN_SPHERICAL 75 /* Ritchey-Chrétien telescope with spherical mirrors */ @@ -136,7 +143,9 @@ #define D_CIRCLE_LATTICE_RHOMBUS 97 /* rhombus-based lattice of connected circles */ #define D_DISC_WAVEGUIDE 98 /* a disc with an attached waveguide */ #define D_DISC_WAVEGUIDE_SHIFTED 981 /* a disc with a shifted attached waveguide */ -#define D_DISC_WAVEGUIDE_ELLIPSE 982 /* a disc with a shifted attached waveguide */ +#define D_DISC_WAVEGUIDE_ELLIPSE 982 /* an ellipse with a shifted attached waveguide */ +#define D_WAVEGUIDE_ELLIPSE_OBLIQUE 983 /* an ellipse with an angled attached waveguide */ +#define D_CIRCLE_PAIRS 99 /* pairs of circles connected by channels */ /* for wave_sphere.c */ @@ -163,6 +172,7 @@ #define IM_HAPPY_NEW_YEAR 3 /* "Happy New Year 2025" */ #define IM_DICKSON 4 /* Dickson fjord and others nearbys in Greenland */ #define IM_DICKSON_ZOOM 5 /* zoom on Dickson fjord in Greenland */ +#define IM_HAPPY_NEW_YEAR_TWOSIX 6 /* "Happy New Year 2026" */ #define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */ diff --git a/lennardjones.c b/lennardjones.c index 9f23c6a..4130157 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -58,34 +58,34 @@ #define YMIN 0.0 #define YMAX 3.141592654 /* y interval for 9/16 aspect ratio */ -#define INITXMIN 0.025 -#define INITXMAX 6.255 /* x interval for initial condition */ -#define INITYMIN 0.4 -#define INITYMAX 2.74 /* y interval for initial condition */ +#define INITXMIN 2.95 +#define INITXMAX 3.4 /* x interval for initial condition */ +#define INITYMIN 1.37 +#define INITYMAX 1.77 /* y interval for initial condition */ #define THERMOXMIN -1.25 #define THERMOXMAX 1.25 /* x interval for initial condition */ #define THERMOYMIN 0.0 #define THERMOYMAX 0.75 /* y interval for initial condition */ -#define ADDXMIN -1.9 -#define ADDXMAX 1.9 /* x interval for adding particles */ -#define ADDYMIN 1.2 -#define ADDYMAX 1.3 /* y interval for adding particles */ +#define ADDXMIN 0.0 +#define ADDXMAX 0.1 /* x interval for adding particles */ +#define ADDYMIN 1.57 +#define ADDYMAX 1.57 /* y interval for adding particles */ #define ADDRMIN 2.0 #define ADDRMAX 2.1 /* r interval for adding particles */ #define BCXMIN 0.0 #define BCXMAX 6.283185307 /* x interval for boundary condition */ -#define BCYMIN 0.05 -#define BCYMAX 3.091592654 /* y interval for boundary condition */ +#define BCYMIN 0.3 +#define BCYMAX 2.841592654 /* y interval for boundary condition */ #define OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ #define OBSYMIN -1.125 #define OBSYMAX 1.125 /* x interval for motion of obstacle */ -#define CIRCLE_PATTERN 0 /* pattern of circles, see list in global_ljones.c */ +#define CIRCLE_PATTERN 20 /* pattern of circles, see list in global_ljones.c */ #define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */ #define CIRCLE_PATTERN_B 0 /* pattern of circles for additional particles */ @@ -119,9 +119,9 @@ #define OBSTACLE_OMEGA 300.0 /* obstacle rotation speed */ #define TWO_TYPES 0 /* set to 1 to have two types of particles */ -#define TYPE_PROPORTION 0.75 /* proportion of particles of first type */ +#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 0 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ +#define SYMMETRIZE_FORCE 1 /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */ #define CENTER_PX 0 /* set to 1 to center horizontal momentum */ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ #define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ @@ -139,10 +139,10 @@ #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.022 /* parameter controlling radius of particles */ -#define MU_B 0.022 /* parameter controlling radius of particles of second type */ -#define MU_ADD 0.022 /* parameter controlling radius of added particles */ -#define MU_ADD_B 0.022 /* parameter controlling radius of added particles */ +#define MU 0.02 /* parameter controlling radius of particles */ +#define MU_B 0.02 /* parameter controlling radius of particles of second type */ +#define MU_ADD 0.022 /* parameter controlling radius of added particles */ +#define MU_ADD_B 0.022 /* parameter controlling radius of added particles */ #define NPOLY 3 /* number of sides of polygon */ #define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define AWEDGE 0.5 /* opening angle of wedge, in units of Pi/2 */ @@ -151,8 +151,8 @@ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 40 /* number of grid point for grid of disks */ -#define NGRIDY 20 /* number of grid point for grid of disks */ +#define NGRIDX 15 /* number of grid point for grid of disks */ +#define NGRIDY 15 /* 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 */ @@ -160,6 +160,9 @@ #define NOBSX 24 #define NOBSY 14 /* obstacles for O_HEX obstacle pattern */ #define NTREES 15 /* number of trees in S_TREES */ +#define OFSSET_TREES 0.5 /* vertical offset in S_TREES_B */ +#define SLOPE_TREES 0.015 /* slope in S_TREES_B (default: 0.3) */ +#define SLOPE_TREES_B 0.015 /* slope in S_TREES_B (default: 0.25) */ #define X_SHOOTER -0.2 #define Y_SHOOTER -0.6 @@ -168,8 +171,8 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2300 /* number of frames of movie */ -#define NVID 100 /* number of iterations between images displayed on screen */ +#define NSTEPS 3800 /* number of frames of movie */ +#define NVID 100 /* number of iterations between images displayed on screen */ #define NSEG 25 /* number of segments of boundary of circles */ #define INITIAL_TIME 0 /* time after which to start saving frames */ #define OBSTACLE_INITIAL_TIME 0 /* time after which to start moving obstacle */ @@ -190,13 +193,13 @@ /* Plot type, see list in global_ljones.c */ -#define PLOT 5 +#define PLOT 13 #define PLOT_B 13 /* plot type for second movie */ /* Background color depending on particle properties */ #define COLOR_BACKGROUND 0 /* set to 1 to color background */ -#define BG_COLOR 10 /* type of background coloring, see list in global_ljones.c */ +#define BG_COLOR 1 /* type of background coloring, see list in global_ljones.c */ #define BG_COLOR_B 3 /* type of background coloring, see list in global_ljones.c */ #define OBSTACLE_COLOR 0 /* type of obstacle, see OC_ in global_ljones.c */ @@ -250,8 +253,8 @@ #define HUEAMP -50.0 /* amplitude of variation of hue for color scheme C_HUE */ #define COLOR_HUESHIFT 1.0 /* shift in color hue (for some cyclic palettes) */ -#define PRINT_PARAMETERS 0 /* set to 1 to print certain parameters */ -#define PRINT_TEMPERATURE 1 /* set to 1 to print current temperature */ +#define PRINT_PARAMETERS 1 /* set to 1 to print certain parameters */ +#define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */ #define PRINT_ANGLE 0 /* set to 1 to print obstacle orientation */ #define PRINT_OMEGA 0 /* set to 1 to print angular speed */ #define PRINT_PARTICLE_SPEEDS 0 /* set to 1 to print average speeds/momenta of particles */ @@ -259,6 +262,7 @@ #define PRINT_SEGMENTS_FORCE 0 /* set to 1 to print force on segments */ #define PRINT_NPARTICLES 0 /* print number of active particles */ #define PRINT_TYPE_PROP 0 /* print type proportion */ +#define PRINT_NABSORBED 1 /* print number of absorbed particles */ #define FORCE_FACTOR 0.1 /* factor controlling length of force vector */ /* particle properties */ @@ -268,7 +272,7 @@ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ #define PARTICLE_EMIN 0.0 /* energy of particle with coolest color */ -#define PARTICLE_EMAX 2000.0 /* energy of particle with hottest color */ +#define PARTICLE_EMAX 250000.0 /* energy of particle with hottest color */ #define PARTICLE_DMIN 200.0 /* energy of particle with largest local density */ #define PARTICLE_DMAX 500.0 /* energy of particle with largest local density */ #define SEGMENT_HUE_MIN 275.0 /* color of original segment */ @@ -277,15 +281,17 @@ #define OBSTACLE_VMAX 4.0 /* speed of obstacle with largest luminosity */ #define HUE_TYPE0 320.0 /* hue of particles of type 0 */ #define HUE_TYPE1 60.0 /* hue of particles of type 1 */ -#define HUE_TYPE2 320.0 /* hue of particles of type 2 */ +#define HUE_TYPE2 100.0 /* hue of particles of type 2 */ #define HUE_TYPE3 140.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 */ -#define HUE_TYPE6 130.0 /* hue of particles of type 6 */ -#define HUE_TYPE7 200.0 /* hue of particles of type 7 */ +#define HUE_TYPE4 180.0 /* hue of particles of type 4 */ +#define HUE_TYPE5 220.0 /* hue of particles of type 5 */ +#define HUE_TYPE6 260.0 /* hue of particles of type 6 */ +#define HUE_TYPE7 300.0 /* hue of particles of type 7 */ +#define HUE_TYPE8 330.0 /* hue of particles of type 7 */ #define BG_LOG_EKIN_SHIFT 1.0 /* constant in BG_LOG_EKIN background color scheme */ #define BG_FORCE_SLOPE 1.0e-6 /* constant in BG_FORCE backgound color scheme */ -#define BG_CHARGE_SLOPE 0.75 /* constant in BG_CHARGE backgound color scheme (default: 0.5) */ +#define BG_CHARGE_SLOPE 1.0 /* constant in BG_CHARGE backgound color scheme (default: 0.5) */ +#define CHARGE_HUE_RANGE 0.5 /* range of charge colors */ #define PARTICLE_LMAX 1.5e4 /* angular momentum particle with brightest color */ #define RANDOM_RADIUS 0 /* set to 1 for random particle radius */ @@ -295,36 +301,36 @@ #define ADAPT_DAMPING_TO_RADIUS 0.0 /* set to positive value to for friction prop to power of radius */ #define ADAPT_DAMPING_FACTOR 0.0 /* factor by which damping is adapted to radius */ #define DT_PARTICLE 1.0e-6 /* time step for particle displacement */ -#define KREPEL 50.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 3.5 /* Lennard-Jones equilibrium distance */ -#define EQUILIBRIUM_DIST_B 3.5 /* Lennard-Jones equilibrium distance for second type of particle */ +#define KREPEL 40.0 /* constant in repelling force between particles */ +#define EQUILIBRIUM_DIST 2.5 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST_B 2.5 /* Lennard-Jones equilibrium distance for second type of particle */ #define SEGMENT_FORCE_EQR 1.0 /* equilibrium distance factor for force from segments (default 1.5) */ #define REPEL_RADIUS 25.0 /* radius in which repelling force acts (in units of particle radius) */ #define DAMPING 0.0 /* damping coefficient of particles */ #define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */ #define DAMPING_ROT 0.0 /* damping coefficient for rotation of particles */ #define DAMPING_PAIRS 0.0 /* damping between paired particles */ -#define PARTICLE_MASS 1.0 /* mass of particle of radius MU */ -#define PARTICLE_MASS_B 1.0 /* mass of particle of radius MU_B */ -#define PARTICLE_ADD_MASS 1.0 /* mass of added particles */ +#define PARTICLE_MASS 2.0 /* mass of particle of radius MU */ +#define PARTICLE_MASS_B 2.0 /* mass of particle of radius MU_B */ +#define PARTICLE_ADD_MASS 2.0 /* mass of added particles */ #define PARTICLE_ADD_MASS_B 1.0 /* mass of added particles */ #define PARTICLE_INERTIA_MOMENT 0.1 /* moment of inertia of particle */ #define PARTICLE_INERTIA_MOMENT_B 0.1 /* moment of inertia of second type of particle */ -#define V_INITIAL 25.0 /* initial velocity range */ -#define V_INITIAL_ADD 0.0 /* initial velocity range for added particles */ +#define V_INITIAL 0.0 /* initial velocity range */ +#define V_INITIAL_ADD 5500.0 /* initial velocity range for added particles */ #define OMEGA_INITIAL 100.0 /* initial angular velocity range */ #define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */ #define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */ #define COULOMB_LJ_FACTOR 1.0 /* relative intensity of LJ interaction in I_COULOMB_LJ interaction (default: 0.01) */ -#define KCOULOMB_FACTOR 50.0 /* relative intensity of Coulomb interaction in I_COULOMB_LJ (default: 100.0) */ +#define KCOULOMB_FACTOR 500.0 /* relative intensity of Coulomb interaction in I_COULOMB_LJ (default: 100.0) */ #define OBSTACLE_DAMPING 0.0 /* damping of oscillating obstacles */ #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.0002 /* initial inverse temperature */ +#define BETA 0.00005 /* 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 5.0e11 /* harmonic potential of obstacles */ @@ -376,7 +382,7 @@ #define WIND_FORCE 1.35e6 /* force of wind */ #define WIND_YMIN -0.6 /* min altitude of region with wind */ -#define ROTATION 1 /* set to 1 to include rotation of particles */ +#define ROTATION 0 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ #define DIMENSION_FACTOR 0.25 /* scaling factor taking into account number of degrees of freedom */ #define KTORQUE 2.0e3 /* force constant in angular dynamics */ @@ -392,12 +398,12 @@ #define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */ #define INCREASE_BETA 0 /* set to 1 to increase BETA during simulation */ -#define BETA_SCHEDULE 5 /* type of temperature schedule, see TS_* in global_ljones */ -#define BETA_FACTOR 50.0 /* factor by which to change BETA during simulation */ +#define BETA_SCHEDULE 0 /* type of temperature schedule, see TS_* in global_ljones */ +#define BETA_FACTOR 50000.0 /* 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.5 /* number of temperature oscillations in BETA schedule */ +#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 */ -#define INITIAL_CONSTANT_PHASE 500 /* initial phase in which temperature is constant */ +#define INITIAL_CONSTANT_PHASE 200 /* initial phase in which temperature is constant */ #define MIDDLE_CONSTANT_PHASE 0 /* middle phase in which temperature is constant */ #define FINAL_DECREASE_PHASE 1 /* final phase in which temperature decreases */ #define FINAL_CONSTANT_PHASE 200 /* final phase in which temperature is constant */ @@ -437,24 +443,24 @@ #define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */ #define NPART_BOTTOM 100 /* number of particles at the bottom */ -#define ADD_PARTICLES 0 /* set to 1 to add particles */ -#define ADD_REGION 1 /* shape of add regions, cf ADD_* in global_ljones */ -#define ADD_TIME 0 /* time at which to add first particle */ -#define ADD_PERIOD 6 /* time interval between adding further particles */ +#define ADD_PARTICLES 1 /* set to 1 to add particles */ +#define ADD_REGION 0 /* shape of add regions, cf ADD_* in global_ljones */ +#define ADD_TIME 20 /* time at which to add first particle */ +#define ADD_PERIOD 10000 /* time interval between adding further particles */ #define ADD_TYPE 1 /* type of added particles */ #define N_ADD_PARTICLES 1 /* number of particles to add */ #define FINAL_NOADD_PERIOD 1800 /* final period where no particles are added */ #define SAFETY_FACTOR 10.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ -#define ADD_ALTERNATE_CHARGE 1 /* set to 1 to randomly select sign of added charge */ +#define ADD_ALTERNATE_CHARGE 0 /* set to 1 to randomly select sign of added charge */ #define TIME_DEPENDENT_ADD_CHARGE 0 /* set to 1 to have added charge depend on time */ #define ALTERNATE_CHARGE_PROPORTION 0.5 /* proportion of particles of opposite charge */ #define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */ -#define N_TRACER_PARTICLES 3000 /* number of tracer particles */ +#define N_TRACER_PARTICLES 1000 /* number of tracer particles */ #define TRACER_STEPS 5 /* number of tracer steps recorded between images */ #define TRAJECTORY_LENGTH 7000 /* length of recorded trajectory */ -#define TRAJECTORY_DRAW_LENGTH 250 /* length of drawn trajectory */ -#define TRACER_LUM_FACTOR 100.0 /* controls luminosity decrease of trajectories with time */ +#define TRAJECTORY_DRAW_LENGTH 1000 /* length of drawn trajectory */ +#define TRACER_LUM_FACTOR 50.0 /* controls luminosity decrease of trajectories with time */ #define TRACER_PARTICLE_MASS 1.0 /* relative mass of tracer particle */ #define TRAJECTORY_WIDTH 2 /* width of tracer particle trajectory */ @@ -506,28 +512,29 @@ #define POSITION_DEP_X -0.625 /* threshold value for position-dependent type */ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ -#define SPECIAL_IC 1 /* set to 1 for choosing specaial initial condition RD_INITIAL_COND */ -#define REACTION_DIFFUSION 1 /* set to 1 to simulate a chemical reaction (particles may change type) */ +#define SPECIAL_IC 0 /* set to 1 for choosing special initial condition RD_INITIAL_COND */ +#define REACTION_DIFFUSION 0 /* set to 1 to simulate a chemical reaction (particles may change type) */ #define REACTION_MAX_TIME 100000 /* time after which no reactions take place */ -#define RD_REACTION 11 /* type of reaction, see list in global_ljones.c */ -#define RD_TYPES 5 /* number of types in reaction-diffusion equation */ -#define RD_PLOT_TYPES 4 /* number of types shown in graph */ -#define RD_INITIAL_COND 5 /* initial condition of particles */ -#define REACTION_DIST 2.0 /* maximal distance for reaction to occur */ -#define REACTION_PROB 0.5 /* probability controlling reaction term */ +#define RD_REACTION 22 /* type of reaction, see list in global_ljones.c */ +#define RD_TYPES 8 /* number of types in reaction-diffusion equation */ +#define RD_PLOT_TYPES 8 /* number of types shown in graph */ +#define RD_INITIAL_COND 2 /* initial condition of particles */ +#define REACTION_DIST 2.8 /* 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 DELTAMAX 0.1 /* max orientation difference for pairing polygons */ #define CENTER_COLLIDED_PARTICLES 1 /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */ -#define EXOTHERMIC 1 /* set to 1 to make reaction exo/endothermic */ +#define EXOTHERMIC 0 /* set to 1 to make reaction exo/endothermic */ #define DELTA_EKIN -2.0e3 /* change of kinetic energy in reaction */ #define CORRECT_EQUILIBRIUM_POSITION 1 /* set to 1 to nudge particle dist towards eq dist */ -#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 4 /* maximal number of partners for CHEM_AGGREGATION reaction */ +#define NUDGE_FACTOR 0.0005 /* factor by which to correct particle distance */ +#define COLLISION_TIME 35 /* time during which collisions are shown */ +#define COLLISION_RADIUS 3.0 /* radius of discs showing collisions, in units of MU */ +#define DELTAVMAX 500.0 /* maximal deltav allowed for pairing molecules */ +#define AGREGMAX 3 /* maximal number of partners for CHEM_AGGREGATION reaction */ #define AGREG_DECOUPLE 12 /* minimal number of partners to decouple from thermostat */ -#define NEUTRALIZE_REACTING_PARTICLES 0 /* set to 1 for reacting particles to become neutral */ +#define NEUTRALIZE_REACTING_PARTICLES 1 /* set to 1 for reacting particles to become neutral */ #define CLUSTER_PARTICLES 0 /* set to 1 for particles to form rigid clusters */ #define CLUSTER_MAXSIZE 2 /* max size of clusters */ #define SMALL_CLUSTER_MAXSIZE 2 /* size limitation on smaller cluster */ @@ -540,7 +547,7 @@ #define MU_RATIO 0.666666667 /* ratio by which to increase radius */ #define PRINT_PARTICLE_NUMBER 0 /* set to 1 to print total number of particles */ -#define PLOT_PARTICLE_NUMBER 1 /* set to 1 to make of plot of particle number over time */ +#define PLOT_PARTICLE_NUMBER 0 /* set to 1 to make of plot of particle number over time */ #define PARTICLE_NB_PLOT_FACTOR 1.0 /* expected final number of particles over initial number */ #define PRINT_LEFT 0 /* set to 1 to print certain parameters at the top left instead of right */ #define PLOT_SPEEDS 0 /* set to 1 to add a plot of obstacle speeds (e.g. for rockets) */ @@ -573,27 +580,33 @@ #define KSPRING_PAIRS 5.0e9 /* 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 1 /* number of partners of particles - for DNA, set NPARTNERS_DNA */ +#define NPARTNERS 2 /* 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 99 /* type of pairing, see POLY_ in global_ljones.c */ +#define NARMS 4 /* number of "arms" for certain paring types */ +#define PAIRING_TYPE 99 /* 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 1.0 /* ratio between equilibrium distance and radius (default: 1.0) */ -#define MU_C 0.035 /* radius of partner particle */ +#define MU_C 0.022 /* radius of partner particle */ #define PARTICLE_MASS_C 1.0 /* mass or partner particle */ -#define CHARGE_C 1.0 /* charge of 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 ALTERNATE_POLY_CHARGE 0 /* 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 0 /* set to 1 to pair particle of type 1 */ -#define NPARTNERS_B 5 /* number of partners of particles */ +#define NPARTNERS_B 18 /* number of partners of particles */ #define NARMS_B 1 /* number of "arms" for certain paring types */ -#define PAIRING_TYPE_B 81 /* type of pairing, see POLY_ in global_ljones.c */ -#define MU_D 0.035 /* radius of partner particle */ +#define PAIRING_TYPE_B 5 /* type of pairing, see POLY_ in global_ljones.c */ +#define MU_D 0.022 /* radius of partner particle */ #define PARTICLE_MASS_D 1.0 /* mass or partner particle */ -#define CHARGE_D 1.0 /* charge of partner particle */ +#define CHARGE_D -1.0 /* charge of partner particle */ + +#define ADD_ABSORBERS 1 /* set to 1 to add absorbing discs */ +#define ABSORBER_PATTERN 4 /* pattern of absorbers, see AP_* in global_ljones */ +#define ABSORBER_X 0.0 +#define ABSORBER_Y 0.0 /* coordinates of first absorber */ +#define ABSORBER_R 0.16 /* radius of absorber */ #define NXMAZE 16 /* width of maze */ #define NYMAZE 16 /* height of maze */ @@ -621,25 +634,26 @@ /* constants related to evolution on a sphere */ #define SPHERE 1 /* set to 1 to compute evolution in spherical geometry */ -#define SIN_THETA_REG 0.05 /* regularization of sin(theta) for motion on sphere */ -#define POLAR_PADDING 0.1 /* region around poles that belong to the same hashcell */ -#define DRAW_SPHERE 1 /* set to 1 to draw 3D sphere */ +#define SIN_THETA_REG 0.01 /* regularization of sin(theta) for motion on sphere */ +#define POLAR_PADDING 0.05 /* region around poles that belong to the same hashcell */ +#define DRAW_SPHERE 1 /* set to 1 to draw 3D sphere */ #define NX_SPHERE 3000 #define NY_SPHERE 1500 /* number of points on sphere */ #define Z_SCALING_FACTOR 0.75 /* overall scaling factor of z axis for REP_PROJ_3D representation */ #define XY_SCALING_FACTOR 1.9 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define FLIPX -1.0 /* set to -1 to flip left/right */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ #define XSHIFT_3D -0.0 /* overall x shift for REP_PROJ_3D representation */ #define YSHIFT_3D -0.0 /* overall y shift for REP_PROJ_3D representation */ #define COS_VISIBLE -0.35 /* limit on cosine of normal to shown facets */ #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 750.0 /* total angle of rotation during simulation */ #define VIEWPOINT_TRAJ 1 /* type of viewpoint trajectory */ #define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */ -double light[3] = {0.40824829, -0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {3.0, -3.0, -2.5}; /* 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] = {-3.0, 0.0, 1.5}; /* location of observer for REP_PROJ_3D representation */ #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)) @@ -863,7 +877,8 @@ double temperature_schedule(int i) // else if ((i < INITIAL_TIME + t3)&&(t3 > t2)) beta = BETA*factor2; // else if ((i < INITIAL_TIME + t4)&&(t4 > t3)) // beta = BETA*exp(bexp2*(double)(i - INITIAL_TIME - t3) + 1.0e-10); - else beta = BETA*1.0e10; + else beta = BETA*factor2; +// else beta = BETA*1.0e10; // else beta = BETA; printf("beta = %.3lg\n", beta); return(beta); @@ -1082,7 +1097,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM double a, totalenergy = 0.0, damping, damping1, damping_rot1, direction, dmean, dratio; double ctheta, stheta, stheta_reg, ffx, ffy; static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; - int j, move, ncoup; + int j, move, ncoup, k, p; if (initial_phase) damping = INITIAL_DAMPING; else damping = DAMPING; @@ -1098,16 +1113,20 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM ctheta = cos(particle[j].yc); stheta = sin(particle[j].yc); stheta_reg = stheta + SIN_THETA_REG; - ffx = -2.0*ctheta*px[j]*py[j]; - ffy = stheta*ctheta*px[j]*px[j]; +// ffx = -2.0*ctheta*px[j]*py[j]; +// ffy = stheta*ctheta*px[j]*px[j]; + ffx = -2.0*ctheta*px[j]*py[j]*particle[j].mass_inv; + ffy = stheta*ctheta*px[j]*px[j]*particle[j].mass_inv; particle[j].vx = px[j] + 0.5*DT_PARTICLE*(particle[j].fx + ffx)/stheta_reg; particle[j].vy = py[j] + 0.5*DT_PARTICLE*(particle[j].fy + ffy); particle[j].omega = pangle[j] + 0.5*DT_PARTICLE*particle[j].torque; // particle[j].omega = pangle[j] + 0.5*DT_PARTICLE*particle[j].torque/stheta_reg; /* TODO: check/fix angular evolution */ - ffx = -2.0*ctheta*particle[j].vx*particle[j].vy; - ffy = stheta*ctheta*particle[j].vx*particle[j].vx; +// ffx = -2.0*ctheta*particle[j].vx*particle[j].vy; +// ffy = stheta*ctheta*particle[j].vx*particle[j].vx; + ffx = -2.0*ctheta*particle[j].vx*particle[j].vy*particle[j].mass_inv; + ffy = stheta*ctheta*particle[j].vx*particle[j].vx*particle[j].mass_inv; px[j] = particle[j].vx + 0.5*DT_PARTICLE*(particle[j].fx + ffx)/stheta_reg; py[j] = particle[j].vy + 0.5*DT_PARTICLE*(particle[j].fy + ffy); pangle[j] = particle[j].omega + 0.5*DT_PARTICLE*particle[j].torque; @@ -1252,6 +1271,8 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NM // } } +// wrap_particles(particle, molecule, cluster); + *ncoupled = ncoup; return(totalenergy); } @@ -1266,6 +1287,7 @@ double evolve_clusters(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXC 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; + double ctheta, stheta, stheta_reg, ffx, ffy; static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; int j, k, move, ncoup, mol, cl; short int moved[NMAXCIRCLES]; @@ -1277,16 +1299,41 @@ double evolve_clusters(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXC #pragma omp parallel for private(j,xi,totalenergy,a,move) for (j=0; j= 2)) @@ -2317,6 +2379,22 @@ void animation() } } /* end of for (n=0; n INITIAL_TIME)&&(REACTION_DIFFUSION)&&(i < REACTION_MAX_TIME)) { - ncollisions = update_types(particle, molecule, cluster, 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, tracer_n); if (EXOTHERMIC) params.beta *= 1.0/(1.0 + delta_energy/totalenergy); params.nactive = 0; for (j=0; j= dpoisson*dpoisson); + if (sphere) + { + far = far*(dist_sphere(x, y, point[k].xc, point[k].yc) >= dpoisson); + } + else + { + far = far*((x - point[k].xc)*(x - point[k].xc) + (y - point[k].yc)*(y - point[k].yc) >= dpoisson*dpoisson); + } /* new circle is in domain */ far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin); } @@ -612,9 +630,10 @@ int generate_poisson_discs(t_point *point, int nmax, double xmin, double xmax, d void init_particle_config(t_particle particles[NMAXCIRCLES]) /* initialise particle configuration */ { - int i, j, k, n, ncirc0, n_p_active, ncandidates = PDISC_CANDIDATES, naccepted; + int i, j, k, n, ncirc0, n_p_active, ncandidates = PDISC_CANDIDATES, naccepted, nx; double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = PDISC_DISTANCE*MU, xx[4], yy[4]; short int active_poisson[NMAXCIRCLES], far; + t_point *point; switch (CIRCLE_PATTERN) { case (C_SQUARE): @@ -635,12 +654,7 @@ void init_particle_config(t_particle particles[NMAXCIRCLES]) particles[n].active = 1; printf("(x,y) = (%.3lg,%.3lg)\n", particles[n].xc, particles[n].yc); - /* activate only circles that intersect the domain */ -// if ((particles[n].yc < INITYMAX + MU)&&(particles[n].yc > INITYMIN - MU)&&(particles[n].xc < INITXMAX + MU)&&(particles[n].xc > INITXMIN - MU)) particles[n].active = 1; -// else particles[n].active = 0; } - -// sleep(5); break; } case (C_HEX): @@ -848,6 +862,23 @@ void init_particle_config(t_particle particles[NMAXCIRCLES]) printf("Generated %i circles\n", ncircles); break; } + case (C_POISSON_DISC_SPHERE): + { + point = (t_point *)malloc(NMAXOBSTACLES*sizeof(t_point)); + ncircles = generate_poisson_discs(point, NMAXCIRCLES, INITXMIN, INITXMAX, INITYMIN, INITYMAX, PDISC_DISTANCE*MU, 1); + + for (n=0; n 0) { lum = (double)collisions[i].time/(double)COLLISION_TIME; - if (collisions[i].color == 0.0) for (j=0; j<3; j++) rgb[j] = lum; - else hsl_to_rgb_palette(collisions[i].color, 0.9, lum, rgb, COLOR_PALETTE); x = collisions[i].x; y = collisions[i].y; + if (COLOR_BACKGROUND) + { + lum1 = 1.0-lum; + cell = hash_cell(x1, y1); + rgb_bg[0] = hashgrid[cell].r; + rgb_bg[1] = hashgrid[cell].g; + rgb_bg[2] = hashgrid[cell].b; + if (collisions[i].color == 0.0) for (j=0; j<3; j++) rgb[j] = lum1*rgb_bg[j] + lum; + else + { + hsl_to_rgb_palette(collisions[i].color, 0.9, lum, rgb, COLOR_PALETTE); + for (j=0; j<3; j++) rgb[j] += lum1*rgb_bg[j]; + } + } + else + { + if (collisions[i].color == 0.0) for (j=0; j<3; j++) rgb[j] = lum; + else hsl_to_rgb_palette(collisions[i].color, 0.9, lum, rgb, COLOR_PALETTE); + } // printf("collision %i rgb = (%.3lf, %.3lf, %.3lf)\n", i, rgb[0], rgb[1], rgb[2]); @@ -8118,11 +8259,11 @@ void draw_collisions(t_collision *collisions, int ncollisions) } -void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, int *tracer_n, int plot) +void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_hashgrid hashgrid[HASHX*HASHY], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, int *tracer_n, int plot) /* draw tracer particle trajectory */ { - int i, j, time, p, width, imin, k; - double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum, lum1; + int i, j, time, p, width, imin, k, cell; + double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum, lum1, rgb_bg[3]; // blank(); glLineWidth(TRAJECTORY_WIDTH); @@ -8155,6 +8296,15 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], lum1 = 0.5*(1.0-lum); glColor3f(lum1 + lum*particle[tracer_n[j]].rgb[0], lum1 + lum*particle[tracer_n[j]].rgb[1], lum1 + lum*particle[tracer_n[j]].rgb[2]); } + else if (COLOR_BACKGROUND) + { + lum1 = 1.0-lum; + cell = hash_cell(x1, y1); + rgb_bg[0] = hashgrid[cell].r; + rgb_bg[1] = hashgrid[cell].g; + rgb_bg[2] = hashgrid[cell].b; + glColor3f(lum1*rgb_bg[0] + lum*particle[tracer_n[j]].rgb[0], lum1*rgb_bg[1] + lum*particle[tracer_n[j]].rgb[1], lum1*rgb_bg[2] + lum*particle[tracer_n[j]].rgb[2]); + } else glColor3f(lum*particle[tracer_n[j]].rgb[0], lum*particle[tracer_n[j]].rgb[1], lum*particle[tracer_n[j]].rgb[2]); if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) @@ -8174,7 +8324,21 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], time = traj_position + traj_length - i; lum = 0.9 - TRACER_LUM_FACTOR*(double)time/(double)(TRAJECTORY_LENGTH*TRACER_STEPS); if (lum < 0.0) lum = 0.0; - glColor3f(lum*particle[tracer_n[j]].rgb[0], lum*particle[tracer_n[j]].rgb[1], lum*particle[tracer_n[j]].rgb[2]); + if (FILL_OBSTACLE_TRIANGLES) + { + lum1 = 0.5*(1.0-lum); + glColor3f(lum1 + lum*particle[tracer_n[j]].rgb[0], lum1 + lum*particle[tracer_n[j]].rgb[1], lum1 + lum*particle[tracer_n[j]].rgb[2]); + } + else if (COLOR_BACKGROUND) + { + lum1 = 1.0-lum; + cell = hash_cell(x1, y1); + rgb_bg[0] = hashgrid[cell].r; + rgb_bg[1] = hashgrid[cell].g; + rgb_bg[2] = hashgrid[cell].b; + glColor3f(lum1*rgb_bg[0] + lum*particle[tracer_n[j]].rgb[0], lum1*rgb_bg[1] + lum*particle[tracer_n[j]].rgb[1], lum1*rgb_bg[2] + lum*particle[tracer_n[j]].rgb[2]); + } + else glColor3f(lum*particle[tracer_n[j]].rgb[0], lum*particle[tracer_n[j]].rgb[1], lum*particle[tracer_n[j]].rgb[2]); if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); @@ -8191,7 +8355,21 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], time = traj_position - i; lum = 0.9 - TRACER_LUM_FACTOR*(double)time/(double)(TRAJECTORY_LENGTH*TRACER_STEPS); if (lum < 0.0) lum = 0.0; - glColor3f(lum*particle[tracer_n[j]].rgb[0], lum*particle[tracer_n[j]].rgb[1], lum*particle[tracer_n[j]].rgb[2]); + if (FILL_OBSTACLE_TRIANGLES) + { + lum1 = 0.5*(1.0-lum); + glColor3f(lum1 + lum*particle[tracer_n[j]].rgb[0], lum1 + lum*particle[tracer_n[j]].rgb[1], lum1 + lum*particle[tracer_n[j]].rgb[2]); + } + else if (COLOR_BACKGROUND) + { + lum1 = 1.0-lum; + cell = hash_cell(x1, y1); + rgb_bg[0] = hashgrid[cell].r; + rgb_bg[1] = hashgrid[cell].g; + rgb_bg[2] = hashgrid[cell].b; + glColor3f(lum1*rgb_bg[0] + lum*particle[tracer_n[j]].rgb[0], lum1*rgb_bg[1] + lum*particle[tracer_n[j]].rgb[1], lum1*rgb_bg[2] + lum*particle[tracer_n[j]].rgb[2]); + } + else glColor3f(lum*particle[tracer_n[j]].rgb[0], lum*particle[tracer_n[j]].rgb[1], lum*particle[tracer_n[j]].rgb[2]); if ((lum > 0.1)&&(module2(x2 - x1, y2 - y1) < 0.25*(YMAX - YMIN))) draw_line(x1, y1, x2, y2); @@ -8601,7 +8779,7 @@ void color_background(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAX counter = 1 - counter; } -void draw_background(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY]) +void draw_background(t_particle particle[NMAXCIRCLES], t_segment segment[NMAXSEGMENTS], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY]) /* color background according to particle properties */ { int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs; @@ -8974,6 +9152,33 @@ void draw_particles(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRC } } + +void draw_absorbers(t_absorber absorber[NMAX_ABSORBERS]) +/* draw the absorbing discs */ +{ + int i, p; + double rgb[3], s; + + rgb[0] = 0.5; + rgb[1] = 0.5; + rgb[2] = 0.5; + + for (i=0; i 0.1) + draw_colored_ellipse_precomp(absorber[i].xc + p*DPI, absorber[i].yc, absorber[i].radius/s, absorber[i].radius, rgb); + else + { + if (absorber[i].yc > 0.0) + draw_colored_rectangle(XMIN, YMAX - absorber[i].radius, XMAX, YMAX, rgb); + else + draw_colored_rectangle(XMIN, YMIN, XMAX, YMIN + absorber[i].radius, rgb); + } + } +} + void obstacle_hue(t_obstacle obstacle, double rgb[3]) /* compute obstacle color hue */ { @@ -9775,6 +9980,13 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_ if (PRINT_LEFT) write_text(xtext, y, message); else write_text(xmidtext, y, message); } + if (PRINT_NABSORBED) + { + erase_area_hsl(xbox, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "%i in pocket", params.nabsorbed); + write_text(xtext + 0.1, y, message); + } if (DECREASE_CONTAINER_SIZE) /* print density */ { density = (double)ncircles/((lengthcontainer)*(INITYMAX - INITYMIN)); @@ -11140,15 +11352,16 @@ void compute_partner_force(int j, int n, double eq_distance, double eq_angle, do double x1, x2, y1, y2, rmax, alpha, cosa, rj, rn, phi, dvx, dvy, pscal, newdistance, xg, yg, vxg, vyg; int p, q; - dx = particle[n].xc - particle[j].xc; - dy = particle[n].yc - particle[j].yc; - - if (dx > 0.5*(BCXMAX - BCXMIN)) dx -= (BCXMAX-BCXMIN); - else if (dx < -0.5*(BCXMAX - BCXMIN)) dx += (BCXMAX-BCXMIN); - if (dy > 0.5*(BCYMAX - BCYMIN)) dy -= (BCYMAX-BCYMIN); - else if (dy < -0.5*(BCYMAX - BCYMIN)) dy += (BCYMAX-BCYMIN); - + if (BOUNDARY_COND != BC_SPHERE) + { + dx = particle[n].xc - particle[j].xc; + dy = particle[n].yc - particle[j].yc; + if (dx > 0.5*(BCXMAX - BCXMIN)) dx -= (BCXMAX-BCXMIN); + else if (dx < -0.5*(BCXMAX - BCXMIN)) dx += (BCXMAX-BCXMIN); + if (dy > 0.5*(BCYMAX - BCYMIN)) dy -= (BCYMAX-BCYMIN); + else if (dy < -0.5*(BCYMAX - BCYMIN)) dy += (BCYMAX-BCYMIN); + } if ((PAIR_PARTICLES)&&(PAIRING_TYPE == POLY_SEG_POLYGON)) { @@ -11174,12 +11387,21 @@ void compute_partner_force(int j, int n, double eq_distance, double eq_angle, do } // 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 (BOUNDARY_COND == BC_SPHERE) + { + r = distance_sphere(j, n, particle, &ca, &sa); + if (r < 1.0e-10) r = 1.0e-10; + if (r > 2.0*eq_distance) r = 2.0*eq_distance; + } + 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; + ca = dx/r; + sa = dy/r; + } /* TODO: adjust max distance */ // if (r < 2.0*eq_distance) force = KSPRING_PAIRS*(r - eq_distance); @@ -11208,13 +11430,26 @@ void compute_partner_force(int j, int n, double eq_distance, double eq_angle, do /* TEST - experimental */ if (CORRECT_EQUILIBRIUM_POSITION) { - newdistance = 0.995*r + 0.005*eq_distance; - xg = 0.5*(particle[n].xc + particle[j].xc); - yg = 0.5*(particle[n].yc + particle[j].yc); - particle[n].xc = xg + 0.5*ca*newdistance; - particle[n].yc = yg + 0.5*sa*newdistance; - particle[j].xc = xg - 0.5*ca*newdistance; - particle[j].yc = yg - 0.5*sa*newdistance; + newdistance = (1.0 - NUDGE_FACTOR)*r + NUDGE_FACTOR*eq_distance; + if (BOUNDARY_COND == BC_SPHERE) + { +// printf("phi1 = %.3lg, theta1 = %.3lg\n", particle[j].xc, particle[j].yc); +// printf("phi2 = %.3lg, theta2 = %.3lg\n", particle[n].xc, particle[n].yc); + + compute_convex_combination_sphere(particle[j].xc, particle[j].yc, particle[n].xc, particle[n].yc, &particle[j].xc, &particle[j].yc, &particle[n].xc, &particle[n].yc, newdistance/r); +// particle[j].xc = xg; +// particle[j].yc = yg; +// printf("xg = %.3lg, yg = %.3lg\n", xg, yg); + } + else + { + xg = 0.5*(particle[n].xc + particle[j].xc); + yg = 0.5*(particle[n].yc + particle[j].yc); + particle[n].xc = xg + 0.5*ca*newdistance; + particle[n].yc = yg + 0.5*sa*newdistance; + particle[j].xc = xg - 0.5*ca*newdistance; + particle[j].yc = yg - 0.5*sa*newdistance; + } // vxg = 0.5*(particle[j].vx + particle[n].vx); // vyg = 0.5*(particle[j].vy + particle[n].vy); // particle[j].vx = vxg; @@ -11715,6 +11950,11 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) /* add copies in case of particle pairing */ if (PAIR_PARTICLES) init_particle_pairs(particle, molecule); + /* TEST */ +// if (PAIRING_TYPE == POLY_PLUSMINUS) +// for (i=0; i NGRIDX/2)||(moli-molk > NGRIDX/2))) triangle_condition = 0; } + /* also limit number of partners of second particle */ + if (no_triangles == -2) + { + if (nq > maxpartners + 1) triangle_condition = 0; + if (moli == molk) triangle_condition = 0; + for (q=0; q 2)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) - ncollisions = chem_split(i, 1, particle[i].type - 1, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 1, particle[i].type - 1, particle, collisions, ncollisions, inv_masses, radii, tracer_n); } printf("%i collisions\n", ncollisions); @@ -15695,7 +15967,7 @@ int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCL { k = rand()%(type-1) + 1; // printf("Splitting type %i into type %i and type %i\n", type, k, type - k); - ncollisions = chem_split(i, k, type - k, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, k, type - k, particle, collisions, ncollisions, inv_masses, radii, tracer_n); } } @@ -15725,7 +15997,7 @@ int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCL } else if ((particle[i].active)&&(particle[i].type == 3)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) { - ncollisions = chem_split(i, 1, 2, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 1, 2, particle, collisions, ncollisions, inv_masses, radii, tracer_n); if (EXOTHERMIC) *delta_e -= (double)(ncollisions - oldncollisions)*DELTA_EKIN; } } @@ -15749,12 +16021,12 @@ int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCL } else if ((particle[i].active)&&(particle[i].type == 3)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) { - ncollisions = chem_split(i, 1, 2, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 1, 2, particle, collisions, ncollisions, inv_masses, radii, tracer_n); if (EXOTHERMIC) *delta_e -= (double)(ncollisions - oldncollisions)*DELTA_EKIN; } else if ((particle[i].active)&&(particle[i].type == 4)&&((double)rand()/RAND_MAX < DISSOCIATION_PROB)) { - ncollisions = chem_split(i, 1, 3, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 1, 3, particle, collisions, ncollisions, inv_masses, radii, tracer_n); if (EXOTHERMIC) *delta_e -= (double)(ncollisions - oldncollisions)*DELTA_EKIN; } } @@ -15799,12 +16071,12 @@ int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCL } case (7): { - ncollisions = chem_split(i, 3, 3, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 3, 3, particle, collisions, ncollisions, inv_masses, radii, tracer_n); break; } case (8): { - ncollisions = chem_split(i, 5, 5, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 5, 5, particle, collisions, ncollisions, inv_masses, radii, tracer_n); break; } } @@ -15919,7 +16191,29 @@ int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCL { if (particle[i].npartners < AGREGMAX) for (k=0; k<3; k++) { - ncollisions = chem_multi_glue_molecule(i, k, AGREGMAX, 0, 0, 0, 0, particle, molecule, cluster, collisions, ncollisions, REACTION_PROB, 0); + ncollisions = chem_multi_glue_molecule(i, k, AGREGMAX, 2, 0, 0, 0, particle, molecule, cluster, collisions, ncollisions, REACTION_PROB, 0); + } + + /* decouple particles with several partners from thermostat */ + if (particle[i].npartners >= 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_AGGREGATION_MAX): + { + for (i=0; i Sr + Xe */ { if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) - ncollisions = chem_split(i, 3, 4, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 3, 4, particle, collisions, ncollisions, inv_masses, radii, tracer_n); break; } case (6): /* B -> n + n */ { if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) - ncollisions = chem_split(i, 1, 1, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 1, 1, particle, collisions, ncollisions, inv_masses, radii, tracer_n); break; } } @@ -16475,13 +16769,13 @@ int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCL case (7): /* A -> Sr + Xe */ { if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) - ncollisions = chem_split(i, 5, 6, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 5, 6, particle, collisions, ncollisions, inv_masses, radii, tracer_n); break; } case (8): /* B -> n + n */ { if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) - ncollisions = chem_split(i, 4, 4, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 4, 4, particle, collisions, ncollisions, inv_masses, radii, tracer_n); break; } } @@ -16516,13 +16810,13 @@ int update_types(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCL case (8): /* A -> Sr + Xe */ { if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) - ncollisions = chem_split(i, 5, 6, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 5, 6, particle, collisions, ncollisions, inv_masses, radii, tracer_n); break; } case (9): /* B -> n + n */ { if ((double)rand()/RAND_MAX < DISSOCIATION_PROB) - ncollisions = chem_split(i, 4, 4, particle, collisions, ncollisions, inv_masses, radii); + ncollisions = chem_split(i, 4, 4, particle, collisions, ncollisions, inv_masses, radii, tracer_n); break; } } @@ -17215,18 +17509,18 @@ int wrap_particles_periodic_bc(t_particle particle[NMAXCIRCLES], t_molecule mole return(0); } - fprintf(lj_log, "%i particles and containing molecules to be moved\n", move); + fprintf(lj_log, "[wrap_particles_periodic_bc] %i particles and containing molecules to be moved\n", move); /* determine molecules to be moved */ for (mol = 0; mol < ncircles; mol++) move_molecule[mol] = 0; for (i=0; i BCXMAX + paddingx) + { + move_part[i] = 1; + deltax[i] = -dx; + } + else if (particle[i].xc < BCXMIN - paddingx) + { + move_part[i] = 1; + deltax[i] = dx; + } +// if (particle[i].yc > BCYMAX + paddingy) +// { +// move_part[i] = 1; +// deltay[i] = -dy; +// } +// else if (particle[i].yc < BCYMIN - paddingy) +// { +// move_part[i] = 1; +// deltay[i] = dy; +// } + move += move_part[i]; +// if (move_part[i]) +// fprintf(lj_log, "Particle %i at (%.5lg, %.5lg) should be moved by (%.5lg, %.5lg)\n", i, particle[i].xc, particle[i].yc, deltax[i], deltay[i]); + } + + if (move == 0) + { + free(move_part); + free(move_molecule); + free(molecule_refparticle); + free(move_cluster); + free(cluster_refparticle); + free(deltax); + free(deltay); + return(0); + } + + fprintf(lj_log, "[wrap_particles_periodic_bc] %i particles and containing molecules to be moved\n", move); + /* determine molecules to be moved */ + for (mol = 0; mol < ncircles; mol++) move_molecule[mol] = 0; + + for (i=0; i DPI) absorber[k].xc -= DPI; + } + break; + } + case (AP_ICOSA): + { + alpha = acos(1.0/sqrt(5.0)); + + absorber[nabsorb].xc = 0.0; + absorber[nabsorb].yc = 0.0; + nabsorb++; + + for (k=0; k<5; k++) + { + absorber[nabsorb].xc = (double)k*DPI/5.0; + absorber[nabsorb].yc = alpha; + nabsorb++; + } + + for (k=0; k<5; k++) + { + absorber[nabsorb].xc = (double)k*DPI/5.0 + PI/5.0; + absorber[nabsorb].yc = PI - alpha; + nabsorb++; + } + + absorber[nabsorb].xc = 0.0; + absorber[nabsorb].yc = PI; + nabsorb++; + + for (k=0; k DPI) absorber[k].xc -= DPI; + } + break; + } + } + return(nabsorb); +} + int wrap_particles(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES]) /* wrap particles, molecules, clusters, for certain boundary conditions */ { @@ -17286,6 +17850,10 @@ int wrap_particles(t_particle particle[NMAXCIRCLES], t_molecule molecule[NMAXCIR { return(wrap_particles_periodic_bc(particle, molecule, cluster)); } + case (BC_SPHERE): + { + return(wrap_particles_sphere_bc(particle, molecule, cluster)); + } default: /* TODO */ return(0); } @@ -17303,7 +17871,7 @@ void draw_frame_2d(int i, int plot, int bg_color, int ncollisions, int traj_posi t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], t_group_data *group_speeds, t_group_segments *segment_group, t_belt *conveyor_belt, int *tracer_n, - t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES]) + t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES], t_absorber absorber[NMAX_ABSORBERS]) /* draw a movie frame */ { printf("Drawing frame\n"); @@ -17311,19 +17879,16 @@ void draw_frame_2d(int i, int plot, int bg_color, int ncollisions, int traj_posi compute_all_particle_colors(particle, cluster, plot); if ((COLOR_BACKGROUND)&&(bg_color > 0)) { - /* old version */ -// color_background(particle, obstacle, bg_color, hashgrid); - - /* new version */ - compute_background_color(particle, obstacle, bg_color, hashgrid); - draw_background(particle, obstacle, bg_color, hashgrid); + compute_background_color(particle, segment, obstacle, bg_color, hashgrid); + draw_background(particle, segment, obstacle, bg_color, hashgrid); } // else if (!TRACER_PARTICLE) blank(); - if ((REACTION_DIFFUSION)&&(ncollisions > 0)) draw_collisions(collisions, ncollisions); + if ((REACTION_DIFFUSION)&&(ncollisions > 0)) draw_collisions(collisions, ncollisions, hashgrid); if (DRAW_OBSTACLE_LINKS) draw_obstacle_triangles(obstacle, otriangle, ofacet); // if (DRAW_OBSTACLE_LINKS) old_draw_obstacle_triangles(obstacle); - if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length, particle, cluster, tracer_n, plot); + if (TRACER_PARTICLE) draw_trajectory(trajectory, hashgrid, traj_position, traj_length, particle, cluster, tracer_n, plot); + if (ADD_ABSORBERS) draw_absorbers(absorber); draw_particles(particle, cluster, plot, params.beta, collisions, ncollisions, bg_color, hashgrid, params); draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, conveyor_belt, wall); if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, refresh); @@ -17355,17 +17920,26 @@ void draw_frame_3d(int i, int plot, int bg_color, int ncollisions, int traj_posi t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], t_group_data *group_speeds, t_group_segments *segment_group, t_belt *conveyor_belt, int *tracer_n, - t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE]) + t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], + t_absorber absorber[NMAX_ABSORBERS]) /* draw a movie frame */ { +// static int first; +// +// if (first) +// { +// if (ADD_ABSORBERS) draw_absorbers_sphere(absorber, wsphere); +// first = 0; +// } + /* initialize background colors */ - if (COLOR_BACKGROUND) - compute_background_color(particle, obstacle, bg_color, hashgrid); + if ((COLOR_BACKGROUND)&&(bg_color > 0)) + compute_background_color(particle, segment, obstacle, bg_color, hashgrid); /* compute all particle colors */ compute_all_particle_colors(particle, cluster, plot); /* initialize sphere radius and color */ - init_sphere_radius(particle, hashgrid, cluster, trajectory, traj_position, traj_length, wsphere, tracer_n, plot); + init_sphere_radius(particle, hashgrid, cluster, trajectory, segment, traj_position, traj_length, wsphere, tracer_n, plot, bg_color, absorber); /* initialize light angle (cf compute_light_angle_sphere_rde() in sub_rde.c) */ compute_light_angle_sphere(wsphere); /* draw sphere with particles (cf draw_wave_sphere_3d_rde() in sub_rde.c) */ @@ -17389,11 +17963,12 @@ void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_positio t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], t_group_data *group_speeds, t_group_segments *segment_group, t_belt *conveyor_belt, int *tracer_n, - t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE]) + t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], + t_absorber absorber[NMAX_ABSORBERS]) /* draw a movie frame */ { - if (DRAW_SPHERE) draw_frame_3d(i, plot, bg_color, ncollisions, traj_position, traj_length, - wall, pressure, pleft, pright, currents, particle_numbers, refresh, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere); + if (DRAW_SPHERE) draw_frame_3d(i, plot, bg_color, ncollisions, traj_position, + traj_length, wall, pressure, pleft, pright, currents, particle_numbers, refresh, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere, absorber); else draw_frame_2d(i, plot, bg_color, ncollisions, traj_position, traj_length, - wall, pressure, pleft, pright, currents, particle_numbers, refresh, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet); + wall, pressure, pleft, pright, currents, particle_numbers, refresh, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, absorber); } diff --git a/sub_lj_sphere.c b/sub_lj_sphere.c index aa9ce6f..7e529b5 100644 --- a/sub_lj_sphere.c +++ b/sub_lj_sphere.c @@ -92,6 +92,64 @@ double dist_sphere(double phi1, double theta1, double phi2, double theta2) } +void compute_midpoint_sphere(double phi1, double theta1, double phi2, double theta2, double *phi, double *theta) +/* compute midpoint between two points on the sphere */ +{ + double x1, y1, z1, x2, y2, z2, x3, y3, z3, r; + + x1 = cos(phi1)*sin(theta1); + y1 = sin(phi1)*sin(theta1); + z1 = -cos(theta1); + + x2 = cos(phi2)*sin(theta2); + y2 = sin(phi2)*sin(theta2); + z2 = -cos(theta2); + + x3 = x1 + x2; + y3 = y1 + y2; + z3 = z1 + z2; + + r = sqrt(x3*x3 + y3*y3 + z3*z3); + *theta = acos(-z3/r); + *phi = argument(x3, y3); + if (*phi < 0.0) *phi += DPI; +} + + +void compute_convex_combination_sphere(double phi1, double theta1, double phi2, double theta2, double *phia, double *thetaa, double *phib, double *thetab, double lambda) +/* compute convex combination of two points on the sphere */ +{ + double x1, y1, z1, x2, y2, z2, x3, y3, z3, r, lam1; + + lam1 = 1.0 - lambda; + + x1 = cos(phi1)*sin(theta1); + y1 = sin(phi1)*sin(theta1); + z1 = -cos(theta1); + + x2 = cos(phi2)*sin(theta2); + y2 = sin(phi2)*sin(theta2); + z2 = -cos(theta2); + + x3 = lambda*x1 + lam1*x2; + y3 = lambda*y1 + lam1*y2; + z3 = lambda*z1 + lam1*z2; + + r = sqrt(x3*x3 + y3*y3 + z3*z3); + *thetaa = acos(-z3/r); + *phia = argument(x3, y3); + if (*phia < 0.0) *phia += DPI; + + x3 = lam1*x1 + lambda*x2; + y3 = lam1*y1 + lambda*y2; + z3 = lam1*z1 + lambda*z2; + + r = sqrt(x3*x3 + y3*y3 + z3*z3); + *thetab = acos(-z3/r); + *phib = argument(x3, y3); + if (*phib < 0.0) *phib += DPI; +} + int in_polygon(double x, double y, double r, int npoly, double apoly) /* test whether (x,y) is in regular polygon of npoly sides inscribed in circle of radius r, turned by apoly Pi/2 */ { @@ -176,6 +234,7 @@ double type_hue(int type) if (RD_REACTION == CHEM_BZ) return(HUE_TYPE2); else return(HUE_TYPE7); } + case (8): return(HUE_TYPE8); default: { if (RD_REACTION == CHEM_BZ) @@ -350,7 +409,7 @@ void compute_particle_colors(t_particle particle, t_cluster cluster[NMAXCIRCLES] } case (P_CHARGE): { - hue = (-0.6*tanh(particle.charge)+1.0)*180.0; + hue = (-CHARGE_HUE_RANGE*tanh(BG_CHARGE_SLOPE*particle.charge)+1.0)*180.0; break; } case (P_MOL_ANGLE): @@ -631,7 +690,7 @@ void compute_all_particle_colors(t_particle particle[NMAXCIRCLES], t_cluster clu } -void compute_background_color(t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY]) +void compute_background_color(t_particle particle[NMAXCIRCLES], t_segment segment[NMAXSEGMENTS], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY]) /* color background according to particle properties */ { int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs; @@ -1139,12 +1198,19 @@ double dist_point_to_particle_sphere(int i, double phi, double psi, t_particle* void init_3d() /* initialisation of window */ { + double width, height; + + width = 2.0*(WINWIDTH/1760.0); + height = WINHEIGHT/990.0; + glLineWidth(3); glClearColor(0.0, 0.0, 0.0, 1.0); glClear(GL_COLOR_BUFFER_BIT); - glOrtho(-2.0, 2.0, -1.0, 1.0 , -1.0, 1.0); +// glOrtho(-2.0, 2.0, -1.0, 1.0 , -1.0, 1.0); + + glOrtho(-width, width, -height, height, -1.0, 1.0); } void xyz_to_xy(double x, double y, double z, double xy_out[2]) @@ -1177,7 +1243,7 @@ void xyz_to_xy(double x, double y, double z, double xy_out[2]) xinter[1] = t*observer[1] + (1.0-t)*y; xinter[2] = t*observer[2] + (1.0-t)*z; - xy_out[0] = XSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*h[0] + xinter[1]*h[1]); + xy_out[0] = XSHIFT_3D + FLIPX*XY_SCALING_FACTOR*(xinter[0]*h[0] + xinter[1]*h[1]); xy_out[1] = YSHIFT_3D + XY_SCALING_FACTOR*(xinter[0]*v[0] + xinter[1]*v[1] + xinter[2]*v[2]); } @@ -1222,6 +1288,8 @@ void init_lj_sphere(t_lj_sphere *wsphere) wsphere[i*NY_SPHERE+j].radius = 1.0; wsphere[i*NY_SPHERE+j].cos_angle_sphere = wsphere[i*NY_SPHERE+j].x*light[0] + wsphere[i*NY_SPHERE+j].y*light[1] + wsphere[i*NY_SPHERE+j].z*light[2]; + + wsphere[i*NY_SPHERE+j].locked = 0; } } } @@ -1296,8 +1364,8 @@ double test_distance(double phi, double psi, int i, t_particle* particle) void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRCLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE]) /* compute the effect at point (i,j) of adding a particle */ { - int draw_normal_particle = 1; - double x, y, r, dist, ca, sa, x1, y1, x2, y2, d1, d2, r2; + int draw_normal_particle = 1, k; + double x, y, r, dist, ca, sa, x1, y1, x2, y2, d1, d2, r2, omega, h; static double dphi, dtheta; static int first = 1; @@ -1335,7 +1403,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC r2 = 0.49*r*r; - if (d1 < r2) + if ((d1 < r2)&&(!wsphere[i*NY_SPHERE+j].locked)) { wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0]; wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1]; @@ -1365,7 +1433,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC r2 = 0.49*r*r; - if (d1 < r2) + if ((d1 < r2)&&(!wsphere[i*NY_SPHERE+j].locked)) { wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0]; wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1]; @@ -1395,7 +1463,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC r2 = 0.9*r*r; - if (d1 < r2) + if ((d1 < r2)&&(!wsphere[i*NY_SPHERE+j].locked)) { wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0]; wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1]; @@ -1421,7 +1489,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC r2 = 0.7*r*r; - if (d1 < r2) + if ((d1 < r2)&&(!wsphere[i*NY_SPHERE+j].locked)) { wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0]; wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1]; @@ -1433,6 +1501,50 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC } break; } + case (CHEM_POLYMER): + { + dist = dist_point_to_particle_sphere(part, x, y, particle, &ca, &sa); + x2 = dist*ca; + y2 = dist*sa; + r = 1.2*MU; + + if ((dist < r)&&(!wsphere[i*NY_SPHERE+j].locked)) + { + wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0]; + wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1]; + wsphere[i*NY_SPHERE+j].b = particle[part].rgb[2]; + wsphere[i*NY_SPHERE+j].radius += 1.5*sqrt(r*r - dist*dist); + } + + if (particle[part].type > 2) + { + omega = DPI/(double)(particle[part].type - 2); + + for (k=0; k wsphere[i*NY_SPHERE+j].radius)&&(!wsphere[i*NY_SPHERE+j].locked)) + { + wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0]; + wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1]; + wsphere[i*NY_SPHERE+j].b = particle[part].rgb[2]; + wsphere[i*NY_SPHERE+j].radius = h; + } + } + } + } + + draw_normal_particle = 0; + break; + } } } @@ -1442,19 +1554,23 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC if (dist < r) { - wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0]; - wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1]; - wsphere[i*NY_SPHERE+j].b = particle[part].rgb[2]; - wsphere[i*NY_SPHERE+j].radius += 1.5*sqrt(r*r - dist*dist); + h = 1.0 + 1.5*sqrt(r*r - dist*dist); + if ((h > wsphere[i*NY_SPHERE+j].radius)&&(!wsphere[i*NY_SPHERE+j].locked)) + { + wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0]; + wsphere[i*NY_SPHERE+j].g = particle[part].rgb[1]; + wsphere[i*NY_SPHERE+j].b = particle[part].rgb[2]; + wsphere[i*NY_SPHERE+j].radius = h; + } } } } -void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot) +void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_hashgrid hashgrid[HASHX*HASHY], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot) /* draw tracer particle trajectory */ { int i, j, i0, j0, time, p, q, width, imin, cell, i1, j1; - double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum, lum1; + double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum, lum1, rgb_bg[3]; static double dphi, dtheta; static int first = 1; @@ -1489,6 +1605,17 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI if (lum < 0.0) lum = 0.0; lum1 = 1.0 - lum; + if (COLOR_BACKGROUND) + { + cell = hash_cell(x1, y1); + if (!wsphere[cell].locked) + { + rgb_bg[0] = hashgrid[cell].r; + rgb_bg[1] = hashgrid[cell].g; + rgb_bg[2] = hashgrid[cell].b; + } + } + if ((x2 != x1)||(y2 != y1)) for (p=-1; p<2; p++) for (q=-1; q<2; q++) { @@ -1499,9 +1626,21 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI if (j1 < 0) j1 = 0; if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1; cell = i1*NY_SPHERE+j1; - wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1; - wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1; - wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1; + if (!wsphere[cell].locked) + { + if (COLOR_BACKGROUND) + { + wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1*rgb_bg[0]; + wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1*rgb_bg[1]; + wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1*rgb_bg[2]; + } + else + { + wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1; + wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1; + wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1; + } + } } } else @@ -1522,6 +1661,17 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI if (lum < 0.0) lum = 0.0; lum1 = 1.0 - lum; + if (COLOR_BACKGROUND) + { + cell = hash_cell(x1, y1); + if (!wsphere[cell].locked) + { + rgb_bg[0] = hashgrid[cell].r; + rgb_bg[1] = hashgrid[cell].g; + rgb_bg[2] = hashgrid[cell].b; + } + } + if ((x2 != x1)||(y2 != y1)) for (p=-1; p<2; p++) for (q=-1; q<2; q++) { @@ -1532,9 +1682,21 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI if (j1 < 0) j1 = 0; if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1; cell = i1*NY_SPHERE+j1; - wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1; - wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1; - wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1; + if (!wsphere[cell].locked) + { + if (COLOR_BACKGROUND) + { + wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1*rgb_bg[0]; + wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1*rgb_bg[1]; + wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1*rgb_bg[2]; + } + else + { + wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1; + wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1; + wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1; + } + } } } for (i=imin; i < traj_position-1; i++) @@ -1550,6 +1712,17 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI if (lum < 0.0) lum = 0.0; lum1 = 1.0 - lum; + if (COLOR_BACKGROUND) + { + cell = hash_cell(x1, y1); + if (!wsphere[cell].locked) + { + rgb_bg[0] = hashgrid[cell].r; + rgb_bg[1] = hashgrid[cell].g; + rgb_bg[2] = hashgrid[cell].b; + } + } + if ((x2 != x1)||(y2 != y1)) for (p=-1; p<2; p++) for (q=-1; q<2; q++) { @@ -1560,16 +1733,119 @@ void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTI if (j1 < 0) j1 = 0; if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1; cell = i1*NY_SPHERE+j1; - wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1; - wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1; - wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1; + if (!wsphere[cell].locked) + { + if (COLOR_BACKGROUND) + { + wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1*rgb_bg[0]; + wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1*rgb_bg[1]; + wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1*rgb_bg[2]; + } + else + { + wsphere[cell].r = particle[tracer_n[j]].rgb[0]*lum + lum1; + wsphere[cell].g = particle[tracer_n[j]].rgb[1]*lum + lum1; + wsphere[cell].b = particle[tracer_n[j]].rgb[2]*lum + lum1; + } + } } } } } -void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_cluster cluster[NMAXCIRCLES], t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], int traj_position, int traj_length, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot) +void draw_segments_sphere(t_segment segment[NMAXSEGMENTS], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE]) +/* draw the repelling segments on the sphere */ +{ + int s, i, cell, npoints, i0, j0, i1, j1, p, q; + double x1, y1, x2, y2, x, y, dt, length; + static double dphi, dtheta; + static int first = 1; + + if (first) + { + dphi = DPI/(double)NX_SPHERE; + dtheta = PI/(double)NY_SPHERE; + first = 0; + } + + for (s=0; s= NX_SPHERE) i1 = 0; + j1 = j0 + q; + if (j1 < 0) j1 = 0; + if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1; + cell = i1*NY_SPHERE+j1; + wsphere[cell].r = 0.0; + wsphere[cell].g = 0.0; + wsphere[cell].b = 0.0; + + wsphere[cell].radius += 0.005; + } + } + } +} + +void draw_absorbers_sphere(t_absorber absorber[NMAX_ABSORBERS], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE]) +/* draw the absorbing discs on the sphere */ +{ + int i, j, n, cell; + double x, y, dist; + static double dphi, dtheta; + static int first = 1; + + if (first) + { + dphi = DPI/(double)NX_SPHERE; + dtheta = PI/(double)NY_SPHERE; + first = 0; + } + + for (i=0; i 0)) { cell = hash_cell(wsphere[i].phi, wsphere[i].theta); wsphere[i].r = hashgrid[cell].r; @@ -1607,7 +1883,13 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA /* add tracer trajectories */ if (TRACER_PARTICLE) - draw_trajectory_sphere(trajectory, traj_position, traj_length, particle, cluster, wsphere, tracer_n, plot); + draw_trajectory_sphere(trajectory, hashgrid, traj_position, traj_length, particle, cluster, wsphere, tracer_n, plot); + + if (ADD_FIXED_SEGMENTS) + draw_segments_sphere(segment, wsphere); + +// if (ADD_ABSORBERS) +// draw_absorbers_sphere(absorber, wsphere); /* draw zero meridian, for debugging */ // for (j=0; j 0.0) + { + if (vabs(y) > 0.5*LAMBDA) return(1); + x1 = -y*y/LAMBDA + 0.25*LAMBDA; + if ((x < x1)||(x > x1 + WALL_WIDTH)) return(1); + if (vabs(y) < 0.1*LAMBDA) return(0); + return(2); + } + else + { + for (i=0; i polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(2); + if (vabs(y) > 0.5*MU) return(1); + x1 = y*y/MU - 0.25*MU; + if ((x > x1)||(x < x1 - WALL_WIDTH)) return(1); + return(2); + } + + } + case (D_TWO_PARABOLAS_LENS): + { + if (x > 0.0) + { + if (vabs(y) > 0.5*LAMBDA) return(1); + if (vabs(y) > 0.1*LAMBDA) + { + x1 = -y*y/LAMBDA + 0.25*LAMBDA; + if ((x < x1)||(x > x1 + WALL_WIDTH)) return(1); + return(2); + } + x1 = 0.24*LAMBDA; + x2 = x1 + WALL_WIDTH; + r2 = 1.5*0.25*LAMBDA; + r2 = r2*r2; + x3 = sqrt(r2 - 0.01*LAMBDA*LAMBDA) - x2; + +// r2 = x2*x2 + 0.01*LAMBDA*LAMBDA; +// if ((x*x + y*y < r2)&&((x-x1-x2)*(x-x1-x2) + y*y < r2)) return(0); + if (((x+x3)*(x+x3) + y*y < r2)&&(x > x1)) return(0); + +// if (vabs(y) < 0.1*LAMBDA) return(0); + return(1); + } + else + { + for (i=0; i polyrect[i].x1)&&(x < polyrect[i].x2)&&(y > polyrect[i].y1)&&(y < polyrect[i].y2)) return(2); + if (vabs(y) > 0.5*MU) return(1); + x1 = y*y/MU - 0.25*MU; + if ((x > x1)||(x < x1 - WALL_WIDTH)) return(1); + return(2); + } + + } + case (D_TWO_PARABOLAS_CLOSED_SEMITRANS): + { + if (vabs(y) > 0.75 + WALL_WIDTH) return(1); + x1 = 0.5*LAMBDA + 0.25*MU - y*y/MU; + if ((vabs(y) > 0.75)&&(vabs(x) < x1 + WALL_WIDTH)) return(2); + if ((vabs(x) < x1)||(vabs(x) > x1 + WALL_WIDTH)) return(1); + if ((x > 0.0)&&(vabs(y) < 0.2)) return(0); + return(2); + } + case (D_TWO_PARABOLAS_CLOSED_LENS): + { + if (vabs(y) > 0.75 + WALL_WIDTH) return(1); + if ((x > 0.0)&&(vabs(y) < 0.2)) + { + /* lens */ + y1 = 0.2; + x1 = 0.5*LAMBDA + 0.25*MU - y1*y1/MU; + x2 = x1 + WALL_WIDTH; + r2 = 1.5*0.25*LAMBDA; + r2 = r2*r2; + x3 = sqrt(r2 - y1*y1) - x2; + if (((x+x3)*(x+x3) + y*y < r2)&&(x > x1)) return(0); + return(1); + } + x1 = 0.5*LAMBDA + 0.25*MU - y*y/MU; + if ((vabs(y) > 0.75)&&(vabs(x) < x1 + WALL_WIDTH)) return(2); + if ((vabs(x) < x1)||(vabs(x) > x1 + WALL_WIDTH)) return(1); + return(2); + } + case (D_TWO_PARABOLAS_CLOSED_LENS_R): + { + if (vabs(y) > 0.75 + WALL_WIDTH) return(1); + if ((x > 0.0)&&(vabs(y) < WALL_WIDTH_B)) + { + /* lens */ + y1 = WALL_WIDTH_B; + x1 = 0.5*LAMBDA + 0.25*MU - y1*y1/MU; + x2 = x1 + WALL_WIDTH; + r2 = MU_B; + r2 = r2*r2; + x3 = sqrt(r2 - y1*y1) - x2; + if (((x+x3)*(x+x3) + y*y < r2)&&(x > x1)) return(0); + return(1); + } + x1 = 0.5*LAMBDA + 0.25*MU - y*y/MU; + if ((vabs(y) > 0.75)&&(vabs(x) < x1 + WALL_WIDTH)) return(2); + if ((vabs(x) < x1)||(vabs(x) > x1 + WALL_WIDTH)) return(1); + return(2); + } case (D_FOUR_PARABOLAS): { x1 = MU + LAMBDA - 0.25*y*y/MU; @@ -5812,6 +6102,13 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ condition += xy_in_triangle_tvertex(x, y, polyline[NPOLY], polyline[i], polyline[i+1]); return(condition >= 1); } + case (D_STAR_CHANNEL): + { + condition = ((x < 0.0)&&(vabs(y) < 0.5*WALL_WIDTH)); + for (i = 0; i < NPOLY+3; i++) + condition += xy_in_triangle_tvertex(x, y, polyline[NPOLY+3], polyline[i+1], polyline[i]); + return(condition >= 1); + } case (D_FRESNEL): { if (vabs(y) > 0.9*vabs(LAMBDA)) return(1); @@ -6352,6 +6649,45 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ return(1); } + case (D_TREE_BEAM): + { + /* upper triangle */ + h = 1.5*LAMBDA; + r = 1.5*LAMBDA; + x1 = vabs(x); + y1 = y - h; + zz[0][0] = r*cos(PI/6.0); zz[0][1] = -0.5*r; + zz[1][0] = 0.0; zz[1][1] = r; + zz[2][0] = 0.0; zz[2][1] = -0.5*r; + + /* middle triangle */ + h = 0.25*LAMBDA; + r = 2.0*LAMBDA; + y2 = y - h; + zz[3][0] = r*cos(PI/6.0); zz[3][1] = -0.5*r; + zz[4][0] = 0.0; zz[4][1] = r; + zz[5][0] = 0.0; zz[5][1] = -0.5*r; + + /* bottom triangle */ + h = -1.5*LAMBDA; + r = 3.0*LAMBDA; + y3 = y - h; + zz[6][0] = r*cos(PI/6.0); zz[6][1] = -0.5*r; + zz[7][0] = 0.0; zz[7][1] = r; + zz[8][0] = 0.0; zz[8][1] = -0.5*r; + + if (xy_in_triangle(x1, y1, zz[0], zz[1], zz[2])) return(1); + if (xy_in_triangle(x1, y2, zz[3], zz[4], zz[5])) return(1); + if (xy_in_triangle(x1, y3, zz[6], zz[7], zz[8])) return(1); + + if ((x < 0.0)&&(y > OSCIL_YMID - 0.5*WALL_WIDTH)&&(y < OSCIL_YMID + 0.5*WALL_WIDTH)) + return(1); + + if ((x < -zz[6][0] + 0.5*WALL_WIDTH)&&(y > OSCIL_YMID - 0.75*WALL_WIDTH)&&(y < OSCIL_YMID + 0.75*WALL_WIDTH)) + return(2); + + return(0); + } case (D_MICHELSON): { if ((vabs(x) > LAMBDA)&&(vabs(y) > LAMBDA)) return(2); @@ -6850,6 +7186,28 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_ if (vabs(y1) < 0.75*WALL_WIDTH) return(2); return(0); } + case (D_WAVEGUIDE_ELLIPSE_OBLIQUE): + { + if (x*x/(LAMBDA*LAMBDA) + y*y < 1.0) return(1); + if (x > 0.0) return(0.0); + y1 = y + 1.0 - MU - 0.5*WALL_WIDTH - (x - XMIN)*tan(APOLY*PID); + if (vabs(y1) < 0.5*WALL_WIDTH) return(1); + if (vabs(y1) < 0.75*WALL_WIDTH) return(2); + return(0); + } + case (D_CIRCLE_PAIRS): + { + dx = (XMAX - XMIN)/(double)NGRIDX; + for (k=0; k 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; phase = t - a*t*t; // if (time%1000 == 0) printf("time = %i, phase = %.4lg\n", time, phase); return(AMPLITUDE*cos(phase)*exp(-phase*DAMPING)); @@ -10882,6 +11646,8 @@ double oscillating_bc(int time, int j) case (OSC_WAVE_PACKET): { t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; a = sqrt(INITIAL_VARIANCE)/OMEGA; phase = AMPLITUDE*cos(t); envelope = exp(-(t-INITIAL_SHIFT)*(t-INITIAL_SHIFT)/(a*a))*sqrt(DPI/a); @@ -10890,6 +11656,8 @@ double oscillating_bc(int time, int j) case (OSC_WAVE_PACKETS): { t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; t1 = t - (double)((int)(t/WAVE_PACKET_SHIFT + 0.5))*WAVE_PACKET_SHIFT; a = sqrt(INITIAL_VARIANCE)/OMEGA; phase = AMPLITUDE*cos(t); @@ -10900,6 +11668,8 @@ double oscillating_bc(int time, int j) { // a = 0.25; t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; phase = t + ACHIRP*t*t; return(AMPLITUDE*sin(phase)*exp(-phase*DAMPING)); } @@ -10913,6 +11683,8 @@ double oscillating_bc(int time, int j) else dist2 = 0.0; dist2 = 500.0*dist2*dist2; t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING); amp *= exp(-dist2/INITIAL_VARIANCE); return(amp); @@ -10922,6 +11694,8 @@ double oscillating_bc(int time, int j) dist2 = (double)(j - NY/2)*(YMAX-YMIN)/(double)NY; dist2 = dist2*dist2; t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING); amp *= exp(-dist2/INITIAL_VARIANCE); return(amp); @@ -10937,6 +11711,8 @@ double oscillating_bc(int time, int j) if (j > jmax) return(0.0); if (j < jmin) return(0.0); t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; amp = AMPLITUDE*cos(t)*exp(-(double)t*DAMPING); amp *= cos(PID*dist2); return(amp); @@ -10952,6 +11728,8 @@ double oscillating_bc(int time, int j) if (j > jmax) return(0.0); if (j < jmin) return(0.0); t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; a = INITIAL_VARIANCE/(OMEGA*OMEGA); amp = AMPLITUDE*cos(t)*exp(-(double)(t-INITIAL_SHIFT)*(t-INITIAL_SHIFT)/a); amp *= cos(PID*dist2); @@ -10967,6 +11745,8 @@ double oscillating_bc(int time, int j) else dist2 = 0.0; dist2 = 500.0*dist2*dist2; t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; amp = AMPLITUDE*(0.5*cos(t) + 0.5*cos(sqrt(7.0)*t))*exp(-(double)t*DAMPING); amp *= exp(-dist2/INITIAL_VARIANCE); return(amp); @@ -10982,6 +11762,8 @@ double oscillating_bc(int time, int j) if (j > jmax) return(0.0); if (j < jmin) return(0.0); t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; amp = AMPLITUDE*(0.5*cos(t) + 0.5*cos(sqrt(7.0)*t))*exp(-(double)t*DAMPING); amp *= cos(PID*dist2); return(amp); @@ -10997,6 +11779,8 @@ double oscillating_bc(int time, int j) if (j > jmax) return(0.0); if (j < jmin) return(0.0); t = (double)time*OMEGA; + if (OSCIL_LEFT_YSHIFT > 0.0) + t -= (double)j*OSCIL_LEFT_YSHIFT/(double)NY; phase = t + ACHIRP*t*t; // if (phase > DPI) phase = DPI; amp = AMPLITUDE*sin(phase)*exp(-phase*DAMPING); diff --git a/sub_wave_3d.c b/sub_wave_3d.c index 16c53a0..4c618e4 100644 --- a/sub_wave_3d.c +++ b/sub_wave_3d.c @@ -1665,8 +1665,26 @@ void init_speed_dissipation(short int xy_in[NX*NY], double tc[NX*NY], double tcc } case (IOR_MICHELSON): { - printf("Case IOR_MICHELSON in init_speed_dissipation of sub_wave_3d needs to be updated\n"); - exit(1); + for (i=0; i zfloor)&&(*wave[(i+1)*NY+j].p_zfield[movie] > zfloor)&&(*wave[i*NY+j+1].p_zfield[movie] > zfloor)&&(*wave[(i+1)*NY+j+1].p_zfield[movie] > zfloor); diff --git a/wave_billiard.c b/wave_billiard.c index eaf5023..7233f8e 100644 --- a/wave_billiard.c +++ b/wave_billiard.c @@ -48,7 +48,7 @@ #define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ -#define VARIABLE_IOR 1 /* set to 1 for a variable index of refraction */ +#define VARIABLE_IOR 0 /* set to 1 for a variable index of refraction */ #define IOR 17 /* choice of index of refraction, see list in global_pdes.c */ #define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */ #define MANDEL_IOR_SCALE -0.05 /* parameter controlling dependence of IoR on Mandelbrot escape speed */ @@ -62,8 +62,8 @@ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ -#define YMIN -1.197916667 -#define YMAX 1.197916667 /* y interval for 9/16 aspect ratio */ +#define YMIN -1.297916667 +#define YMAX 1.097916667 /* y interval for 9/16 aspect ratio */ #define HIGHRES 1 /* set to 1 if resolution of grid is double that of displayed image */ #define HRES 1 /* dummy, only used by rde.c */ @@ -72,7 +72,7 @@ /* Choice of the billiard table */ -#define B_DOMAIN 982 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 421 /* choice of domain shape, see list in global_pdes.c */ #define CIRCLE_PATTERN 202 /* pattern of circles or polygons, see list in global_pdes.c */ #define IMAGE_FILE 5 /* for option D_IMAGE */ @@ -87,11 +87,11 @@ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ #define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ -#define LAMBDA 1.4 /* parameter controlling the dimensions of domain */ -#define MU 0.05 /* parameter controlling the dimensions of domain */ -#define MU_B 1.0 /* parameter controlling the dimensions of domain */ -#define NPOLY 3 /* number of sides of polygon */ -#define APOLY 1.0 /* angle by which to turn polygon, in units of Pi/2 */ +#define LAMBDA 0.9 /* parameter controlling the dimensions of domain */ +#define MU -0.2124612 /* parameter controlling the dimensions of domain */ +#define MU_B 0.3 /* parameter controlling the dimensions of domain */ +#define NPOLY 10 /* number of sides of polygon */ +#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 6 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ @@ -99,8 +99,8 @@ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ #define NGRIDX 6 /* number of grid point for grid of disks */ #define NGRIDY 14 /* number of grid point for grid of disks */ -#define WALL_WIDTH 0.1 /* width of wall separating lenses */ -#define WALL_WIDTH_B 0.05 /* width of wall separating lenses */ +#define WALL_WIDTH 0.075 /* width of wall separating lenses */ +#define WALL_WIDTH_B 0.1 /* width of wall separating lenses */ #define WALL_WIDTH_RND 0.0 /* proportion of width of width for random arrangements */ #define RADIUS_FACTOR 0.3 /* controls inner radius for C_RING arrangements */ #define WALL_WIDTH_ASYM 0.75 /* asymmetry of wall width (D_CIRCLE_LATTICE_NONISO) */ @@ -121,22 +121,22 @@ /* xy_in_billiard and draw_billiard below */ /* Physical parameters of wave equation */ -#define TWOSPEEDS 1 /* set to 1 to replace hardcore boundary by medium with different speed */ +#define TWOSPEEDS 0 /* set to 1 to replace hardcore boundary by medium with different speed */ #define OSCILLATE_LEFT 1 /* 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 42 /* oscillation schedule, see list in global_pdes.c */ -#define OSCIL_YMAX 0.05 /* defines oscilling beam range */ -#define OSCIL_YMID -0.9 /* defines oscilling beam midpoint */ +#define OSCIL_YMAX 0.0375 /* defines oscilling beam range */ +#define OSCIL_YMID 0.0 /* defines oscilling beam midpoint */ #define INITIAL_SHIFT 20.0 /* time shift of initial wave packet (in oscillation periods) */ #define WAVE_PACKET_SHIFT 200.0 /* time shift between wave packets (in oscillation periods) */ -#define OMEGA 0.01 /* frequency of periodic excitation */ +#define OMEGA 0.006 /* frequency of periodic excitation */ #define AMPLITUDE 2.0 /* amplitude of periodic excitation */ #define ACHIRP 0.25 /* acceleration coefficient in chirp */ #define DAMPING 0.0 /* damping of periodic excitation */ -#define COURANT 0.25 /* Courant number in medium B */ -#define COURANTB 0.1 /* Courant number */ -#define GAMMA 5.0e-6 /* damping factor in wave equation */ +#define COURANT 0.08 /* Courant number in medium B */ +#define COURANTB 0.25 /* Courant number */ +#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 */ #define GAMMA_TOPBOT 1.0e-7 /* damping factor on boundary */ @@ -150,11 +150,11 @@ /* For similar wave forms, COURANT^2*GAMMA should be kept constant */ #define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */ -#define OSCILLATING_SOURCE_PERIOD 25.0 /* period of oscillating source */ +#define OSCILLATING_SOURCE_PERIOD 6.0 /* period of oscillating source */ #define ALTERNATE_OSCILLATING_SOURCE 1 /* set to 1 to alternate sign of oscillating source */ #define N_SOURCES 1 /* number of sources, for option draw_sources */ #define ALTERNATE_SOURCE_PHASES 0 /* set to 1 to alternate initial phases of sources */ -#define MAX_PULSING_TIME 750 /* max time for adding pulses */ +#define MAX_PULSING_TIME 5000 /* max time for adding pulses */ #define ADD_WAVE_PACKET_SOURCES 0 /* set to 1 to add several sources emitting wave packets */ #define WAVE_PACKET_SOURCE_TYPE 3 /* type of wave packet sources */ @@ -169,10 +169,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2500 /* number of frames of movie */ +#define NSTEPS 2700 /* number of frames of movie */ #define NVID 12 /* number of iterations between images displayed on screen */ #define NSEG 1000 /* number of segments of boundary */ -#define INITIAL_TIME 150 /* time after which to start saving frames */ +#define INITIAL_TIME 100 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PRINT_SPEED 0 /* print speed of moving source */ #define PRINT_FREQUENCY 0 /* print frequency (for phased array) */ @@ -199,8 +199,8 @@ /* 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 18 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 18 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* background */ @@ -208,17 +208,17 @@ #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 COLOR_RANGE 1.0 /* max range of color (default: 1.0) */ +#define COLOR_RANGE 0.75 /* max range of color (default: 1.0) */ #define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define VSHIFT_AMPLITUDE -0.1 /* additional shift for wave amplitude */ +#define VSHIFT_AMPLITUDE -0.5 /* additional shift for wave amplitude */ #define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for wave amplitude */ -#define E_SCALE 50.0 /* scaling factor for energy representation */ +#define E_SCALE 300.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 AVRG_E_FACTOR 0.95 /* controls time window size in P_AVERAGE_ENERGY scheme */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ #define FADE_IN_OBSTACLE 1 /* set to 1 to fade color inside obstacles */ #define SHADE_2D 1 /* set to 1 to add pseudo-3d shading effect */ @@ -233,16 +233,16 @@ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ #define COLORBAR_RANGE 1.7 /* scale of color scheme bar */ -#define COLORBAR_RANGE_B 0.8 /* scale of color scheme bar for 2nd part */ +#define COLORBAR_RANGE_B 4.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 */ #define CIRC_COLORBAR_B 0 /* set to 1 to draw circular color scheme */ -#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 DRAW_WAVE_PROFILE 1 /* set to 1 to draw a profile of the wave */ +#define HORIZONTAL_WAVE_PROFILE 1 /* set to 1 to draw wave profile vertically */ +#define VERTICAL_WAVE_PROFILE 0 /* 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.12 /* value of y to sample wave profile */ +#define WAVE_PROFILE_Y 0.0 /* value of y to sample wave profile */ #define PROFILE_AT_BOTTOM 1 /* draw wave profile at bottom instead of top */ #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 */ @@ -691,9 +691,9 @@ void animation() /* initialise polyline and similar for drawing some domains */ npolyline = init_poly(MDEPTH, polyline, polyrect, polyrectrot, polyarc, circles, &npolyrect, &npolyrect_rot, &npolyarc, &ncircles, 1); - for (i=0; i