diff --git a/global_3d.c b/global_3d.c index bdcea5f..aa842eb 100644 --- a/global_3d.c +++ b/global_3d.c @@ -16,6 +16,8 @@ #define E_KURAMOTO_SPHERE 91 /* Kuramoto on the sphere, with cubic grid */ #define E_KELLER_SEGEL 10 /* Keller-Segel model */ #define E_KELLER_SEGEL_SPHERE 11 /* Keller-Segel model on the sphere, cubic grid */ +#define E_GRAY_SCOTT 12 /* Gray-Scott model */ +#define E_GRAY_SCOTT_SPHERE 121 /* Gray-Scott model on the sphere, cubic grid */ /* Choice of potential */ @@ -102,7 +104,7 @@ #define RDE_PLANET (((ADAPT_STATE_TO_BC)&&((OBSTACLE_GEOMETRY == D_SPHERE_EARTH)||(OBSTACLE_GEOMETRY == D_SPHERE_MARS)||(OBSTACLE_GEOMETRY == D_SPHERE_VENUS)))) // #define RDE_PLANET (((ADAPT_STATE_TO_BC)&&((OBSTACLE_GEOMETRY == D_SPHERE_EARTH)||(OBSTACLE_GEOMETRY == D_SPHERE_MARS)||(OBSTACLE_GEOMETRY == D_SPHERE_VENUS)))||((RDE_EQUATION == E_SHALLOW_WATER)&&(SWATER_DEPTH == SH_EARTH))) -#define NMAXCIRC_SPHERE 100 /* max number of circles on sphere */ +#define NMAXCIRC_SPHERE 1000 /* max number of circles on sphere */ #define NMAX_TRACER_PTS 20 /* max number of tracer points recorded per cell */ #define NMAX_SPHERE_NEIGHB 10 /* max number of neighbours in Poisson simulation grid */ @@ -157,6 +159,7 @@ typedef struct double depth; /* water depth */ double cos_depth_angle; /* cos of angle of depth profile */ double gradx, grady; /* gradient of water depth */ + double time; /* time when threshold is reached */ short int tracer; /* has value 1 if cell contains a tracer */ short int n_tracer_pts; /* number of recorded tracer points per cell */ double tracerx[NMAX_TRACER_PTS], tracery[NMAX_TRACER_PTS]; /* coordinates of tracer */ @@ -195,6 +198,7 @@ typedef struct 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 */ + double laplace_weight; /* weight for adjusted Laplacian */ } t_wave_sphere; diff --git a/global_ljones.c b/global_ljones.c index 886c8f6..e636c29 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -15,7 +15,7 @@ #define NMAXCIRCLES 200000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define MAXNEIGH 20 /* max number of neighbours kept in memory */ -#define NMAXOBSTACLES 5000 /* max number of obstacles */ +#define NMAXOBSTACLES 10000 /* max number of obstacles */ #define NMAXSEGMENTS 1200 /* max number of repelling segments */ #define NMAXGROUPS 100 /* max number of groups of segments */ #define NMAXCOLLISIONS 200000 /* max number of collisions */ @@ -26,7 +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 NMAX_ABSORBERS 100 /* max number of absorbing discs */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -88,6 +88,7 @@ #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_VERTICAL 111 /* vertical double line */ #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 */ @@ -347,6 +348,12 @@ #define IP_X 0 /* color depends on x coordinate of initial position */ #define IP_Y 1 /* color depends on y coordinate of initial position */ +/* Position-dependent initial conditions (for POSITION_DEPENDENT_TYPE) */ + +#define PDIC_TWO 1 /* two layers */ +#define PDIC_THREE 2 /* three layers */ +#define PDIC_MASSY 3 /* mass depends on y coordinate*/ + /* Space dependence of electric field */ #define EF_CONST 0 /* constant magnetic field */ @@ -396,6 +403,7 @@ #define BG_POS_OBSTACLES 8 /* background color depends on displacement of obstacles */ #define BG_CURRENTX 9 /* background color depends on x-component of current */ #define BG_DIRECTION 10 /* background color depends on particle orientation */ +#define BG_POTENTIAL 12 /* potential added to sphere */ /* Obstacle color schemes */ @@ -422,6 +430,18 @@ #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 */ +#define AP_PDISC 9 /* Poisson disc process */ +#define AP_POLES 10 /* absorbers at the North and South pole */ + +/* Type of sphere potential */ + +#define SPP_WELLS 0 /* potential wells located at points in table wells[] */ +#define SPP_CUBE 1 /* potential having shape of a cube */ + +/* Type of sphere gravity */ + +#define SG_POLAR 1 /* gravity moving from pole to pole */ +#define SG_EQUATOR 2 /* gravity moving along equator */ /* Color schemes */ @@ -706,6 +726,7 @@ typedef struct double fboundary; /* boundary force */ double pressure; /* pressure */ double gravity; /* gravity */ + double g_angle; /* angle of gravity (for sphere) */ double radius; /* particle radius */ double angle; /* orientation of obstacle */ double omega; /* angular speed of obstacle */ @@ -714,6 +735,7 @@ typedef struct double prop; /* proportion of types */ double thermo_radius; /* radius of thermostat region */ double current; /* electrical current */ + double omega_sphere; /* angular frequency of rotating sphere */ } t_lj_parameters; @@ -725,12 +747,16 @@ 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 initial_radius; /* radius before particles are added */ double r, g, b; /* RGB values for image */ double cos_angle; /* cosine of light angle */ double cos_angle_sphere; /* cosine of light angle for perfect sphere */ + double cos_angle_pot; /* cosine of light angle for 2D potential */ + double potential; /* value of potential */ + double nablax, nablay; /* gradient of potential */ 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, nabsorbers = 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, nwells = 0; FILE *lj_log; diff --git a/global_pdes.c b/global_pdes.c index 65b956d..7866c28 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -146,6 +146,7 @@ #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 */ +#define D_CIRCLE_POLYGON 991 /* circles connected to polygons by channels */ /* for wave_sphere.c */ @@ -208,6 +209,8 @@ #define C_SPH_ICOSA 31 /* icosahedron (on sphere) */ #define C_SPH_OCTA 32 /* octahedron (on sphere) */ #define C_SPH_CUBE 33 /* cube (on sphere) */ +#define C_SPH_CUBE_B 34 /* rotated cube (on sphere) */ +#define C_SPH_POISSON 35 /* Poisson process (on sphere) */ #define C_ONE 97 /* one single circle, as for Sinai */ #define C_TWO 98 /* two concentric circles of different type */ @@ -415,6 +418,9 @@ #define Z_SWATER_SPEED 71 /* speed */ #define Z_SWATER_DIRECTION_SPEED 74 /* hue for direction of velocity, luminosity for speed */ +/* for Gray-Scott model */ +#define Z_TIME_SLICE 80 /* first time field 0 reaches threshold */ + /* special boundary conditions for Euler equation */ #define BCE_TOPBOTTOM 1 /* special flow at top and bottom */ #define BCE_TOPBOTTOMLEFT 2 /* special flow at top, bottom and left side */ diff --git a/lennardjones.c b/lennardjones.c index 4130157..12a2102 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -58,18 +58,18 @@ #define YMIN 0.0 #define YMAX 3.141592654 /* y interval for 9/16 aspect ratio */ -#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 INITXMIN 0.0 +#define INITXMAX 6.28 /* x interval for initial condition */ +#define INITYMIN 0.0 +#define INITYMAX 3.14 /* 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 0.0 -#define ADDXMAX 0.1 /* x interval for adding particles */ +#define ADDXMIN 0.1 +#define ADDXMAX 0.2 /* x interval for adding particles */ #define ADDYMIN 1.57 #define ADDYMAX 1.57 /* y interval for adding particles */ #define ADDRMIN 2.0 @@ -85,7 +85,7 @@ #define OBSYMIN -1.125 #define OBSYMAX 1.125 /* x interval for motion of obstacle */ -#define CIRCLE_PATTERN 20 /* pattern of circles, see list in global_ljones.c */ +#define CIRCLE_PATTERN 81 /* 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 */ @@ -108,7 +108,7 @@ #define COUPLE_MINLENGTH 0.5 /* length at which bonds decouple */ #define ADD_FIXED_SEGMENTS 0 /* set to 1 to add fixed segments as obstacles */ -#define SEGMENT_PATTERN 100 /* pattern of repelling segments, see list in global_ljones.c */ +#define SEGMENT_PATTERN 111 /* pattern of repelling segments, see list in global_ljones.c */ #define ROCKET_SHAPE 3 /* shape of rocket combustion chamber, see list in global_ljones.c */ #define ROCKET_SHAPE_B 3 /* shape of second rocket */ #define NOZZLE_SHAPE 6 /* shape of nozzle, see list in global_ljones.c */ @@ -126,21 +126,21 @@ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ #define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ -#define INTERACTION 1 /* particle interaction, see list in global_ljones.c */ -#define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */ +#define INTERACTION 1 /* particle interaction, see list in global_ljones.c */ +#define INTERACTION_B 1 /* particle interaction for second type of particle, see list in global_ljones.c */ #define SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */ #define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */ #define MOL_ANGLE_FACTOR 1.0 /* rotation angle for P_MOL_ANGLE color scheme */ #define P_PERCOL 0.25 /* probability of having a circle in C_RAND_PERCOL arrangement */ #define NPOISSON 100 /* number of points for Poisson C_RAND_POISSON arrangement */ -#define PDISC_DISTANCE 5.0 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 1.8 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ -#define MU 0.02 /* parameter controlling radius of particles */ -#define MU_B 0.02 /* parameter controlling radius of particles of second type */ +#define MU 0.018 /* parameter controlling radius of particles */ +#define MU_B 0.018 /* 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 */ @@ -171,7 +171,7 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 3800 /* number of frames of movie */ +#define NSTEPS 1750 /* 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 */ @@ -193,15 +193,17 @@ /* Plot type, see list in global_ljones.c */ -#define PLOT 13 +#define PLOT 14 #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 1 /* type of background coloring, see list in global_ljones.c */ +#define BG_COLOR 2 /* 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 */ +#define SHADE_BG_COLOR_2D 1 /* set to 1 to shade BG color, for option BG_POTENTIAL */ +#define SHADE_SCALE_BG_2D 0.1 /* controls 2D shading */ #define DRAW_BONDS 0 /* set to 1 to draw bonds between neighbours */ #define COLOR_BONDS 1 /* set to 1 to color bonds according to length */ @@ -212,7 +214,7 @@ #define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */ #define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */ #define N_PARTICLE_COLORS 5 /* number of colors for P_NUMBER color scheme */ -#define INITIAL_POS_TYPE 0 /* type of initial position dependence */ +#define INITIAL_POS_TYPE 1 /* type of initial position dependence */ #define ERATIO 0.995 /* ratio for time-averaging in P_EMEAN color scheme */ #define DRATIO 0.999 /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */ #define OBSTACLE_AREA_SHADE_FACTOR 75.0 /* controls sensitivity of triangle shade for option FILL_OBSTACLE_TRIANGLES */ @@ -222,7 +224,7 @@ #define COLOR_PALETTE 10 /* Color palette, see list in global_ljones.c */ #define COLOR_PALETTE_EKIN 10 /* Color palette for kinetic energy */ -#define COLOR_PALETTE_ANGLE 0 /* Color palette for angle representation */ +#define COLOR_PALETTE_ANGLE 17 /* Color palette for angle representation */ #define COLOR_PALETTE_DIRECTION 0 /* Color palette for direction representation */ #define COLOR_PALETTE_INITIAL_POS 10 /* Color palette for initial position representation */ #define COLOR_PALETTE_DIFFNEIGH 10 /* Color palette for different neighbours representation */ @@ -262,17 +264,17 @@ #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 PRINT_NABSORBED 0 /* print number of absorbed particles */ #define FORCE_FACTOR 0.1 /* factor controlling length of force vector */ /* particle properties */ -#define ENERGY_HUE_MIN 330.0 /* color of original particle */ +#define ENERGY_HUE_MIN 350.0 /* color of original particle */ #define ENERGY_HUE_MAX 50.0 /* color of saturated particle */ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ #define PARTICLE_EMIN 0.0 /* energy of particle with coolest color */ -#define PARTICLE_EMAX 250000.0 /* energy of particle with hottest color */ +#define PARTICLE_EMAX 4000.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 */ @@ -280,7 +282,7 @@ #define OBSTACLE_EMAX 1000000.0 /* energy of obstacle with hottest color */ #define OBSTACLE_VMAX 4.0 /* speed of obstacle with largest luminosity */ #define HUE_TYPE0 320.0 /* hue of particles of type 0 */ -#define HUE_TYPE1 60.0 /* hue of particles of type 1 */ +#define HUE_TYPE1 40.0 /* hue of particles of type 1 */ #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 180.0 /* hue of particles of type 4 */ @@ -291,6 +293,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 1.0 /* constant in BG_CHARGE backgound color scheme (default: 0.5) */ +#define BG_POTENTIAL_SLOPE 50.0 /* constant in BG_POTENTIAL background color scheme */ #define CHARGE_HUE_RANGE 0.5 /* range of charge colors */ #define PARTICLE_LMAX 1.5e4 /* angular momentum particle with brightest color */ @@ -300,14 +303,14 @@ #define ADAPT_MASS_TO_RADIUS 0 /* set to positive value to for mass prop to power of radius */ #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 DT_PARTICLE 2.0e-6 /* time step for particle displacement */ #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 EQUILIBRIUM_DIST 3.5 /* Lennard-Jones equilibrium distance */ +#define EQUILIBRIUM_DIST_B 4.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 25.0 /* damping coefficient of particles */ +#define INITIAL_DAMPING 0.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 2.0 /* mass of particle of radius MU */ @@ -316,39 +319,43 @@ #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 0.0 /* initial velocity range */ -#define V_INITIAL_ADD 5500.0 /* initial velocity range for added particles */ +#define V_INITIAL 50.0 /* initial velocity range */ +#define V_INITIAL_ADD 4500.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 500.0 /* relative intensity of Coulomb interaction in I_COULOMB_LJ (default: 100.0) */ +#define COULOMB_ALWAYS_REPEL 1 /* set to 1 to always repel with I_COULOMB_IMAGINARY */ #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 0 /* set to 1 to switch on thermostat */ +#define THERMOSTAT 1 /* set to 1 to switch on thermostat */ #define VARY_THERMOSTAT 0 /* set to 1 for time-dependent thermostat schedule */ #define SIGMA 5.0 /* noise intensity in thermostat */ -#define BETA 0.00005 /* initial inverse temperature */ +#define BETA 0.002 /* 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 */ -#define NBH_DIST_FACTOR 6.0 /* radius in which to count neighbours */ +#define NBH_DIST_FACTOR 6.0 /* radius in which to count neighbours */ #define BOND_DIST_FACTOR 6.0 /* radius in which to draw bonds */ #define GRAVITY 0.0 /* gravity acting on all particles */ #define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */ #define CIRCULAR_GRAVITY 0 /* set to 1 to have gravity directed to center */ +#define SPHERE_GRAVITY 0 /* set to 1 to have gravity at constant angle wrt sphere */ #define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */ #define GRAVITY_SCHEDULE 1 /* type of gravity schedule, see list in global_ljones.c */ -#define GRAVITY_FACTOR 10.0 /* factor by which to increase gravity */ -#define GRAVITY_INITIAL_TIME 250 /* time at start of simulation with constant gravity */ -#define GRAVITY_RESTORE_TIME 500 /* time at end of simulation with gravity restored to initial value */ +#define GRAVITY_FACTOR 2000.0 /* factor by which to increase gravity */ +#define GRAVITY_INITIAL_TIME 100 /* time at start of simulation with constant gravity */ +#define GRAVITY_RESTORE_TIME 750 /* time at end of simulation with gravity restored to initial value */ +#define GRAVITY_INITIAL_ANGLE 0.0 /* initial angle for SPHERE_GRAVITY */ +#define GRAVITY_DELTA_ANGLE 1440.0 /* increase of angle for SPHERE_GRAVITY */ #define KSPRING_VICSEK 0.2 /* spring constant for I_VICSEK_SPEED interaction */ #define VICSEK_REPULSION 10.0 /* repulsion between particles in Vicsek model */ #define ADD_EFIELD 0 /* set to 1 to add an electric field */ -#define EFIELD 150000.0 /* value of electric field */ +#define EFIELD 20000.0 /* value of electric field */ #define EFIELD_Y 0.0 /* value of electric field */ #define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ #define BFIELD 20000.0 /* value of magnetic field */ @@ -359,7 +366,7 @@ #define INCREASE_E 0 /* set to 1 to increase electric field */ #define OSCILLATE_E 0 /* set to 1 for oscillating electric field */ #define E_PERIOD 1000 /* period of oscillating electric field */ -#define EFIELD_FACTOR 180000.0 /* factor by which to increase electric field */ +#define EFIELD_FACTOR 1000.0 /* factor by which to increase electric field */ #define INCREASE_B 0 /* set to 1 to increase magnetic field */ #define BFIELD_FACTOR 1000.0 /* factor by which to increase magnetic field */ #define CHARGE_OBSTACLES 1 /* set to 1 for obstacles to be charged */ @@ -382,6 +389,11 @@ #define WIND_FORCE 1.35e6 /* force of wind */ #define WIND_YMIN -0.6 /* min altitude of region with wind */ +#define ROTATE_SPHERE 1 /* set to 1 to add Coriolis and centripetal force */ +#define OMEGA_SPHERE 4.0 /* angular frequency of rotating sphere */ +#define CHANGE_OMEGA_SPHERE 1 /* set to 1 to change sphere rotation frequency */ +#define OMEGA_SPHERE_FACTOR 200.0 /* change factor of sphere rotation frequency */ + #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 */ @@ -398,15 +410,15 @@ #define QUADRUPOLE_RATIO 0.6 /* anisotropy in quadrupole potential */ #define INCREASE_BETA 0 /* set to 1 to increase 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 BETA_SCHEDULE 3 /* type of temperature schedule, see TS_* in global_ljones */ +#define BETA_FACTOR 0.002 /* factor by which to change BETA during simulation */ #define TS_SLOPE 8.5 /* controls speed of change of BETA for TS_TANH schedule (default 1.0) */ #define N_TOSCILLATIONS 1.0 /* number of temperature oscillations in BETA schedule */ #define NO_OSCILLATION 0 /* set to 1 to have exponential BETA change only */ #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 */ +#define FINAL_CONSTANT_PHASE 400 /* final phase in which temperature is constant */ #define DECREASE_CONTAINER_SIZE 0 /* set to 1 to decrease size of container */ #define SMOOTH_CONTAINER_DECREASE 1 /* set to 1 to decrease size smoothly at each simulation step */ @@ -429,10 +441,10 @@ #define N_T_AVERAGE 1 /* size of temperature averaging window */ #define MAX_PRESSURE 3.0e10 /* pressure shown in "hottest" color */ #define PARTIAL_THERMO_COUPLING 0 /* set to 1 to couple only some particles to thermostat */ -#define PARTIAL_THERMO_REGION 9 /* region for partial thermostat coupling (see list in global_ljones.c) */ +#define PARTIAL_THERMO_REGION 2 /* region for partial thermostat coupling (see list in global_ljones.c) */ #define PARTIAL_THERMO_SHIFT 0.2 /* distance from obstacle at the right of which particles are coupled to thermostat */ -#define PARTIAL_THERMO_WIDTH 1.0 /* vertical size of partial thermostat coupling */ -#define PARTIAL_THERMO_HEIGHT 0.0 /* vertical size of partial thermostat coupling */ +#define PARTIAL_THERMO_WIDTH 0.3 /* vertical size of partial thermostat coupling */ +#define PARTIAL_THERMO_HEIGHT 0.2 /* vertical size of partial thermostat coupling */ #define PARTIAL_THERMO_RIN 0.5 /* initial radius of region without coupling */ #define PARTIAL_THERMO_RFIN 1.3 /* final radius of region without coupling */ @@ -443,20 +455,20 @@ #define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */ #define NPART_BOTTOM 100 /* number of particles at the bottom */ -#define ADD_PARTICLES 1 /* set to 1 to add particles */ +#define ADD_PARTICLES 0 /* 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 FINAL_NOADD_PERIOD 100 /* 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 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 1000 /* number of tracer particles */ +#define TRACER_PARTICLE 0 /* set to 1 to have a tracer particle */ +#define N_TRACER_PARTICLES 5500 /* 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 1000 /* length of drawn trajectory */ @@ -506,10 +518,12 @@ #define TRACK_SEGMENT_GROUPS 0 /* set to 1 for view to track group of segments */ #define TRACK_X_PADDING 2.0 /* distance from x boundary where tracking starts */ -#define POSITION_DEPENDENT_TYPE 0 /* set to 1 to make particle type depend on initial position */ -#define POSITION_Y_DEPENDENCE 0 /* set to 1 for the separation between particles to be horizontal */ +#define POSITION_DEPENDENT_TYPE 0 /* set to PDIC_* to make particle type depend on initial position */ +#define POSITION_Y_DEPENDENCE 1 /* set to 1 for the separation between particles to be horizontal */ #define POSITION_DEP_SIGN -1.0 /* sign in position dependence condition */ -#define POSITION_DEP_X -0.625 /* threshold value for position-dependent type */ +#define POSITION_DEP_X 1.5 /* threshold value for position-dependent type */ +#define POSITION_DEP_Y 1.5 /* threshold value for position-dependent type */ +#define POSITION_DEP_MASS_RATIO 5.0 /* position-dependent mass factor */ #define PRINT_ENTROPY 0 /* set to 1 to compute entropy */ #define SPECIAL_IC 0 /* set to 1 for choosing special initial condition RD_INITIAL_COND */ @@ -602,11 +616,21 @@ #define PARTICLE_MASS_D 1.0 /* mass or 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 ADD_ABSORBERS 0 /* set to 1 to add absorbing discs */ +#define ABSORBER_PATTERN 3 /* 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 ABSORBER_R 0.015 /* radius of absorber */ +#define ABSORBER_PDIST 0.4 /* parameter of Poisson disc process */ + +#define ADD_POTENTIAL_SPHERE 0 /* add potential for gradient force field on sphere */ +#define DRAW_POTENTIAL_SPHERE 1 /* draw sphere radius depending on potential */ +#define SPHERE_POTENTIAL 1 /* type of sphere potential */ +#define SPHERE_POT_PATTERN 3 /* pattern of local minma of SPP_WELLS sphere potential */ +#define POT_SPHERE_AMP 1.0 /* amplitude in definition of potential on sphere */ +#define POT_SPHERE_RADIUS 0.1 /* radius in definition of potential on sphere */ +#define POT_SPHERE_SMOOTH 0.5 /* smoothing of potential well */ +#define POT_SPHERE_STRENGTH 1.0e3 /* coefficient of gradient force */ #define NXMAZE 16 /* width of maze */ #define NYMAZE 16 /* height of maze */ @@ -620,8 +644,8 @@ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 60 /* size of hashgrid in x direction */ -#define HASHY 30 /* size of hashgrid in y direction */ +#define HASHX 120 /* size of hashgrid in x direction */ +#define HASHY 60 /* size of hashgrid in y direction */ #define HASHMAX 100 /* maximal number of particles per hashgrid cell */ #define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ @@ -634,27 +658,28 @@ /* constants related to evolution on a sphere */ #define SPHERE 1 /* set to 1 to compute evolution in spherical geometry */ -#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 SIN_THETA_REG 0.01 /* regularization of sin(theta) for motion on sphere */ +#define POLAR_PADDING 0.01 /* region around poles that belong to the same hashcell */ #define DRAW_SPHERE 1 /* set to 1 to draw 3D sphere */ +#define DRAW_ELLIPSES_ON_SPHERE 1 /* set to 1 to draw ellipses instead of circles on sphere in 2D */ #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 XY_SCALING_FACTOR 2.0 /* 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 RSCALE_POTENTIAL 1.0 /* radial scaling of potential */ #define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ -#define ROTATE_ANGLE 750.0 /* total angle of rotation during simulation */ +#define ROTATE_ANGLE 360.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, 0.0, 1.5}; /* location of observer for REP_PROJ_3D representation */ - +double observer[3] = {-3.0, 2.0, 2.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)) #define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE)) @@ -992,6 +1017,18 @@ double gravity_schedule(int i, int j) } } + +double g_angle_schedule(int i, int j) +{ + double time; + + if (i < INITIAL_TIME + GRAVITY_INITIAL_TIME) return(0.0); + time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_INITIAL_TIME); + return(GRAVITY_INITIAL_ANGLE*PI/180.0 + time*GRAVITY_DELTA_ANGLE*PI/180.0); +// return(PI*time); +} + + double rotation_angle(double phase) { /* case of rotating hourglass */ @@ -1851,13 +1888,11 @@ void evolve_obstacles(t_obstacle obstacle[NMAXOBSTACLES]) void animation() { double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, - a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, angle, theta, sum, alpha, bfield, track_x0, track_y0, efield, efieldy, dist; + a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, angle, theta, sum, alpha, bfield, track_x0, track_y0, efield, efieldy, dist, dphi_inv, dtheta_inv, sincos; double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds, *currents; double *cqx, *cqy, *cpx, *cpy, *cqangle, *cpangle; - int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, cl, - min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, - tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0, - group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt, obs, reset_obs; + int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, + total_neighbours = 0, cl, min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0, group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt, obs, reset_obs; int *particle_numbers; static int imin, imax; static short int first = 1; @@ -1876,7 +1911,7 @@ void animation() t_otriangle *otriangle; t_ofacet *ofacet; t_lj_sphere *wsphere; - t_absorber *absorber; + t_absorber *absorber, *wells; char message[100]; ratioc = 1.0 - ratio; @@ -1888,16 +1923,19 @@ void animation() params.xmaxcontainer = BCXMAX; params.fboundary = 0.0; params.gravity = GRAVITY; + params.g_angle = 0.0; params.radius = MU; params.nabsorbed = 0; particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle)); /* particles */ /* sphere data */ - if (DRAW_SPHERE) + if (SPHERE) { wsphere = (t_lj_sphere *)malloc(NX_SPHERE*NY_SPHERE*sizeof(t_lj_sphere)); init_lj_sphere(wsphere); + dphi_inv = (double)NX_SPHERE/DPI; + dtheta_inv = (double)NY_SPHERE/PI; } if (CLUSTER_PARTICLES) @@ -1949,10 +1987,20 @@ void animation() for (i=0; i<2*NMAXCOLLISIONS; i++) collisions[i].time = 0; } + if (ADD_POTENTIAL_SPHERE) + { + if (SPHERE_POTENTIAL == SPP_WELLS) + { + wells = (t_absorber *)malloc(2*NMAX_ABSORBERS*sizeof(t_absorber)); + nwells = init_absorbers(wells, SPHERE_POT_PATTERN, POT_SPHERE_RADIUS); + } + init_potential_sphere(wsphere, wells); + } + if (ADD_ABSORBERS) { absorber = (t_absorber *)malloc(2*NMAX_ABSORBERS*sizeof(t_absorber)); - nabsorbers = init_absorbers(absorber); + nabsorbers = init_absorbers(absorber, ABSORBER_PATTERN, ABSORBER_R); if (DRAW_SPHERE) draw_absorbers_sphere(absorber, wsphere); } @@ -1965,7 +2013,11 @@ void animation() lj_log = fopen("lj_logfile.txt", "w"); if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle, otriangle, ofacet); - if (ADD_FIXED_SEGMENTS) init_segment_config(segment, conveyor_belt); + if (ADD_FIXED_SEGMENTS) + { + init_segment_config(segment, conveyor_belt); + if (DRAW_SPHERE) draw_segments_sphere(segment, wsphere); + } if ((MOVE_SEGMENT_GROUPS)&&(ADD_FIXED_SEGMENTS)) { @@ -2050,6 +2102,8 @@ void animation() else params.efield = EFIELD; if (INCREASE_B) params.bfield = bfield_schedule(i); else params.bfield = BFIELD; + if (CHANGE_OMEGA_SPHERE) params.omega_sphere = omega_sphere_schedule(i); + else params.omega_sphere = OMEGA_SPHERE; if ((PARTIAL_THERMO_COUPLING)&&(PARTIAL_THERMO_REGION == TH_RING_EXPAND)) params.thermo_radius = PARTIAL_THERMO_RIN + (double)i/(double)NSTEPS*(PARTIAL_THERMO_RFIN - PARTIAL_THERMO_RIN); if (DECREASE_CONTAINER_SIZE) @@ -2137,6 +2191,7 @@ void animation() if (INCREASE_GRAVITY) params.gravity = gravity_schedule(i,n); + if (SPHERE_GRAVITY > 0) params.g_angle = g_angle_schedule(i,n); if ((BOUNDARY_COND == BC_RECTANGLE_WALL)&&(i < INITIAL_TIME + WALL_TIME)) wall = 1; else wall = 0; @@ -2244,7 +2299,18 @@ void animation() // } /* add gravity */ - if (INCREASE_GRAVITY) + if (SPHERE_GRAVITY == SG_POLAR) + { + particle[j].fx += params.gravity*sin(params.g_angle)*sin(particle[j].xc)/particle[j].mass_inv; + particle[j].fy -= params.gravity*sin(params.g_angle)*cos(particle[j].xc)*cos(particle[j].yc)/particle[j].mass_inv; + particle[j].fy -= params.gravity*cos(params.g_angle)*sin(particle[j].yc)/particle[j].mass_inv; + } + else if (SPHERE_GRAVITY == SG_EQUATOR) + { + particle[j].fx += params.gravity*sin(params.g_angle - particle[j].xc)/particle[j].mass_inv; + particle[j].fy += params.gravity*cos(params.g_angle - particle[j].xc)*cos(particle[j].yc)/particle[j].mass_inv; + } + else if (INCREASE_GRAVITY) { if (CIRCULAR_GRAVITY) { @@ -2268,7 +2334,11 @@ void animation() if (ADD_EFIELD) { if ((INCREASE_E)||(OSCILLATE_E)) - particle[j].fx += params.efield*particle[j].charge; + { + if (SPHERE) /* zero E field at poles */ + particle[j].fx += params.efield*particle[j].charge*sin(particle[j].yc); + else particle[j].fx += params.efield*particle[j].charge; + } else { efield = EFIELD; @@ -2278,7 +2348,9 @@ void animation() efield *= (double)partial_efield(particle[j].xc, particle[j].yc); efieldy *= (double)partial_efield(particle[j].xc, particle[j].yc); } - particle[j].fx += efield*particle[j].charge; + if (SPHERE) /* zero E field at poles */ + particle[j].fx += efield*particle[j].charge*sin(particle[j].yc); + else particle[j].fx += efield*particle[j].charge; particle[j].fy += efieldy*particle[j].charge; } } @@ -2299,6 +2371,26 @@ void animation() particle[j].fx += WIND_FORCE*particle[j].radius*particle[j].mass_inv; } + if ((SPHERE)&&(ROTATE_SPHERE)) + { + /* centrifugal force */ + sincos = cos(particle[j].yc)*sin(particle[j].yc); + particle[j].fy += params.omega_sphere*params.omega_sphere*sincos; + /* Coriolis force */ + particle[j].fx += 2.0*params.omega_sphere*cos(particle[j].yc)*particle[j].vy; + particle[j].fy += 2.0*params.omega_sphere*sincos*particle[j].vx; + } + + /* add potential force on sphere */ + if ((SPHERE)&&(ADD_POTENTIAL_SPHERE)) + { + i0 = (int)(particle[j].xc*dphi_inv); + j0 = (int)(particle[j].yc*dtheta_inv); +// printf("(fx, fy) = (%.3lg, %.3lg)\n", POT_SPHERE_STRENGTH*wsphere[i0+NY_SPHERE+j0].nablax, POT_SPHERE_STRENGTH*wsphere[i0+NY_SPHERE+j0].nablay); + particle[j].fx -= POT_SPHERE_STRENGTH*wsphere[i0*NY_SPHERE+j0].nablax; + particle[j].fy -= POT_SPHERE_STRENGTH*wsphere[i0*NY_SPHERE+j0].nablay; + } + if (FLOOR_FORCE) { if (particle[j].fx > FMAX) particle[j].fx = FMAX; @@ -2814,13 +2906,15 @@ void animation() if (ADD_ABSORBERS) free(absorber); + if ((ADD_POTENTIAL_SPHERE)&&(SPHERE_POTENTIAL == SPP_WELLS)) free(wells); + if (SAVE_TIME_SERIES) { fclose(lj_time_series); fclose(lj_final_position); } - if (DRAW_SPHERE) free(wsphere); + if (SPHERE) free(wsphere); fclose(lj_log); } diff --git a/percolation.c b/percolation.c index 47060ca..e8352d2 100644 --- a/percolation.c +++ b/percolation.c @@ -111,6 +111,7 @@ #define SCALE 0 /* set to 1 to adjust color scheme to variance of field */ #define SLOPE 1.0 /* sensitivity of color on wave amplitude */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ +#define COLOR_RANGE 1.0 /* max range of color (default: 1.0) */ #define HUE_CLOSED 350.0 /* color hue of closed cells */ #define HUE_OPEN 200.0 /* color hue of open (dry) cells */ diff --git a/rde.c b/rde.c index d475ea0..27f5c8e 100644 --- a/rde.c +++ b/rde.c @@ -40,7 +40,7 @@ #include #define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -48,8 +48,8 @@ #define WINWIDTH 1920 /* window width */ #define WINHEIGHT 1150 /* window height */ -#define NX 870 /* number of grid points on x axis */ -#define NY 520 /* number of grid points on y axis */ +#define NX 1600 /* number of grid points on x axis */ +#define NY 960 /* number of grid points on y axis */ #define HRES 1 /* factor for high resolution plots */ #define XMIN -2.0 @@ -59,7 +59,7 @@ /* Choice of simulated equation */ -#define RDE_EQUATION 11 /* choice of reaction term, see list in global_3d.c */ +#define RDE_EQUATION 121 /* choice of reaction term, see list in global_3d.c */ #define NFIELDS 3 /* number of fields in reaction-diffusion equation */ #define NLAPLACIANS 0 /* number of fields for which to compute Laplacian */ @@ -141,7 +141,7 @@ /* Physical parameters of wave equation */ -#define DT 0.000002 +#define DT 0.000007 #define VISCOSITY 0.02 #define POISSON_STIFFNESS 1.0 /* stiffness of Poisson equation solver for incompressible Euler */ @@ -161,6 +161,8 @@ #define K_KURAMOTO 20.0 /* constant in Kuramoto model */ #define NU_KURAMOTO 0.0 /* viscosity in Kuramoto model */ #define OMEGA_KURAMOTO 0.02 /* frequency in Kuramoto model */ + +/* Keller-Segel model */ #define A_KS 0.15 /* coupling constant in Keller-Segel model */ #define C_KS 3.5 /* coupling constant in Keller-Segel model */ #define D_KSU 1.0 /* organism diffusion in Keller-Segel model */ @@ -168,6 +170,17 @@ #define KS_RSCALE 0.07 /* scaling factor for reaction term of Keller-Segel model */ #define KS_UMAX 100.0 /* max value for u component of Keller-Segel model */ #define KS_CHIMAX_FACTOR 30.0 /* limiter on the value of chi/D_KSU */ + +/* Gray-Scott model */ +#define A_GS 0.029 /* coupling constant in Gray-Scott model */ +#define B_GS 0.06 /* coupling constant in Gray-Scott model */ +#define D_GSU 1.0 /* u diffusion in Gray-Scott model */ +#define D_GSV 2.0 /* v diffusion in Gray-Scott model */ +#define GS_RSCALE 0.1 /* scaling factor for reaction term of Gray-Scott model */ +#define GS_DSCALE 0.1 /* scaling factor for diffusion term of Gray-Scott model */ + +#define PDISC_DISTANCE 0.5 /* constant for Poisson disc initial condition */ + #define V_MAZE 0.4 /* potential in walls of maze */ #define BZQ 0.0008 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */ @@ -217,6 +230,16 @@ #define RPSLZB_INITIAL_TIME 0 /* initial time during which rpslzb remains constant */ #define RPSLZB_FINAL_TIME 500 /* final time during which rpslzb remains constant */ +#define CHANGE_A_GS 0 /* set to 1 to change first parameter in Gray-Scott model */ +#define A_GS_CHANGE -0.02 /* amount by which to decrease a_gs parameter */ +#define A_GS_INITIAL_TIME 0 /* initial time during which a_gs remains constant */ +#define A_GS_FINAL_TIME 100 /* final time during which a_gs remains constant */ + +#define CHANGE_B_GS 0 /* set to 1 to change second parameter in Gray-Scott model */ +#define B_GS_CHANGE -0.01 /* amount by which to decrease b_gs parameter */ +#define B_GS_INITIAL_TIME 0 /* initial time during which b_gs remains constant */ +#define B_GS_FINAL_TIME 100 /* final time during which b_gs remains constant */ + #define CHANGE_FLOW_SPEED 0 /* set to 1 to change speed of laminar flow */ #define IN_OUT_FLOW_BC 0 /* type of in-flow/out-flow boundary conditions for Euler equation, 0 for no b.c. */ #define IN_OUT_BC_FACTOR 0.001 /* factor of convex combination between old and new flow */ @@ -251,7 +274,7 @@ /* Parameters for length and speed of simulation */ #define NSTEPS 1400 /* number of frames of movie */ -#define NVID 12 /* number of iterations between images displayed on screen */ +#define NVID 30 /* number of iterations between images displayed on screen */ #define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */ #define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation */ #define MAX_DT 0.024 /* maximal value of integration step */ @@ -273,6 +296,8 @@ #define PLOT_SPHERE 1 /* draws fields on a sphere */ #define SIMULATION_GRID 1 /* type of simulation grid, for certain equations */ #define OPTIMIZE_GRID 1 /* set to 1 to optimize grid by evolving grid points */ +#define OPTIMIZE_STEPS 200 /* number of grid optimization steps */ +#define WEIGHT_GRID 1 /* set to 1 to compute correcting weights for Laplacian */ #define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ #define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ @@ -287,12 +312,12 @@ /* Plot type - color scheme */ #define CPLOT 0 -#define CPLOT_B 0 +#define CPLOT_B 80 /* Plot type - height of 3D plot */ -#define ZPLOT 0 /* z coordinate in 3D plot */ -#define ZPLOT_B 0 /* z coordinate in second 3D plot */ +#define ZPLOT 80 /* z coordinate in 3D plot */ +#define ZPLOT_B 80 /* z coordinate in second 3D plot */ #define AMPLITUDE_HIGH_RES 1 /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */ #define NON_DIRICHLET_BC 0 /* set to 1 to draw only facets in domain, if field is not zero on boundary */ @@ -308,12 +333,15 @@ #define FLOODING 0 /* set to 1 for drawing water when higher than continents */ #define FLOODING_VSHIFT -10.0 /* controls when wave is considered higher than land */ #define FLOODING_VSHIFT_2D 0.56 /* controls when wave is considered higher than land */ +#define TIME_SLICE_THRESHOLD 0.1 /* threshold of field 0 for Z_TIME_SLICE */ #define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */ #define PRINT_TIME 0 /* set to 1 to print running time */ #define PRINT_VISCOSITY 0 /* set to 1 to print viscosity */ #define PRINT_RPSLZB 0 /* set to 1 to print rpslzb parameter */ +#define PRINT_A_GS 0 /* set to 1 to print a_gs parameter */ +#define PRINT_B_GS 0 /* set to 1 to print b_gs parameter */ #define PRINT_PROBABILITIES 0 /* set to 1 to print probabilities (for Ehrenfest urn configuration) */ #define PRINT_NOISE 0 /* set to 1 to print noise intensity */ #define PRINT_FLOW_SPEED 0 /* set to 1 to print speed of flow */ @@ -337,8 +365,8 @@ /* Color schemes, see list in global_pdes.c */ -#define COLOR_PALETTE 15 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 17 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 14 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ #define COLOR_OUT_R 1.0 /* color outside domain */ @@ -352,8 +380,9 @@ #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 VSHIFT_AMPLITUDE -0.75 /* additional shift for wave amplitude */ -#define VSCALE_AMPLITUDE 1.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSHIFT_AMPLITUDE -0.4 /* additional shift for wave amplitude */ +#define VSCALE_AMPLITUDE 2.5 /* additional scaling factor for color scheme P_3D_AMPLITUDE */ +#define VSCALE_TIMESLICE 1.6 /* additional scaling factor for color scheme Z_TIMESLICE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ #define CURL_SCALE 1.25 /* scaling factor for curl representation */ #define RESCALE_COLOR_IN_CENTER 0 /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */ @@ -395,8 +424,8 @@ #define MAZE_XSHIFT 0.0 /* horizontal shift of maze */ #define MAZE_WIDTH 0.04 /* half width of maze walls */ -#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */ +#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ +#define COLORBAR_RANGE 1.5 /* scale of color scheme bar */ #define COLORBAR_RANGE_B 2.5 /* scale of color scheme bar for 2nd part */ #define ROTATE_COLOR_SCHEME 0 /* set to 1 to draw color scheme horizontally */ #define CIRC_COLORBAR 0 /* set to 1 to draw circular color scheme */ @@ -478,16 +507,16 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO #define VENUS_NODATA_FACTOR 0.5 /* altitude to assign to DEM points without data (fraction of mean altitude) */ #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 XY_SCALING_FACTOR 1.7 /* overall scaling factor for on-screen (x,y) coordinates after projection */ #define ZMAX_FACTOR 1.0 /* max value of z coordinate for REP_PROJ_3D representation */ -#define XSHIFT_3D -0.0 /* overall x shift for REP_PROJ_3D representation */ -#define YSHIFT_3D -0.0 /* overall y shift for REP_PROJ_3D representation */ +#define XSHIFT_3D 0.0 /* overall x shift for REP_PROJ_3D representation */ +#define YSHIFT_3D 0.0 /* overall y shift for REP_PROJ_3D representation */ #define BORDER_PADDING 0 /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */ -#define DRAW_ARROW 1 /* set to 1 to draw arrow above sphere */ +#define DRAW_ARROW 0 /* set to 1 to draw arrow above sphere */ #define ZMAX 3.0 /* maximal value of z coordinate */ #define ZMIN -3.0 /* maximal value of z coordinate */ -#define RSCALE 0.02 /* scaling factor of radial component */ +#define RSCALE 0.25 /* scaling factor of radial component */ #define RSCALE_B 1.00 /* experimental, additional radial scaling factor in ij_to_sphere */ #define RSHIFT 0.0 /* shift in radial component */ #define RMAX 4.0 /* max value of radial component */ @@ -496,13 +525,12 @@ int reset_view = 0; /* switch to reset 3D view parameters (for option RO /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ -#define VMAX 100.0 /* max value of wave amplitude */ +#define VMAX 10.0 /* max value of wave amplitude */ #define TEST_GRADIENT 0 /* print norm squared of gradient */ - #define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT) /* to save computing time, to be improved */ #define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_ANGLE))) -#define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)||(PRINT_AVERAGE_SPEED)) +#define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_A_GS)||(PRINT_B_GS)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)||(PRINT_AVERAGE_SPEED)) #define COMPUTE_PRESSURE ((ZPLOT == Z_EULER_PRESSURE)||(CPLOT == Z_EULER_PRESSURE)||(ZPLOT_B == Z_EULER_PRESSURE)||(CPLOT_B == Z_EULER_PRESSURE)) #define ASYM_SPEED_COLOR (VMEAN_SPEED == 0.0) @@ -1301,6 +1329,8 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) first = 0; } +// printf("[evolve_wave_half], equation = %i\n", RDE_EQUATION); + /* smooth vector fields at poles */ // if ((SPHERE)&&(!SMOOTHBLOCKS)) for (k=0; k 1280) { boxheight = 0.035; - boxwidth = 0.21; + boxwidth = 0.15; +// boxwidth = 0.21; if (left) { xbox = XMIN + 0.4; - xtext = XMIN + 0.2; +// xtext = XMIN + 0.2; + xtext = XMIN + 0.3; } else { @@ -2159,6 +2235,7 @@ void print_parameters(double *phi[NFIELDS], t_rde rde[NX*NY], short int xy_in[NX if (PRINT_TIME) sprintf(message, "Time %.3f", time); else if (PRINT_VISCOSITY) sprintf(message, "Viscosity %.3f", viscosity); else if (PRINT_RPSLZB) sprintf(message, "b = %.3f", rpslzb); + else if (PRINT_A_GS) sprintf(message, "a = %.3f", a_gs); else if (PRINT_NOISE) sprintf(message, "noise %.3f", noise); else if (PRINT_FLOW_SPEED) sprintf(message, "Speed %.3f", flow_speed); else if (PRINT_AVERAGE_SPEED) @@ -2181,6 +2258,22 @@ void print_parameters(double *phi[NFIELDS], t_rde rde[NX*NY], short int xy_in[NX write_text(pos[0], pos[1], message2); } } + if (PRINT_B_GS) + { + y -= 0.1; + erase_area_hsl(xbox, y + 0.02, boxwidth, boxheight, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + sprintf(message, "b = %.3f", b_gs); + if (PLOT_3D) + { + write_text(xtext, y, message); + } + else + { + xy_to_pos(xtext, y, pos); + write_text(pos[0], pos[1], message); + } + } } } @@ -2247,6 +2340,32 @@ double rpslzb_schedule(int i) } } +double a_gs_schedule(int i) +{ + double ratio; + + if (i < A_GS_INITIAL_TIME) return (A_GS); + else if (i > NSTEPS - A_GS_FINAL_TIME) return(A_GS - A_GS_CHANGE); + else + { + ratio = (double)(i - A_GS_INITIAL_TIME)/(double)(NSTEPS - A_GS_INITIAL_TIME - A_GS_FINAL_TIME); + return (A_GS - ratio*A_GS_CHANGE); + } +} + +double b_gs_schedule(int i) +{ + double ratio; + + if (i < B_GS_INITIAL_TIME) return (B_GS); + else if (i > NSTEPS - B_GS_FINAL_TIME) return(B_GS - B_GS_CHANGE); + else + { + ratio = (double)(i - B_GS_INITIAL_TIME)/(double)(NSTEPS - B_GS_INITIAL_TIME - B_GS_FINAL_TIME); + return (B_GS - ratio*B_GS_CHANGE); + } +} + double flow_speed_schedule(int i) { double ratio; @@ -2395,7 +2514,7 @@ void animation() // if (ADD_TRACERS) tracers = (double *)malloc(2*NSTEPS*N_TRACERS*sizeof(double)); if (ADD_TRACERS) tracers = (double *)malloc(4*NSTEPS*N_TRACERS*sizeof(double)); - if ((RDE_EQUATION == E_KURAMOTO_SPHERE)||(RDE_EQUATION == E_KELLER_SEGEL_SPHERE)) + if ((RDE_EQUATION == E_KURAMOTO_SPHERE)||(RDE_EQUATION == E_KELLER_SEGEL_SPHERE)||(RDE_EQUATION == E_GRAY_SCOTT_SPHERE)) ngridpoints = initialize_simulation_grid_sphere(wsphere); dx = (XMAX-XMIN)/((double)NX); @@ -2474,8 +2593,27 @@ void animation() // init_linear_blob_sphere(0, 1.3*PI, 0.65*PI, 0.0, 0.0, 0.3, 0.1, 0.1, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); - init_random(0.5, 0.4, phi, xy_in, wsphere); -// init_onefield_random(1, 0.2, 0.1, phi, xy_in, wsphere); +// init_random(0.5, 0.1, phi, xy_in, wsphere); + +// init_random_zerosum(0.4, 0.4, phi, xy_in, wsphere); + +// init_gs_blob(0.0, 0.0, 0.0, 1.0, 0.25, phi, xy_in, wsphere); +// init_gs_blob(0.0, 0.0, 0.0, 1.0, 0.05, phi, xy_in, wsphere); + +// init_gs_blob(0.0, 0.0, 0.0, 1.0, 0.25, phi, xy_in, wsphere); +// init_gs_blobs_hex(6, 4, 1.0, 0.01, 0.0025, phi, xy_in, wsphere); +// init_gs_blobs_square(1, 1, 1.0, 0.01, 0.002, phi, xy_in, wsphere); + +// init_gs_blobs_poisson(0.5, 100, 1.0, 0.007, 0.0015, phi, xy_in, wsphere); + +// init_gs_blobs_sphere(1.0, 0.1, 0.1, phi, xy_in, wsphere, C_SPH_CUBE_B); + init_gs_blobs_sphere(1.0, 0.05, 0.05, phi, xy_in, wsphere, C_SPH_POISSON); + + // init_gs_rect_sphere(1.0, PID, PI, PID-0.02, PID+0.02, phi, xy_in, wsphere, 0); + + +// init_onefield_random(0, 0.5, 0.5, phi, xy_in, wsphere); +// init_onefield_random(1, 0.5, 0.5, phi, xy_in, wsphere); // init_random_smoothed(0.0, 0.4, phi, xy_in, wsphere); // init_expanding_blob_sphere(0, (279.0/180.0)*PI, (115.0/180.0)*PI, 0.75, 0.04, 0.06, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); @@ -2503,7 +2641,7 @@ void animation() if (ADD_MOON_FORCING) compute_forcing_schedule(0, wsphere); if (RDE_EQUATION == E_KURAMOTO_SPHERE) - convert_fields_from_grids(phi[1], phi[0], wsphere); + convert_fields_from_grids(phi[1], phi[0], rde, wsphere); blank(); glColor3f(0.0, 0.0, 0.0); @@ -2513,7 +2651,7 @@ void animation() printf("Drawing wave\n"); - draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); + draw_wave_rde(0, 0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, VISCOSITY, noise); if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0); @@ -2542,6 +2680,8 @@ void animation() } } if (CHANGE_RPSLZB) rpslzb = rpslzb_schedule(i); + if (CHANGE_A_GS) a_gs = a_gs_schedule(i); + if (CHANGE_B_GS) b_gs = b_gs_schedule(i); if (CHANGE_FLOW_SPEED) flow_speed = flow_speed_schedule(i); else flow_speed = IN_OUT_FLOW_AMP; if (ADD_MOON_FORCING) compute_forcing_schedule(i, wsphere); @@ -2553,7 +2693,7 @@ void animation() } printf("Drawing wave %i\n", i); - draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); + draw_wave_rde(i, 0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); // nvid = (int)((double)NVID*(1.0 + (ACCELERATION_FACTOR - 1.0)*(double)i/(double)NSTEPS)); /* increase integration step */ @@ -2652,7 +2792,7 @@ void animation() if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE)) { - draw_wave_rde(1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B); + draw_wave_rde(i, 1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B); if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, i, 0, 1.0); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise); @@ -2684,7 +2824,7 @@ void animation() { if (DOUBLE_MOVIE) { - draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); + draw_wave_rde(NSTEPS, 0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1); if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, NSTEPS, 0, 1.0); // draw_billiard(); if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise); @@ -2696,7 +2836,7 @@ void animation() else for (i=0; i HASHX-1) hashx_sphere[j] = HASHX-1; dx_sphere[j] = DPI/(double)hashx_sphere[j]; - printf("hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]); +// printf("hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]); fprintf(lj_log, "hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]); } for (j=HASHY/2 + 1; j 0) + { + erase_area_hsl(xmid, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0); + glColor3f(1.0, 1.0, 1.0); + if (SPHERE_GRAVITY == 1) + sprintf(message, "g angle %.2f", params.g_angle*180.0/PI - 90.0); + else sprintf(message, "g angle %.2f", params.g_angle*180.0/PI); + write_text(xmidtext, y, message); + } + if (CHANGE_RADIUS) { erase_area_hsl(xmid, y + 0.025*scale, 0.3*scale, 0.05*scale, 0.0, 0.9, 0.0); @@ -11685,7 +11771,7 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) /* initialize all particles, obstacles, and the hashgrid */ { int i, j, k, n, type, nactive = 0, hashcell; - double x, y, h, xx, yy, rnd, angle, dx, dy; + double x, y, h, xx, yy, rnd, angle, dx, dy, mratio; printf("Initializing configuration\n"); @@ -11939,13 +12025,57 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) // } /* position-dependent particle type */ - if (POSITION_DEPENDENT_TYPE) for (i=0; i POSITION_DEP_Y) + { + particle[i].type = 3; + particle[i].mass_inv = 1.0/PARTICLE_MASS_C; + particle[i].radius = MU_C; + } + else if (particle[i].yc > YMAX - POSITION_DEP_Y) + { + particle[i].type = 2; + particle[i].mass_inv = 1.0/PARTICLE_MASS_B; + particle[i].radius = MU_B; + } + } + else + { + /* TODO */ + } + } + break; + } + case (PDIC_MASSY): + { + for (i=0; i DPI) absorber[k].xc -= DPI; } @@ -17832,12 +17993,31 @@ int init_absorbers(t_absorber absorber[NMAX_ABSORBERS]) for (k=0; k DPI) absorber[k].xc -= DPI; } break; } + case (AP_PDISC): + { + point = (t_point *)malloc(NMAX_ABSORBERS*sizeof(t_point)); + ndiscs = generate_poisson_discs(point, NMAX_ABSORBERS, PI+0.1, DPI-0.1, 0.1, PI-0.1, ABSORBER_PDIST, 1); + + for (k=0; k PI)&&(point[k].xc < DPI)) + { + absorber[nabsorb].xc = point[k].xc; + absorber[nabsorb].yc = point[k].yc; + absorber[nabsorb].radius = radius; + nabsorb++; + } + } + + free(point); + break; + } } return(nabsorb); } @@ -17871,7 +18051,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_absorber absorber[NMAX_ABSORBERS]) + 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 */ { printf("Drawing frame\n"); @@ -17879,8 +18059,9 @@ 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)) { - compute_background_color(particle, segment, obstacle, bg_color, hashgrid); - draw_background(particle, segment, obstacle, bg_color, hashgrid); + compute_background_color(particle, segment, obstacle, bg_color, hashgrid, wsphere); + if (bg_color == BG_POTENTIAL) draw_potential(wsphere); + else draw_background(particle, segment, obstacle, bg_color, hashgrid); } // else if (!TRACER_PARTICLE) blank(); @@ -17934,7 +18115,7 @@ void draw_frame_3d(int i, int plot, int bg_color, int ncollisions, int traj_posi /* initialize background colors */ if ((COLOR_BACKGROUND)&&(bg_color > 0)) - compute_background_color(particle, segment, obstacle, bg_color, hashgrid); + compute_background_color(particle, segment, obstacle, bg_color, hashgrid, wsphere); /* compute all particle colors */ compute_all_particle_colors(particle, cluster, plot); @@ -17970,5 +18151,5 @@ void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_positio 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, absorber); + 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); } diff --git a/sub_lj_sphere.c b/sub_lj_sphere.c index 7e529b5..4bb72ec 100644 --- a/sub_lj_sphere.c +++ b/sub_lj_sphere.c @@ -690,11 +690,11 @@ void compute_all_particle_colors(t_particle particle[NMAXCIRCLES], t_cluster clu } -void compute_background_color(t_particle particle[NMAXCIRCLES], t_segment segment[NMAXSEGMENTS], 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], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE]) /* color background according to particle properties */ { - int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs; - double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue, valx, valy, lum; + int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs, i0, j0, cell, i1, j1, cell1; + double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue, valx, valy, lum, nablax, nablay, norm, pscal, ca; static int first = 1, counter = 0; static double area_factor; @@ -1060,6 +1060,22 @@ void compute_background_color(t_particle particle[NMAXCIRCLES], t_segment segmen hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIRECTION); break; } + case (BG_POTENTIAL): + { + i0 = (int)(0.5*(hashgrid[n].x1+hashgrid[n].x2)*(double)NX_SPHERE/DPI); + j0 = (int)(0.5*(hashgrid[n].y1+hashgrid[n].y2)*(double)NY_SPHERE/PI); + cell = i0*NY_SPHERE+j0; + if (cell < 0) cell = 0; + + value = atan(BG_POTENTIAL_SLOPE*wsphere[cell].potential)/PI + 0.5; + + hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*value; + hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb); + + if (SHADE_BG_COLOR_2D) + for (k=0; k<3; k++) rgb[k] *= wsphere[cell].cos_angle_pot; + break; + } } if (DOUBLE_MOVIE) @@ -1286,6 +1302,7 @@ void init_lj_sphere(t_lj_sphere *wsphere) wsphere[i*NY_SPHERE+j].z = -wsphere[i*NY_SPHERE+j].ctheta; wsphere[i*NY_SPHERE+j].radius = 1.0; + wsphere[i*NY_SPHERE+j].initial_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]; @@ -1361,7 +1378,7 @@ 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]) +void add_particle_to_sphere(int i, int j, int i0, int j0, 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, k; @@ -1554,7 +1571,8 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC if (dist < r) { - h = 1.0 + 1.5*sqrt(r*r - dist*dist); +// h = wsphere[i*NY_SPHERE+j].initial_radius + 1.5*sqrt(r*r - dist*dist); + h = wsphere[i0*NY_SPHERE+j0].initial_radius + 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]; @@ -1564,6 +1582,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC } } } + } 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) @@ -1758,7 +1777,7 @@ void draw_segments_sphere(t_segment segment[NMAXSEGMENTS], t_lj_sphere wsphere[N /* 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; + double x1, y1, x2, y2, x, y, dt, length, ca; static double dphi, dtheta; static int first = 1; @@ -1796,11 +1815,19 @@ void draw_segments_sphere(t_segment segment[NMAXSEGMENTS], t_lj_sphere wsphere[N 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; + if (wsphere[cell].locked == 0) + { + ca = wsphere[cell].cos_angle_sphere; + ca = (ca + 1.0)*0.4 + 0.2; + + wsphere[cell].r = ca; + wsphere[cell].g = ca; + wsphere[cell].b = ca; + + wsphere[cell].radius += 0.05; + wsphere[cell].locked = 1; + } } } } @@ -1840,11 +1867,131 @@ void draw_absorbers_sphere(t_absorber absorber[NMAX_ABSORBERS], t_lj_sphere wsph { cell = i*NY_SPHERE+j; wsphere[cell].locked = 1; - wsphere[cell].radius = 0.95; +// wsphere[cell].radius = 0.95; + wsphere[cell].radius = wsphere[cell].initial_radius - 0.05; } } } +void compute_gradient_potential_sphere(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE]) +/* compute the gradient of the potential on the sphere */ +{ + int i, j, cell; + double dphi_inv, dtheta_inv, stheta_inv; + + dphi_inv = 0.5*(double)NX_SPHERE/DPI; + dtheta_inv = 0.5*(double)NY_SPHERE/PI; + + /* bulk */ + for (j=1; j max) max = x; + if (y > max) max = y; + if (z > max) max = z; + wsphere[i].potential = 1.0/max - 1.0; + } + break; + } + } + + for (i=0; i 0)) { @@ -1885,8 +2032,8 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA if (TRACER_PARTICLE) 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_FIXED_SEGMENTS) +// draw_segments_sphere(segment, wsphere); // if (ADD_ABSORBERS) // draw_absorbers_sphere(absorber, wsphere); @@ -1910,13 +2057,14 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA if (particle[part].active) { deltaj = (int)(2.5*particle[part].radius/dtheta); + i0 = 0; j0 = (int)(particle[part].yc/dtheta); jmax = j0 + deltaj; #pragma omp parallel for private(i,j,x,y,dist,r) for (j=0; j 0)&&(npoints < nmaxpoints)) + { + /* randomly select an active circle */ + i = rand()%(npoints); + while (!active_poisson[i]) i = rand()%(npoints); + /* generate new candidates */ + naccepted = 0; + for (j=0; j= dpoisson); + if (distance < dist1) dist1 = distance; + } + } + if ((far)&&(npoints < nmaxpoints)) /* accept new point */ + { + printf("New grid point at (%.3f,%.3f) accepted\n", x, y); + printf("Minimal distance = %.3lg\n", dist1); + printf("%i grid points generated\n", npoints); + px[npoints] = x; + py[npoints] = y; + active_poisson[npoints] = 1; + npoints++; + n_p_active++; + naccepted++; + } + } + if (naccepted == 0) /* inactivate circle i */ + { + active_poisson[i] = 0; + n_p_active--; + } + + printf("%i active points\n", n_p_active); + } + + printf("Generated %i points\n", npoints); + return(npoints); +} + + + + +void init_gs_blobs_poisson(double dpoisson, double nmax, double amplitude, double scalex, double scaley, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY]) +/* initialise Gray-Scott model with Poisson disc sample of ellipses */ +{ + int i, j, k, npoints, p, q; + double x, y, xy[2], distmin, dist2, module, scalex2, scaley2, dx, dy; + double *px, *py; + + px = (double *)malloc(nmax*sizeof(double)); + py = (double *)malloc(nmax*sizeof(double)); + + scalex2 = scalex*scalex; + scaley2 = scaley*scaley; + printf("Initialising field\n"); + + npoints = generate_poisson_sample(XMIN, XMAX, YMIN, YMAX, dpoisson, nmax, px, py); + + for (i=0; i phimin)&&(phix < phimax)&&(thetay > thetamin)&&(thetay < thetamax)) + { + phi[0][i] = 0.0; + phi[1][i] = amplitude; + phi[2][i] = 1.0 - 0.5*amplitude; + } + } + else for (i=0; i phimin)&&(phix < phimax)&&(thetay > thetamin)&&(thetay < thetamax)) module = 1.0; + else module = 0.0; + { + phi[0][i] = 0.0; + phi[1][i] = module*amplitude; + phi[2][i] = 1.0 - 0.5*module*amplitude; + } + } +} + + double distance_to_segment(double x, double y, double x1, double y1, double x2, double y2) /* distance of (x,y) to segment from (x1,y1) to (x2,y2) */ { @@ -4235,6 +4603,27 @@ void adjust_height(double *phi[NFIELDS], t_rde rde[NX*NY]) } } +void update_time_slice(int time, double *phi[NFIELDS], t_rde rde[NX*NY]) +/* adjust height of field */ +{ + int i, j; + static int first = 1; + + if (first) + { + #pragma omp parallel for private(i) + for (i=0; i wmax) wmax = wsphere[i].laplace_weight; + } + + ratio = 1.0/wmax; + printf("Ratio = %.3lg\n", ratio); + + for (i=0; i 1.5) wsphere[i].laplace_weight = 1.5; +// if (wsphere[i].laplace_weight < 0.5) wsphere[i].laplace_weight = 0.5; + + printf("Grid point %i weight = %.3lg\n", i, wsphere[i].laplace_weight); + } +} + int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY]) /* initialize Poisson random simulation grid on a sphere */ { @@ -9526,7 +9978,8 @@ int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY]) case (GRID_CUBIC): { npoints = generate_cubic_grid_sphere(wsphere); - if (OPTIMIZE_GRID) optimize_cubic_grid_sphere(wsphere, npoints, 600); + if (OPTIMIZE_GRID) optimize_cubic_grid_sphere(wsphere, npoints, OPTIMIZE_STEPS); + if (WEIGHT_GRID) compute_cubic_grid_sphere_weights(wsphere, npoints); break; } } @@ -9541,7 +9994,7 @@ int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY]) return(npoints); } -void convert_fields_from_grids(double phi_in[NX*NY], double phi_out[NX*NY], t_wave_sphere wsphere[NX*NY]) +void convert_fields_from_grids(double phi_in[NX*NY], double phi_out[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) /* convert field phi_in in grid coordinates to field phi_out in latitude-longitude coordinates */ { // int i; diff --git a/sub_wave.c b/sub_wave.c index 3b42efb..a958a6e 100644 --- a/sub_wave.c +++ b/sub_wave.c @@ -353,21 +353,21 @@ double gaussian() -// 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 radious r, turned by apoly Pi/2 */ -// { -// int condition = 1, k; -// double omega, cw, angle; -// -// omega = DPI/((double)npoly); -// cw = cos(omega*0.5); -// for (k=0; k 0)&&(npoints < nmaxpoints)) + { + /* randomly select an active circle */ + i = rand()%(npoints); + while (!active_poisson[i]) i = rand()%(npoints); + /* generate new candidates */ + naccepted = 0; + for (j=0; j DPI) x -= DPI; +// if (x < 0.0) x += DPI; +// +// far = 1; +// if (y > PI - polar_padding) far = 0; +// if (y < polar_padding) far = 0; + + far = 1; + dist1 = 10.0; + for (k=0; k= dpoisson); + if (distance < dist1) dist1 = distance; + } + if ((far)&&(npoints < nmaxpoints)) /* accept new point */ + { + printf("New grid point at (%.3f,%.3f) accepted\n", x, y); + printf("Minimal distance = %.3lg\n", dist1); + printf("%i grid points generated\n", npoints); + px[npoints] = x; + py[npoints] = y; + active_poisson[npoints] = 1; + npoints++; + n_p_active++; + naccepted++; + } + } + if (naccepted == 0) /* inactivate circle i */ + { + active_poisson[i] = 0; + n_p_active--; + } + + printf("%i active points\n", n_p_active); + } + + printf("Generated %i points\n", npoints); + return(npoints); +} + + +int init_circle_sphere(t_circles_sphere *circles, int circle_pattern) +/* initialise circles on sphere */ +/* for billiard shape D_SPHERE_CIRCLES */ +{ + int ncircles, k; + double alpha, beta, gamma; + double px[NMAXCIRC_SPHERE], py[NMAXCIRC_SPHERE]; + + switch (circle_pattern) { + case (C_SPH_DODECA): + { + alpha = acos(sqrt(5.0)/3.0); + beta = acos(1.0/3.0); + gamma = asin(sqrt(3.0/8.0)); + + circles[0].theta = 0.0; + circles[0].phi = 0.0; + + for (k=0; k<3; k++) + { + circles[1+k].theta = alpha; + circles[1+k].phi = (double)k*DPI/3.0; + } + + for (k=0; k<3; k++) + { + circles[4+k].theta = beta; + circles[4+k].phi = (double)k*DPI/3.0 + gamma; + circles[7+k].theta = beta; + circles[7+k].phi = (double)k*DPI/3.0 - gamma; + } + + for (k=0; k<3; k++) + { + circles[10+k].theta = PI - beta; + circles[10+k].phi = (double)k*DPI/3.0 + PI/3.0 + gamma; + circles[13+k].theta = PI - beta; + circles[13+k].phi = (double)k*DPI/3.0 + PI/3.0 - gamma; + } + + for (k=0; k<3; k++) + { + circles[16+k].theta = PI - alpha; + circles[16+k].phi = (double)k*DPI/3.0 + PI/3.0; + } + + circles[19].theta = PI; + circles[19].phi = 0.0; + + for (k=0; k<20; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 20; + break; + } + case (C_SPH_ICOSA): + { + alpha = acos(1.0/sqrt(5.0)); + + circles[0].theta = 0.0; + circles[0].phi = 0.0; + + for (k=0; k<5; k++) + { + circles[1+k].theta = alpha; + circles[1+k].phi = (double)k*DPI/5.0; + } + + for (k=0; k<5; k++) + { + circles[6+k].theta = PI - alpha; + circles[6+k].phi = (double)k*DPI/5.0 + PI/5.0; + } + + circles[11].theta = PI; + circles[11].phi = 0.0; + + for (k=0; k<12; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 12; + break; + } + case (C_SPH_OCTA): + { + circles[0].theta = 0.0; + circles[0].phi = 0.0; + + for (k=0; k<4; k++) + { + circles[1+k].theta = PID; + circles[1+k].phi = (double)k*PID; + } + + circles[5].theta = PI; + circles[5].phi = 0.0; + + for (k=0; k<6; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 6; + break; + } + case (C_SPH_CUBE): + { + alpha = acos(1.0/3.0); + + circles[0].theta = 0.0; + circles[0].phi = 0.0; + + for (k=0; k<3; k++) + { + circles[1+k].theta = alpha; + circles[1+k].phi = (double)k*DPI/3.0; + circles[4+k].theta = PI - alpha; + circles[4+k].phi = (double)k*DPI/3.0 + PI/3.0; + } + + circles[7].theta = PI; + circles[7].phi = 0.0; + + for (k=0; k<8; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 8; + break; + } + case (C_SPH_CUBE_B): + { + alpha = acos(1.0/sqrt(3.0)); + + for (k=0; k<4; k++) + { + circles[k].theta = alpha; + circles[k].phi = ((double)k + 0.5)*PID; + circles[4+k].theta = PI - alpha; + circles[4+k].phi = ((double)k + 0.5)*PID; + } + + for (k=0; k<8; k++) + { + circles[k].radius = MU; + circles[k].x = sin(circles[k].theta)*cos(circles[k].phi); + circles[k].y = sin(circles[k].theta)*sin(circles[k].phi); + circles[k].z = cos(circles[k].theta); + } + + ncircles = 8; + break; + } + case (C_SPH_POISSON): + { + ncircles = generate_poisson_sample_sphere(PDISC_DISTANCE, NMAXCIRC_SPHERE, px, py); + for (k=0; k #define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 1 /* set to 1 to save memory when writing tiff images */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -72,7 +72,7 @@ /* Choice of the billiard table */ -#define B_DOMAIN 421 /* choice of domain shape, see list in global_pdes.c */ +#define B_DOMAIN 991 /* 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,17 +87,17 @@ #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */ #define PDISC_CONNECT_FACTOR 1.5 /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */ -#define LAMBDA 0.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 LAMBDA 0.6 /* parameter controlling the dimensions of domain */ +#define MU 0.15 /* parameter controlling the dimensions of domain */ +#define MU_B 0.45 /* parameter controlling the dimensions of domain */ +#define NPOLY 5 /* number of sides of polygon */ +#define APOLY 0.0 /* angle by which to turn polygon, in units of Pi/2 */ #define MDEPTH 6 /* depth of computation of Menger gasket */ #define MRATIO 3 /* ratio defining Menger gasket */ #define MANDELLEVEL 1000 /* iteration level for Mandelbrot set */ #define MANDELLIMIT 10.0 /* limit value for approximation of Mandelbrot set */ #define FOCI 1 /* set to 1 to draw focal points of ellipse */ -#define NGRIDX 6 /* number of grid point for grid of disks */ +#define NGRIDX 4 /* number of grid point for grid of disks */ #define NGRIDY 14 /* number of grid point for grid of disks */ #define WALL_WIDTH 0.075 /* width of wall separating lenses */ #define WALL_WIDTH_B 0.1 /* width of wall separating lenses */ @@ -122,7 +122,7 @@ /* Physical parameters of wave equation */ #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_LEFT 0 /* set to 1 to add oscilating boundary condition on the left */ #define OSCILLATE_TOPBOT 0 /* set to 1 to enforce a planar wave on top and bottom boundary */ #define OSCILLATION_SCHEDULE 42 /* oscillation schedule, see list in global_pdes.c */ #define OSCIL_YMAX 0.0375 /* defines oscilling beam range */ @@ -134,7 +134,7 @@ #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.08 /* Courant number in medium B */ +#define COURANT 0.12 /* 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 */ @@ -149,10 +149,10 @@ /* Increasing COURANT speeds up the simulation, but decreases accuracy */ /* 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 6.0 /* period of oscillating source */ +#define ADD_OSCILLATING_SOURCE 1 /* set to 1 to add an oscillating wave source */ +#define OSCILLATING_SOURCE_PERIOD 12.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 N_SOURCES 4 /* 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 5000 /* max time for adding pulses */ @@ -169,10 +169,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2700 /* number of frames of movie */ +#define NSTEPS 3200 /* 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 100 /* time after which to start saving frames */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ #define BOUNDARY_WIDTH 2 /* width of billiard boundary */ #define PRINT_SPEED 0 /* print speed of moving source */ #define PRINT_FREQUENCY 0 /* print frequency (for phased array) */ @@ -199,7 +199,7 @@ /* Color schemes */ -#define COLOR_PALETTE 18 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 15 /* 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 */ @@ -212,8 +212,8 @@ #define PHASE_FACTOR 1.0 /* factor in computation of phase in color scheme P_3D_PHASE */ #define PHASE_SHIFT 0.0 /* shift of phase in color scheme P_3D_PHASE */ #define ATTENUATION 0.0 /* exponential attenuation coefficient of contrast with time */ -#define VSHIFT_AMPLITUDE -0.5 /* additional shift for wave amplitude */ -#define VSCALE_AMPLITUDE 1.0 /* additional scaling factor for wave amplitude */ +#define VSHIFT_AMPLITUDE -0.0 /* additional shift for wave amplitude */ +#define VSCALE_AMPLITUDE 0.5 /* additional scaling factor for wave amplitude */ #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 */ @@ -231,14 +231,14 @@ #define HUEMEAN 180.0 /* mean value of hue for color scheme C_HUE */ #define HUEAMP -180.0 /* amplitude of variation of hue for color scheme C_HUE */ -#define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ +#define DRAW_COLOR_SCHEME 0 /* set to 1 to plot the color scheme */ #define COLORBAR_RANGE 1.7 /* scale of color scheme bar */ #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 1 /* set to 1 to draw a profile of the wave */ +#define DRAW_WAVE_PROFILE 0 /* set to 1 to draw a profile of the wave */ #define HORIZONTAL_WAVE_PROFILE 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 */ @@ -880,11 +880,14 @@ void animation() /* add oscillating waves */ for (source = 0; source < N_SOURCES; source++) { - wave_source_x[source] = polyarc[2*source].xc; - wave_source_y[source] = polyarc[2*source].yc; - source_periods[source] = OSCILLATING_SOURCE_PERIOD/(1.0 + 0.5*(double)source); + wave_source_x[source] = polyarc[source].xc; + wave_source_y[source] = polyarc[source].yc; source_amp[source] = INITIAL_AMP; } + source_periods[0] = OSCILLATING_SOURCE_PERIOD; + source_periods[1] = OSCILLATING_SOURCE_PERIOD*3/4; + source_periods[2] = OSCILLATING_SOURCE_PERIOD*2/3; + source_periods[3] = OSCILLATING_SOURCE_PERIOD/2; for (source = 0; source < N_SOURCES; source++) { dperiod = source_periods[source];