From 4d6547bad03d162890ed2855a674c549902f268b Mon Sep 17 00:00:00 2001 From: Nils Berglund <83530463+nilsberglund-orleans@users.noreply.github.com> Date: Sat, 12 Oct 2024 18:19:46 +0200 Subject: [PATCH] Add files via upload --- global_3d.c | 5 +- global_ljones.c | 32 +- global_pdes.c | 5 +- lennardjones.c | 216 +++++++---- rde.c | 238 +++++++++--- sub_lj.c | 952 ++++++++++++++++++++++++++++++++++++++-------- sub_rde.c | 375 +++++++++++++++++- sub_sphere.c | 5 + sub_wave.c | 108 +++++- sub_wave_3d_rde.c | 166 +++++++- wave_billiard.c | 73 ++-- wave_sphere.c | 112 +++--- 12 files changed, 1867 insertions(+), 420 deletions(-) diff --git a/global_3d.c b/global_3d.c index 3bdbba4..7237b39 100644 --- a/global_3d.c +++ b/global_3d.c @@ -51,6 +51,7 @@ #define SH_SPHERE_OCTAHEDRON 6 /* octahedron embedded in sphere */ #define SH_SPHERE_DODECAHEDRON 7 /* dodecahedron embedded in sphere */ #define SH_SPHERE_ICOSAHEDRON 8 /* icosahedron embedded in sphere */ +#define SH_EARTH 10 /* depth of Earth oceans */ /* Type of rotating viewpoint */ @@ -89,7 +90,8 @@ #define PLANET ((B_DOMAIN == D_SPHERE_EARTH)||(B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)||(B_DOMAIN == D_SPHERE_VENUS)||(B_DOMAIN == D_SPHERE_MERCURY)) #define OTHER_PLANET ((B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)||(B_DOMAIN == D_SPHERE_VENUS)||(B_DOMAIN == D_SPHERE_MERCURY)) -#define RDE_PLANET ((ADAPT_STATE_TO_BC)&&((OBSTACLE_GEOMETRY == D_SPHERE_EARTH)||(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)))) +// #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 NMAX_TRACER_PTS 20 /* max number of tracer points recorded per cell */ @@ -168,6 +170,7 @@ typedef struct 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 */ diff --git a/global_ljones.c b/global_ljones.c index 4ed659a..36be704 100644 --- a/global_ljones.c +++ b/global_ljones.c @@ -21,6 +21,7 @@ #define NMAXPARTNERMOLECULES 30 /* max number of partners of a molecule */ #define NMAXPARTINCLUSTER 500 /* max number of particles in cluster */ #define NMAXBELTS 10 /* max number of conveyor belts */ +#define NMAXSHOVELS 50 /* max number of shovels */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ @@ -53,6 +54,10 @@ #define O_HEX 6 /* hexagonal lattice */ #define O_SIDES 7 /* grid along the sides of the simulation rectangle */ #define O_SIDES_B 71 /* finer grid along the sides of the simulation rectangle */ +#define O_SIEVE 8 /* obstacles form a sieve */ +#define O_SIEVE_B 81 /* obstacles form a sieve, v2 */ +#define O_SIEVE_LONG 82 /* obstacles form a long sieve */ +#define O_SIEVE_LONG_B 83 /* obstacles form a long sieve, version with varying spacing */ /* pattern of additional repelling segments */ #define S_RECTANGLE 0 /* segments forming a rectangle */ @@ -95,6 +100,13 @@ #define S_TWO_CONVEYOR_BELTS 31 /* two angled conveyor belts */ #define S_PERIODIC_CONVEYORS 32 /* one wrapping belt, and one short horizontal belt */ #define S_TEST_CONVEYORS 321 /* test */ +#define S_CONVEYOR_SHOVELS 33 /* conveyor belt with shovels */ +#define S_CONVEYOR_MIXED 34 /* multiple conveyor belts with and without shovels */ +#define S_CONVEYOR_SIEVE 35 /* conveyor belts for polygon sieve */ +#define S_CONVEYOR_SIEVE_B 351 /* conveyor belts for polygon sieve, v2 with backward top conveyor */ +#define S_CONVEYOR_SIEVE_LONG 352 /* conveyor belts for long polygon sieve */ +#define S_MASS_SPECTROMETER 36 /* bins for mass spectrometer */ +#define S_WIND_FORCE 361 /* bins for sorting by wind force */ /* particle interaction */ @@ -278,6 +290,7 @@ #define P_CLUSTER_SIZE 20 /* colors depend on size of connected component */ #define P_CLUSTER_SELECTED 21 /* colors show which clusters are slected for growth */ #define P_COLLISION 22 /* colors depend on number of collision/reaction */ +#define P_RADIUS 23 /* colors depend on particle radius */ /* Rotation schedules */ @@ -289,6 +302,11 @@ #define IP_X 0 /* color depends on x coordinate of initial position */ #define IP_Y 1 /* color depends on y coordinate of initial position */ +/* Space dependence of magnetic field */ + +#define BF_CONST 0 /* constant magnetic field */ +#define BF_SQUARE 1 /* magnetic field concentrated in square */ + /* Interaction types for polyatomic molecules */ #define POLY_STAR 0 /* star-shaped graph (central molecule attracts outer ones) */ @@ -368,6 +386,7 @@ typedef struct double fx; /* x component of force on particle */ double fy; /* y component of force on particle */ double torque; /* torque on particle */ + double damping; /* factor in front of damping coefficient */ double dirmean; /* time averaged direction */ int close_to_boundary; /* has value 1 if particle is close to a boundary */ short int thermostat; /* whether particle is coupled to thermostat */ @@ -456,8 +475,14 @@ typedef struct typedef struct { double xc, yc, radius; /* center and radius of circle */ + double xc0, yc0; /* center of oscillation for option RATTLE_OBSTACLES */ short int active; /* circle is active */ double charge; /* charge of obstacle, for EM simulations */ + double omega0, omega; /* speed of rotation */ + double angle; /* angle of obstacle */ + short int oscillate; /* has value 1 if the obstacles oscillates over time */ + int period; /* oscillation period */ + double amplitude, phase; /* amplitude and phase of oscillation */ } t_obstacle; typedef struct @@ -465,6 +490,7 @@ typedef struct double x1, x2, y1, y2; /* extremities of segment */ double xc, yc; /* mid-point of segment */ double nx, ny; /* normal vector */ + double nangle; /* angle of normal vector */ double c; /* constant term in cartesian eq nx*x + ny*y = c */ double length; /* length of segment */ short int concave; /* corner is concave, to add extra repelling force */ @@ -482,6 +508,7 @@ typedef struct short int inactivate; /* set to 1 for segment to become inactive at time SEGMENT_DEACTIVATION_TIME */ short int conveyor; /* set to 1 for segment to exert lateral force */ double conveyor_speed; /* speed of conveyor belt */ + short int align_torque; /* set to 1 for segment to exert aligning torque */ } t_segment; typedef struct @@ -533,6 +560,9 @@ typedef struct double length; /* distance between (x1,x2) and (y1,y2) */ double angle; /* angle of (x1,x2) - (y1,y2) */ double tx, ty; /* coordinates of tangent vector */ + int nshovels; /* number of shovels */ + double shovel_pos[NMAXSHOVELS]; /* position od each shovel */ + int shovel_segment[NMAXSHOVELS]; /* first segment of each shovel */ } t_belt; typedef struct @@ -556,6 +586,6 @@ typedef struct -int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0; +int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0; FILE *lj_log; diff --git a/global_pdes.c b/global_pdes.c index f94121d..3abed13 100644 --- a/global_pdes.c +++ b/global_pdes.c @@ -106,6 +106,9 @@ #define D_IMAGE 77 /* Taken from image file */ #define D_MAGNETRON 78 /* simplified magnetron */ #define D_MAGNETRON_CATHODE 781 /* simplified magnetron with central cathode */ +#define D_TWOCIRCLES 79 /* two circles of different size */ +#define D_POLYCIRCLES 791 /* one large circle and NPOLY small ones */ +#define D_POLYCIRCLES_ANGLED 792 /* variant of D_POLYCIRCLES with angled small cavities */ /* for wave_sphere.c */ @@ -125,7 +128,7 @@ #define NMAXCIRCLES 10000 /* total number of circles/polygons (must be at least NCX*NCY for square grid) */ #define NMAXPOLY 50000 /* maximal number of vertices of polygonal lines (for von Koch et al) */ -#define NMAXSOURCES 10 /* maximal number of sources */ +#define NMAXSOURCES 30 /* maximal number of sources */ #define C_SQUARE 0 /* square grid of circles */ #define C_HEX 1 /* hexagonal/triangular grid of circles */ diff --git a/lennardjones.c b/lennardjones.c index 4a9e907..92792bb 100644 --- a/lennardjones.c +++ b/lennardjones.c @@ -34,10 +34,10 @@ #include #include /* Sam Leffler's libtiff library. */ #include -#include +#include #define MOVIE 0 /* set to 1 to generate movie */ -#define DOUBLE_MOVIE 1 /* set to 1 to produce movies for wave height and energy simultaneously */ +#define DOUBLE_MOVIE 0 /* set to 1 to produce movies for wave height and energy simultaneously */ #define SAVE_MEMORY 1 /* set to 1 to save memory while saving frames */ #define NO_EXTRA_BUFFER_SWAP 1 /* some OS require one less buffer swap when recording images */ @@ -58,47 +58,50 @@ #define YMIN -1.125 #define YMAX 1.125 /* y interval for 9/16 aspect ratio */ -#define INITXMIN -1.8 -#define INITXMAX 1.8 /* x interval for initial condition */ -#define INITYMIN 0.9 -#define INITYMAX 2.2 /* y interval for initial condition */ +#define INITXMIN -2.1 +#define INITXMAX -2.0 /* x interval for initial condition */ +#define INITYMIN -0.7 +#define INITYMAX -0.6 /* y interval for initial condition */ #define 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.5 -#define ADDXMAX 1.5 /* x interval for adding particles */ -#define ADDYMIN -1.0 -#define ADDYMAX 1.0 /* y interval for adding particles */ +#define ADDXMIN -2.1 +#define ADDXMAX -2.0 /* x interval for adding particles */ +#define ADDYMIN -0.6 +#define ADDYMAX -0.5 /* y interval for adding particles */ #define ADDRMIN 4.75 #define ADDRMAX 6.0 /* r interval for adding particles */ -#define BCXMIN -2.0 +#define BCXMIN -2.2 #define BCXMAX 2.0 /* x interval for boundary condition */ #define BCYMIN -1.125 -#define BCYMAX 2.4 /* y interval for boundary condition */ +#define BCYMAX 1.4 /* y interval for boundary condition */ #define OBSXMIN -2.0 #define OBSXMAX 2.0 /* x interval for motion of obstacle */ -#define CIRCLE_PATTERN 8 /* pattern of circles, see list in global_ljones.c */ +#define CIRCLE_PATTERN 1 /* pattern of circles, see list in global_ljones.c */ #define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */ #define CIRCLE_PATTERN_B 0 /* pattern of circles for additional particles */ #define ADD_FIXED_OBSTACLES 0 /* set to 1 do add fixed circular obstacles */ -#define OBSTACLE_PATTERN 71 /* pattern of obstacles, see list in global_ljones.c */ +#define OBSTACLE_PATTERN 82 /* pattern of obstacles, see list in global_ljones.c */ +#define RATTLE_OBSTACLES 1 /* set to 1 to rattle obstacles (for pattern O_SIEVE_B) */ #define ADD_FIXED_SEGMENTS 1 /* set to 1 to add fixed segments as obstacles */ -#define SEGMENT_PATTERN 182 /* pattern of repelling segments, see list in global_ljones.c */ +#define SEGMENT_PATTERN 361 /* pattern of repelling segments, see list in global_ljones.c */ #define ROCKET_SHAPE 3 /* shape of rocket combustion chamber, see list in global_ljones.c */ #define ROCKET_SHAPE_B 3 /* shape of second rocket */ #define NOZZLE_SHAPE 6 /* shape of nozzle, see list in global_ljones.c */ #define NOZZLE_SHAPE_B 6 /* shape of nozzle for second rocket, see list in global_ljones.c */ -#define BELT_SPEED1 20.0 /* speed of first conveyor belt */ -#define BELT_SPEED2 -10.0 /* speed of second conveyor belt */ +#define BELT_SPEED1 10.0 /* speed of first conveyor belt */ +#define BELT_SPEED2 15.0 /* speed of second conveyor belt */ +#define BELT_SPEED3 6.0 /* speed of second conveyor belt */ +#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.5 /* proportion of particles of first type */ @@ -108,22 +111,22 @@ #define CENTER_PY 0 /* set to 1 to center vertical momentum */ #define CENTER_PANGLE 0 /* set to 1 to center angular momentum */ -#define INTERACTION 171 /* particle interaction, 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 6.0 /* angular frequency of spin-spin interaction */ -#define SPIN_INTER_FREQUENCY_B 6.0 /* angular frequency of spin-spin interaction for second particle type */ -#define MOL_ANGLE_FACTOR 6.0 /* rotation angle for P_MOL_ANGLE color scheme */ +#define 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 4.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 2.5 /* minimal distance in Poisson disc process, controls density of particles */ +#define PDISC_DISTANCE 1.0 /* minimal distance in Poisson disc process, controls density of particles */ #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */ #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */ #define LAMBDA 0.2 /* parameter controlling the dimensions of domain */ -#define MU 0.025 /* parameter controlling radius of particles */ +#define MU 0.02 /* parameter controlling radius of particles */ #define MU_B 0.03 /* parameter controlling radius of particles of second type */ -#define NPOLY 6 /* number of sides of polygon */ +#define NPOLY 3 /* number of sides of polygon */ #define APOLY 0.075 /* angle by which to turn polygon, in units of Pi/2 */ #define AWEDGE 0.5 /* opening angle of wedge, in units of Pi/2 */ #define MDEPTH 4 /* depth of computation of Menger gasket */ @@ -131,8 +134,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 36 /* number of grid point for grid of disks */ -#define NGRIDY 36 /* number of grid point for grid of disks */ +#define NGRIDX 1 /* number of grid point for grid of disks */ +#define NGRIDY 1 /* 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 */ @@ -148,10 +151,10 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 2700 /* number of frames of movie */ -#define NVID 200 /* number of iterations between images displayed on screen */ +#define NSTEPS 4800 /* number of frames of movie */ +#define NVID 120 /* number of iterations between images displayed on screen */ #define NSEG 25 /* number of segments of boundary of circles */ -#define INITIAL_TIME 200 /* time after which to start saving frames */ +#define INITIAL_TIME 0 /* time after which to start saving frames */ #define OBSTACLE_INITIAL_TIME 0 /* time after which to start moving obstacle */ #define BOUNDARY_WIDTH 1 /* width of particle boundary */ #define LINK_WIDTH 2 /* width of links between particles */ @@ -170,8 +173,8 @@ /* Plot type, see list in global_ljones.c */ -#define PLOT 11 -#define PLOT_B 14 /* plot type for second movie */ +#define PLOT 23 +#define PLOT_B 13 /* plot type for second movie */ /* Background color depending on particle properties */ @@ -185,7 +188,7 @@ #define DRAW_CLUSTER_LINKS 0 /* set to 1 to draw links between particles in cluster */ #define ALTITUDE_LINES 0 /* set to 1 to add horizontal lines to show altitude */ #define COLOR_SEG_GROUPS 0 /* set to 1 to collor segment groups differently */ -#define N_PARTICLE_COLORS 200 /* number of colors for P_NUMBER color scheme */ +#define N_PARTICLE_COLORS 300 /* number of colors for P_NUMBER color scheme */ #define INITIAL_POS_TYPE 0 /* type of initial position dependence */ #define ERATIO 0.995 /* ratio for time-averaging in P_EMEAN color scheme */ #define DRATIO 0.999 /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */ @@ -240,7 +243,7 @@ #define PARTICLE_HUE_MIN 359.0 /* color of original particle */ #define PARTICLE_HUE_MAX 0.0 /* color of saturated particle */ #define PARTICLE_EMIN 10.0 /* energy of particle with coolest color */ -#define PARTICLE_EMAX 100.0 /* energy of particle with hottest color */ +#define PARTICLE_EMAX 2000.0 /* energy of particle with hottest color */ #define HUE_TYPE0 320.0 /* hue of particles of type 0 */ #define HUE_TYPE1 60.0 /* hue of particles of type 1 */ #define HUE_TYPE2 320.0 /* hue of particles of type 2 */ @@ -251,21 +254,27 @@ #define HUE_TYPE7 150.0 /* hue of particles of type 7 */ #define BG_FORCE_SLOPE 7.5e-8 /* contant in BG_FORCE backgound color scheme*/ -#define RANDOM_RADIUS 0 /* set to 1 for random circle radius */ +#define RANDOM_RADIUS 1 /* set to 1 for random particle radius */ +#define RANDOM_RADIUS_MIN 0.25 /* min of random particle radius (default 0.75) */ +#define RANDOM_RADIUS_RANGE 1.5 /* range of random particle radius (default 0.5) */ +#define ADAPT_MASS_TO_RADIUS 0 /* set to positive value to for mass prop to power of radius */ +#define ADAPT_DAMPING_TO_RADIUS 0.5 /* set to positive value to for friction prop to power of radius */ +#define ADAPT_DAMPING_FACTOR 0.5 /* factor by which damping is adapted to radius */ #define DT_PARTICLE 3.0e-6 /* time step for particle displacement */ #define KREPEL 50.0 /* constant in repelling force between particles */ -#define EQUILIBRIUM_DIST 2.3 /* Lennard-Jones equilibrium distance */ +#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 500.0 /* damping coefficient of particles */ +#define DAMPING 150.0 /* damping coefficient of particles */ #define INITIAL_DAMPING 1000.0 /* damping coefficient of particles during initial phase */ -#define DAMPING_ROT 1000.0 /* damping coefficient for rotation of particles */ +#define DAMPING_ROT 5.0 /* damping coefficient for rotation of particles */ #define PARTICLE_MASS 2.0 /* mass of particle of radius MU */ #define PARTICLE_MASS_B 2.0 /* mass of particle of radius MU_B */ #define PARTICLE_INERTIA_MOMENT 0.5 /* moment of inertia of particle */ #define PARTICLE_INERTIA_MOMENT_B 0.5 /* moment of inertia of second type of particle */ -#define V_INITIAL 50.0 /* initial velocity range */ -#define OMEGA_INITIAL 50.0 /* initial angular velocity range */ +#define V_INITIAL 800.0 /* initial velocity range */ +#define OMEGA_INITIAL 100.0 /* initial angular velocity range */ #define VICSEK_VMIN 1.0 /* minimal speed of particles in Vicsek model */ #define VICSEK_VMAX 40.0 /* minimal speed of particles in Vicsek model */ #define COULOMB_LJ_FACTOR 1.0 /* relative intensity of LJ interaction in I_COULOMB_LJ interaction (default: 0.01) */ @@ -279,8 +288,8 @@ #define MU_XI 0.005 /* friction constant in thermostat */ #define KSPRING_BOUNDARY 2.0e11 /* confining harmonic potential outside simulation region */ #define KSPRING_OBSTACLE 1.0e9 /* harmonic potential of obstacles */ -#define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ -#define GRAVITY 4000.0 /* gravity acting on all particles */ +#define NBH_DIST_FACTOR 4.0 /* radius in which to count neighbours */ +#define GRAVITY 5000.0 /* gravity acting on all particles */ #define GRAVITY_X 0.0 /* horizontal gravity acting on all particles */ #define CIRCULAR_GRAVITY 0 /* set to 1 to have gravity directed to center */ #define INCREASE_GRAVITY 0 /* set to 1 to increase gravity during the simulation */ @@ -292,24 +301,29 @@ #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 100000.0 /* value of electric field */ +#define EFIELD 30000.0 /* value of electric field */ #define ADD_BFIELD 0 /* set to 1 to add a magnetic field */ -#define BFIELD 2.666666667 /* value of magnetic field */ +#define BFIELD 225.0 /* value of magnetic field */ #define CHARGE 1.0 /* charge of particles of first type */ #define CHARGE_B 1.0 /* charge of particles of second type */ #define INCREASE_E 0 /* set to 1 to increase electric field */ #define EFIELD_FACTOR 5000000.0 /* factor by which to increase electric field */ #define INCREASE_B 0 /* set to 1 to increase magnetic field */ -#define BFIELD_FACTOR 20000.0 /* factor by which to increase magnetic field */ +#define BFIELD_FACTOR 1.0 /* factor by which to increase magnetic field */ #define CHARGE_OBSTACLES 0 /* set to 1 for obstacles to be charged */ #define OBSTACLE_CHARGE 3.0 /* charge of obstacles */ #define KCOULOMB_OBSTACLE 1000.0 /* Coulomb force constant for charged obstacles */ +#define BFIELD_REGION 1 /* space-dependent magnetic field (0 for constant) */ + +#define ADD_WIND 1 /* set to 1 to add a "wind" friction force */ +#define WIND_FORCE 1.35e6 /* force of wind */ +#define WIND_YMIN -0.6 /* min altitude of region with wind */ #define ROTATION 1 /* set to 1 to include rotation of particles */ #define COUPLE_ANGLE_TO_THERMOSTAT 0 /* set to 1 to couple angular degrees of freedom to thermostat */ #define DIMENSION_FACTOR 1.0 /* scaling factor taking into account number of degrees of freedom */ -#define KTORQUE 5.0e3 /* force constant in angular dynamics */ -#define KTORQUE_BOUNDARY 5.0e5 /* constant in torque from the boundary */ +#define KTORQUE 1.0e5 /* force constant in angular dynamics */ +#define KTORQUE_BOUNDARY 1.0e5 /* constant in torque from the boundary */ #define KTORQUE_B 10.0 /* force constant in angular dynamics */ #define KTORQUE_DIFF 500.0 /* force constant in angular dynamics for different particles */ #define DRAW_SPIN 0 /* set to 1 to draw spin vectors of particles */ @@ -341,7 +355,7 @@ #define CENTER_VIEW_ON_OBSTACLE 0 /* set to 1 to center display on moving obstacle */ #define RESAMPLE_Y 0 /* set to 1 to resample y coordinate of moved particles (for shock waves) */ #define NTRIALS 2000 /* number of trials when resampling */ -#define OBSTACLE_RADIUS 0.015 /* radius of obstacle for circle boundary conditions */ +#define OBSTACLE_RADIUS 0.018 /* radius of obstacle for circle boundary conditions */ #define FUNNEL_WIDTH 0.25 /* funnel width for funnel boundary conditions */ #define OBSTACLE_XMIN 0.0 /* initial position of obstacle */ #define OBSTACLE_XMAX 3.0 /* final position of obstacle */ @@ -365,19 +379,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 0 /* set to 1 to add 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 50 /* time at which to add first particle */ -#define ADD_PERIOD 2 /* time interval between adding further particles */ -#define N_ADD_PARTICLES 8 /* number of particles to add */ -#define FINAL_NOADD_PERIOD 4700 /* final period where no particles are added */ +#define ADD_TIME 0 /* time at which to add first particle */ +#define ADD_PERIOD 15 /* time interval between adding further particles */ +#define N_ADD_PARTICLES 1 /* number of particles to add */ +#define FINAL_NOADD_PERIOD 500 /* final period where no particles are added */ #define SAFETY_FACTOR 4.0 /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */ -#define TRACER_PARTICLE 0 /* set to 1 to have a tracer particle */ -#define N_TRACER_PARTICLES 3 /* number of tracer particles */ -#define TRAJECTORY_LENGTH 8000 /* length of recorded trajectory */ -#define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */ -#define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */ +#define TRACER_PARTICLE 1 /* set to 1 to have a tracer particle */ +#define N_TRACER_PARTICLES 500 /* number of tracer particles */ +#define TRAJECTORY_LENGTH 4800 /* length of recorded trajectory */ +#define TRACER_LUM_FACTOR 5.0 /* controls luminosity decrease of trajectories with time */ +#define TRACER_PARTICLE_MASS 4.0 /* relative mass of tracer particle */ +#define TRAJECTORY_WIDTH 3 /* width of tracer particle trajectory */ #define ROTATE_BOUNDARY 0 /* set to 1 to rotate the repelling segments */ #define SMOOTH_ROTATION 1 /* set to 1 to update segments at each time step (rather than at each movie frame) */ @@ -463,7 +478,7 @@ #define LID_WIDTH 0.1 /* width of lid for BC_RECTANGLE_LID b.c. */ #define WALL_MASS 2000.0 /* mass of wall for BC_RECTANGLE_WALL b.c. */ #define WALL_FRICTION 0.0 /* friction on wall for BC_RECTANGLE_WALL b.c. */ -#define WALL_WIDTH 0.025 /* width of wall for BC_RECTANGLE_WALL b.c. */ +#define WALL_WIDTH 0.015 /* width of wall for BC_RECTANGLE_WALL b.c. */ #define WALL_VMAX 100.0 /* max speed of wall */ #define WALL_TIME 0 /* time during which to keep wall */ @@ -514,8 +529,8 @@ #define FLOOR_OMEGA 0 /* set to 1 to limit particle momentum to PMAX */ #define PMAX 1000.0 /* maximal force */ -#define HASHX 100 /* size of hashgrid in x direction */ -#define HASHY 88 /* size of hashgrid in y direction */ +#define HASHX 40 /* size of hashgrid in x direction */ +#define HASHY 20 /* size of hashgrid in y direction */ #define HASHMAX 100 /* maximal number of particles per hashgrid cell */ #define HASHGRID_PADDING 0.1 /* padding of hashgrid outside simulation window */ @@ -534,7 +549,10 @@ #define COUNT_PARTNER_TYPE ((RD_REACTION == CHEM_H2O_H_OH)||(RD_REACTION == CHEM_2H2O_H3O_OH)) #define PAIR_FORCE ((PAIR_PARTICLES)||((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_AGGREGATION)||(RD_REACTION == CHEM_AGGREGATION_CHARGE)||(RD_REACTION == CHEM_AGGREGATION_NNEIGH)||(RD_REACTION == CHEM_POLYGON_AGGREGATION)))) #define COMPUTE_PAIR_TORQUE (KTORQUE_PAIR_ANGLE != 0.0) -#define ADD_CONVEYOR_FORCE ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_BELT)||(SEGMENT_PATTERN == S_TWO_CONVEYOR_BELTS)||(SEGMENT_PATTERN == S_PERIODIC_CONVEYORS)||(SEGMENT_PATTERN == S_TEST_CONVEYORS))) +#define ADD_CONVEYOR_FORCE ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_BELT)||(SEGMENT_PATTERN == S_TWO_CONVEYOR_BELTS)||(SEGMENT_PATTERN == S_PERIODIC_CONVEYORS)||(SEGMENT_PATTERN == S_TEST_CONVEYORS)||(SEGMENT_PATTERN == S_CONVEYOR_SHOVELS)||(SEGMENT_PATTERN == S_CONVEYOR_MIXED)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_B)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_LONG))) +#define MOVE_CONVEYOR_BELT ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_SHOVELS)||(SEGMENT_PATTERN == S_CONVEYOR_MIXED)||(SEGMENT_PATTERN == S_CONVEYOR_SIEVE_LONG))) +#define ROTATE_OBSTACLES ((ADD_FIXED_OBSTACLES)&&((OBSTACLE_PATTERN == O_SIEVE)||(OBSTACLE_PATTERN == O_SIEVE_B)||(OBSTACLE_PATTERN == O_SIEVE_LONG))) +#define POLYGON_INTERACTION ((INTERACTION == I_POLYGON)||(INTERACTION == I_POLYGON_ALIGN)) double xshift = 0.0; /* x shift of shown window */ double xspeed = 0.0; /* x speed of obstacle */ @@ -968,7 +986,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase) { - double a, totalenergy = 0.0, damping, direction, dmean; + double a, totalenergy = 0.0, damping, damping1, damping_rot1, direction, dmean, dratio; static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0; int j, move, ncoup; @@ -983,7 +1001,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA particle[j].vx = px[j] + 0.5*DT_PARTICLE*particle[j].fx; particle[j].vy = py[j] + 0.5*DT_PARTICLE*particle[j].fy; particle[j].omega = pangle[j] + 0.5*DT_PARTICLE*particle[j].torque; - + px[j] = particle[j].vx + 0.5*DT_PARTICLE*particle[j].fx; py[j] = particle[j].vy + 0.5*DT_PARTICLE*particle[j].fy; pangle[j] = particle[j].omega + 0.5*DT_PARTICLE*particle[j].torque; @@ -1043,19 +1061,26 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA } move = 0; + damping1 = damping; + damping_rot1 = DAMPING_ROT; for (j=0; j 0.0) + { + damping1 = damping*particle[j].damping; + damping_rot1 = DAMPING_ROT*particle[j].damping; + } if ((THERMOSTAT_ON)&&(particle[j].thermostat)) { px[j] *= exp(- 0.5*DT_PARTICLE*xi); py[j] *= exp(- 0.5*DT_PARTICLE*xi); - if (!COUPLE_ANGLE_TO_THERMOSTAT) pangle[j] *= exp(- DT_PARTICLE*DAMPING_ROT); + if (!COUPLE_ANGLE_TO_THERMOSTAT) pangle[j] *= exp(- DT_PARTICLE*damping_rot1); } else { - px[j] *= exp(- DT_PARTICLE*damping); - py[j] *= exp(- DT_PARTICLE*damping); - pangle[j] *= exp(- DT_PARTICLE*DAMPING_ROT); + px[j] *= exp(- DT_PARTICLE*damping1); + py[j] *= exp(- DT_PARTICLE*damping1); + pangle[j] *= exp(- DT_PARTICLE*damping_rot1); // printf("Damping particle angular velocity\n"); } if ((THERMOSTAT_ON)&&(COUPLE_ANGLE_TO_THERMOSTAT)&&(particle[j].thermostat)) @@ -1542,13 +1567,13 @@ 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; + angle, theta, sum, alpha, bfield; double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds; double *cqx, *cqy, *cpx, *cpy, *cqangle, *cpangle; int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, cl, min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0, - group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt; + group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt, obs; int *particle_numbers; static int imin, imax; static short int first = 1; @@ -1702,7 +1727,9 @@ void animation() if (INCREASE_KREPEL) params.krepel = repel_schedule(i); if (INCREASE_BETA) params.beta = temperature_schedule(i); if (INCREASE_E) params.efield = efield_schedule(i); + else params.efield = EFIELD; if (INCREASE_B) params.bfield = bfield_schedule(i); + else params.bfield = BFIELD; 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) @@ -1918,8 +1945,17 @@ void animation() /* add magnetic force */ if (ADD_BFIELD) { - particle[j].fx += params.bfield*particle[j].charge*particle[j].vy*particle[j].mass_inv; - particle[j].fy -= params.bfield*particle[j].charge*particle[j].vx*particle[j].mass_inv; + bfield = params.bfield; + if (BFIELD_REGION != BF_CONST) + bfield *= (double)partial_bfield(particle[j].xc, particle[j].yc); + particle[j].fx += bfield*particle[j].charge*particle[j].vy*particle[j].mass_inv; + particle[j].fy -= bfield*particle[j].charge*particle[j].vx*particle[j].mass_inv; + } + + /* add wind force */ + if ((ADD_WIND)&&(particle[j].yc > WIND_YMIN)) + { + particle[j].fx += WIND_FORCE*particle[j].radius*particle[j].mass_inv; } if (FLOOR_FORCE) @@ -1968,6 +2004,22 @@ void animation() if (ADD_CONVEYOR_FORCE) for (belt = 0; belt < nbelts; belt++) conveyor_belt[belt].position += conveyor_belt[belt].speed*DT_PARTICLE; + if (MOVE_CONVEYOR_BELT) + update_conveyor_belts(segment, conveyor_belt); + + if (ROTATE_OBSTACLES) for (obs = 0; obs < nobstacles; obs++) + { + /* TEST */ + obstacle[obs].omega = obstacle[obs].omega0*(1.0 + 0.5*cos(i*DPI/200.0)); + obstacle[obs].angle -= obstacle[obs].omega*DT_PARTICLE; + if ((RATTLE_OBSTACLES)&&(obstacle[obs].oscillate)) + { + theta = obstacle[obs].phase + DPI*(double)i/(double)obstacle[obs].period; + obstacle[obs].xc = obstacle[obs].xc0 + obstacle[obs].amplitude*cos(theta); + obstacle[obs].yc = obstacle[obs].yc0 + obstacle[obs].amplitude*sin(theta); + } + } + if ((MOVE_SEGMENT_GROUPS)&&(i > INITIAL_TIME + SEGMENT_DEACTIVATION_TIME)) evolve_segment_groups(segment, i, segment_group); // if ((MOVE_SEGMENT_GROUPS)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segment_groups(segment, i, segment_group); @@ -2049,7 +2101,7 @@ void animation() /* update tracer particle trajectory */ if ((TRACER_PARTICLE)&&(i > INITIAL_TIME)) { - for (j=0; j #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 */ @@ -50,7 +50,7 @@ #define WINHEIGHT 1150 /* window height */ #define NX 960 /* number of grid points on x axis */ #define NY 575 /* number of grid points on y axis */ -#define HRES 1 /* factor for high resolution plots */ +#define HRES 4 /* factor for high resolution plots */ #define XMIN -2.0 #define XMAX 2.0 /* x interval */ @@ -64,7 +64,7 @@ #define NLAPLACIANS 0 /* number of fields for which to compute Laplacian */ #define SPHERE 1 /* set to 1 to simulate equation on sphere */ -#define DPOLE 1 /* safety distance to poles */ +#define DPOLE 0 /* safety distance to poles */ #define DSMOOTH 1 /* size of neighbourhood of poles that are smoothed */ #define SMOOTHPOLE 0.05 /* smoothing coefficient at poles */ #define SMOOTHCOTPOLE 0.05 /* smoothing coefficient of cotangent at poles */ @@ -72,21 +72,21 @@ #define SMOOTHBLOCKS 1 /* set to 1 to use blocks of points near the poles */ #define BLOCKDIST 64 /* distance to poles where points are blocked */ #define ZERO_MERIDIAN 190.0 /* choice of zero meridian (will be at left/right boundary of 2d plot) */ -#define POLE_NODRAW 10 /* distance around poles where wave is not drawn */ +#define POLE_NODRAW 2 /* distance around poles where wave is not drawn */ #define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */ #define ADD_MAGNETIC_FIELD 0 /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */ -#define ADD_FORCE_FIELD 0 /* set to 1 to add a foce field (for compressible Euler equation) */ +#define ADD_FORCE_FIELD 1 /* set to 1 to add a foce field (for compressible Euler equation) */ #define POTENTIAL 7 /* type of potential or vector potential, see list in global_3d.c */ #define FORCE_FIELD 6 /* type of force field, see list in global_3d.c */ #define ADD_CORIOLIS_FORCE 1 /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */ #define VARIABLE_DEPTH 1 /* set to 1 for variable depth in shallow water equation */ -#define SWATER_DEPTH 5 /* variable depth in shallow water equation */ +#define SWATER_DEPTH 10 /* variable depth in shallow water equation */ #define ANTISYMMETRIZE_WAVE_FCT 0 /* set tot 1 to make wave function antisymmetric */ -#define ADAPT_STATE_TO_BC 0 /* to smoothly adapt initial state to obstacles */ +#define ADAPT_STATE_TO_BC 1 /* to smoothly adapt initial state to obstacles */ #define OBSTACLE_GEOMETRY 84 /* geometry of obstacles, as in B_DOMAIN */ -#define BC_STIFFNESS 50.0 /* controls region of boundary condition control */ +#define BC_STIFFNESS 0.25 /* controls region of boundary condition control */ #define CHECK_INTEGRAL 1 /* set to 1 to check integral of first field */ #define JULIA_SCALE 0.5 /* scaling for Julia sets */ @@ -140,7 +140,8 @@ #define VISCOSITY 0.02 #define POISSON_STIFFNESS 1.0 /* stiffness of Poisson equation solver for incompressible Euler */ -#define DISSIPATION 1.0e-8 +#define DISSIPATION 0.0 +#define DISSIPATION_EXT 1.0e-3 /* dissipation of shallow water eq. outside domain */ #define RPSA 0.75 /* parameter in Rock-Paper-Scissors-type interaction */ #define RPSLZB 0.0 /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */ @@ -156,25 +157,29 @@ #define BZQ 0.0008 /* parameter in BZ equation */ #define BZF 1.2 /* parameter in BZ equation */ #define B_FIELD 10.0 /* magnetic field */ -#define G_FIELD 0.002 /* gravity/Coriolis force */ -#define BC_FIELD 1.0e-5 /* constant in repulsive field from obstacles */ +#define G_FIELD 0.025 /* gravity/Coriolis force */ +#define BC_FIELD 1.0e-6 /* constant in repulsive field from obstacles */ #define AB_RADIUS 0.2 /* radius of region with magnetic field for Aharonov-Bohm effect */ #define K_EULER 50.0 /* constant in stream function integration of Euler equation */ #define K_EULER_INC 0.5 /* constant in incompressible Euler equation */ #define C_EULER_COMP 0.1 /* constant in compressible Euler equation */ +#define SWATER_VARDEPTH_FACTOR 0.1 /* force constant in front of variable depth term */ +#define SWATER_CORIOLIS_FORCE 0.002 /* Coriolis force for shallow water equation */ #define SMOOTHEN_VORTICITY 0 /* set to 1 to smoothen vorticity field in Euler equation */ #define SMOOTHEN_VELOCITY 1 /* set to 1 to smoothen velocity field in Euler equation */ #define SMOOTHEN_PERIOD 7 /* period between smoothenings */ -#define SMOOTH_FACTOR 0.15 /* factor by which to smoothen */ +#define SMOOTH_FACTOR 0.2 /* factor by which to smoothen */ #define ADD_OSCILLATING_SOURCE 0 /* set to 1 to add an oscillating wave source */ #define OSCILLATING_SOURCE_PERIOD 1 /* period of oscillating source */ #define OSCILLATING_SOURCE_OMEGA 0.2 /* frequency of oscillating source */ -#define ADD_TRACERS 1 /* set to 1 to add tracer particles (for Euler equations) */ +#define ADD_TRACERS 1 /* set to 1 tof add tracer particles (for Euler equations) */ #define N_TRACERS 2000 /* number of tracer particles */ #define TRACERS_STEP 0.1 /* step size in tracer evolution */ +#define RESPAWN_TRACERS 1 /* set to 1 to randomly move tracer position */ +#define RESPAWN_PROBABILTY 0.005 /* probability of moving tracer */ #define T_OUT 2.0 /* outside temperature */ #define T_IN 0.0 /* inside temperature */ @@ -208,11 +213,16 @@ #define LAMINAR_FLOW_YPERIOD 1.0 /* period of laminar flow in y direction */ #define PRESSURE_GRADIENT 0.2 /* amplitude of pressure gradient for Euler equation */ -#define SWATER_MIN_HEIGHT 0.025 /* min height of initial condition for shallow water equation */ -#define DEPTH_FACTOR 0.075 /* proportion of min height in variable depth */ -#define TANH_FACTOR 1.0 /* steepness of variable depth */ +#define SWATER_MIN_HEIGHT 0.05 /* min height of initial condition for shallow water equation */ +#define DEPTH_FACTOR 0.00002 /* proportion of min height in variable depth */ +#define TANH_FACTOR 5.0 /* steepness of variable depth */ #define EULER_GRADIENT_YSHIFT 0.0 /* y-shift in computation of gradient in Euler equation */ +#define ADD_MOON_FORCING 1 /* set to 1 to simulate tidal forces from Moon */ +#define FORCING_AMP 5.0e-6 /* amplitude of periodic forcing */ +#define FORCING_CONST_AMP 0.0 /* amplitude of periodic forcing */ +#define FORCING_PERIOD 1600 /* period of forcing */ +#define FORCING_SHIFT 1.5 /* phase shift of forcing in units of Pi */ /* Boundary conditions, see list in global_pdes.c */ @@ -225,8 +235,8 @@ /* Parameters for length and speed of simulation */ -#define NSTEPS 950 /* number of frames of movie */ -#define NVID 85 /* number of iterations between images displayed on screen */ +#define NSTEPS 1400 /* number of frames of movie */ +#define NVID 80 /* 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 */ @@ -248,13 +258,14 @@ #define PLOT_SPHERE 1 /* draws fields on a sphere */ #define ROTATE_VIEW 1 /* set to 1 to rotate position of observer */ -#define ROTATE_ANGLE 360.0 /* total angle of rotation during simulation */ +#define ROTATE_ANGLE -45.0 /* total angle of rotation during simulation */ #define SHADE_3D 1 /* set to 1 to change luminosity according to normal vector */ #define SHADE_2D 0 /* set to 1 to change luminosity according to normal vector */ #define VIEWPOINT_TRAJ 1 /* type of viewpoint trajectory */ #define MAX_LATITUDE 45.0 /* maximal latitude for viewpoint trajectory VP_ORBIT2 */ #define DRAW_PERIODICISED 0 /* set to 1 to repeat wave periodically in x and y directions */ +#define DRAW_MOON_POSITION 0 /* set to 1 to draw vertical lines for tidal sim */ /* Plot type - color scheme */ @@ -275,8 +286,12 @@ #define ADD_POTENTIAL_TO_Z 0 /* set to 1 to add the external potential to z-coordinate of plot */ #define ADD_POT_CONSTANT 0.35 /* constant added to wave height */ #define DRAW_DEPTH 0 /* set to 1 to draw water depth */ -#define DEPTH_SCALE 0.75 /* vertical scaling of depth plot */ +#define DEPTH_SCALE 0.75 /* vertical scaling of depth plot */ #define DEPTH_SHIFT -0.015 /* vertical shift of depth plot */ +#define FLOODING 1 /* set to 1 for drawing water when higher than continents */ +// #define FLOODING_VSHIFT 0.51 /* controls when wave is considered higher than land */ +#define FLOODING_VSHIFT 0.56 /* controls when wave is considered higher than land */ +#define FLOODING_VSHIFT_2D 0.61 /* controls when wave is considered higher than land */ #define PLOT_SCALE_ENERGY 0.05 /* vertical scaling in energy plot */ @@ -306,8 +321,8 @@ /* Color schemes, see list in global_pdes.c */ -#define COLOR_PALETTE 10 /* Color palette, see list in global_pdes.c */ -#define COLOR_PALETTE_B 12 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE 11 /* Color palette, see list in global_pdes.c */ +#define COLOR_PALETTE_B 16 /* Color palette, see list in global_pdes.c */ #define BLACK 1 /* black background */ #define COLOR_OUT_R 1.0 /* color outside domain */ @@ -330,8 +345,8 @@ #define VSCALE_PRESSURE 2.0 /* additional scaling factor for color scheme Z_EULER_PRESSURE */ #define PRESSURE_SHIFT 10.0 /* shift for color scheme Z_EULER_PRESSURE */ #define PRESSURE_LOG_SHIFT -2.5 /* shift for color scheme Z_EULER_PRESSURE */ -#define VSCALE_WATER_HEIGHT 30.0 /* vertical scaling of water height */ -#define ADD_HEIGHT_CONSTANT -0.025 /* constant added to wave height */ +#define VSCALE_WATER_HEIGHT 15.0 /* vertical scaling of water height */ +#define ADD_HEIGHT_CONSTANT -0.016 /* constant added to wave height */ #define SHADE_SCALE_2D 0.25 /* controls "depth" of 2D shading */ #define COLORHUE 260 /* initial hue of water color for scheme C_LUM */ @@ -354,7 +369,7 @@ #define ZSCALE_SPEED 300.0 /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ #define ZSHIFT_SPEED 0.0 /* additional shift of z-coord Z_EULER_SPEED and Z_SWATER_SPEED */ #define ZSCALE_NORMGRADIENT -0.0001 /* vertical scaling for Z_NORM_GRADIENT */ -#define VSCALE_SWATER 50.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ +#define VSCALE_SWATER 60.0 /* additional scaling factor for color scheme Z_EULER_DENSITY */ #define NXMAZE 7 /* width of maze */ #define NYMAZE 7 /* height of maze */ @@ -364,7 +379,7 @@ #define MAZE_WIDTH 0.04 /* half width of maze walls */ #define DRAW_COLOR_SCHEME 1 /* set to 1 to plot the color scheme */ -#define COLORBAR_RANGE 2.5 /* scale of color scheme bar */ +#define COLORBAR_RANGE 3.0 /* 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 */ @@ -427,35 +442,35 @@ double u_3d[2] = {0.75, -0.45}; /* projections of basis vectors for REP_AXO_3D representation */ double v_3d[2] = {-0.75, -0.45}; double w_3d[2] = {0.0, 0.015}; -double light[3] = {0.40824829, -0.816496581, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ -double observer[3] = {8.0, 4.0, 2.5}; /* location of observer for REP_PROJ_3D representation */ +double light[3] = {-0.816496581, 0.40824829, 0.40824829}; /* vector of "light" direction for P_3D_ANGLE color scheme */ +double observer[3] = {-6.0, -6.0, 4.5}; /* location of observer for REP_PROJ_3D representation */ int reset_view = 0; /* switch to reset 3D view parameters (for option ROTATE_VIEW) */ /* constants for simulations on planets */ #define ADD_DEM 1 /* add DEM (digital elevation model) */ -#define ADD_NEGATIVE_DEM 0 /* add DEM with bathymetric data */ -#define RSCALE_DEM 0.1 /* scaling factor of radial component for DEM */ -#define SMOOTH_DEM 0 /* set to 1 to smoothen DEM (to make altitude less constant) */ -#define DEM_SMOOTH_STEPS 1 /* number of smoothening steps */ +#define ADD_NEGATIVE_DEM 1 /* add DEM with bathymetric data */ +#define RSCALE_DEM 0.075 /* scaling factor of radial component for DEM */ +#define SMOOTH_DEM 1 /* set to 1 to smoothen DEM (to make altitude less constant) */ +#define DEM_SMOOTH_STEPS 10 /* number of smoothening steps */ #define DEM_SMOOTH_HEIGHT 2.0 /* relative height below which to smoothen */ #define DEM_MAXHEIGHT 9000.0 /* max height of DEM (estimated from Everest/Olympus Mons) */ #define DEM_MAXDEPTH -10000 /* max depth of DEM */ -#define PLANET_SEALEVEL 2500.0 /* sea level for flooded planet */ +#define PLANET_SEALEVEL 0.0 /* sea level for flooded planet */ #define VENUS_NODATA_FACTOR 0.5 /* altitude to assign to DEM points without data (fraction of mean altitude) */ #define Z_SCALING_FACTOR 0.8 /* overall scaling factor of z axis for REP_PROJ_3D representation */ -#define XY_SCALING_FACTOR 2.0 /* overall scaling factor for on-screen (x,y) coordinates after projection */ +#define XY_SCALING_FACTOR 2.1 /* 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 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 RSCALE 0.2 /* scaling factor of radial component */ -#define RSHIFT -0.01 /* shift in radial component */ +#define RSCALE 0.15 /* scaling factor of radial component */ +#define RSHIFT -0.075 /* shift in radial component */ #define RMAX 2.0 /* max value of radial component */ #define RMIN 0.5 /* min value of radial component */ -#define COS_VISIBLE -0.75 /* limit on cosine of normal to shown facets */ +#define COS_VISIBLE -0.35 /* limit on cosine of normal to shown facets */ /* For debugging purposes only */ #define FLOOR 1 /* set to 1 to limit wave amplitude to VMAX */ @@ -670,6 +685,8 @@ void compute_gfield(int i, int j, double *gx, double *gy) } } } + + void initialize_potential(double potential_field[NX*NY]) /* initialize the potential field, e.g. for the Schrödinger equation */ { @@ -744,9 +761,18 @@ void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc #pragma omp parallel for private(i,j) for (i=1; i 0.0) + { + gfield[i*NY+j] = -BC_FIELD*(wsphere[(i+1)*NY+j].altitude - wsphere[(i-1)*NY+j].altitude)/dx; + gfield[NX*NY+i*NY+j] = -BC_FIELD*(wsphere[i*NY+j+1].altitude - wsphere[i*NY+j-1].altitude)/dy; + } + else + { + gfield[i*NY+j] = 0.0; + gfield[NX*NY+i*NY+j] = 0.0; + } } } @@ -846,7 +872,7 @@ void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc void compute_water_depth(int i, int j, t_rde *rde, t_wave_sphere wsphere[NX*NY], int ncircles) /* initialize the vector potential, for Schrodinger equation in a magnetic field */ { - double x, y, xy[2], r0, r, z, h, h0, hh, x1, y1, z1, d, rmax, rmin, phi; + double x, y, xy[2], r0, r, z, h, h0, hh, x1, y1, z1, d, rmax, rmin, phi, vshift; int n, nx, ny; static double phi1, phi2, phi3, sq3, rico; static int first = 1; @@ -998,16 +1024,30 @@ void compute_water_depth(int i, int j, t_rde *rde, t_wave_sphere wsphere[NX*NY], rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; break; } + case (SH_EARTH): + { + /* TODO */ + vshift = PLANET_SEALEVEL/(DEM_MAXHEIGHT - DEM_MAXDEPTH); + h = 0.5 - RSCALE_DEM*(wsphere[i*NY+j].altitude - vshift); + rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h; + break; + } } } -double initialize_water_depth(t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) +double initialize_water_depth(t_rde rde[NX*NY], t_wave_sphere *wsphere, t_wave_sphere *wsphere_hr) { int i, j, ncircles; double dx, dy, min, max, pscal, norm, vz = 0.01; if (SWATER_DEPTH == SH_CIRCLES) ncircles = init_circle_config_pattern(circles, CIRCLE_PATTERN); + if (SWATER_DEPTH == SH_EARTH) + { + init_earth_map_rde(wsphere, 1); + init_earth_map_rde(wsphere_hr, HRES); + } + #pragma omp parallel for private(i,j) for (i=0; i= NX) moon_position -= NX; + +// i = NX/2; + i = moon_position; + printf("Phase = %.5lg, Forcing at i = %i: %.5lg, Moon position = %i\n", phase, i, wsphere[i*NY+NY/2].force, moon_position); +} + + void evolve_wave_half(double *phi_in[NFIELDS], double *phi_out[NFIELDS], short int xy_in[NX*NY], double potential_field[NX*NY], double vector_potential_field[2*NX*NY], double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) /* time step of field evolution */ @@ -1485,6 +1549,18 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) // printf("etax = %.3lg, etay = %.3lg, ux = %.3lg, uy = %.3lg, vx = %.3lg, vy = %.3lg\n", etax, etay, ux, uy, vx, vy); + /* do not evolve wave at high altitude */ + if (FLOODING) + { + if (!wsphere[i*NY+j].evolve_wave) + { + phi_out[0][i*NY+j] = eta; + phi_out[1][i*NY+j] = 0.0; + phi_out[2][i*NY+j] = 0.0; + break; + } + } + if (SPHERE) { /* TODO */ @@ -1492,11 +1568,12 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) { phi_out[0][i*NY+j] = eta - intstep*(u*etax + v*etay + eta*(ux + vy)); phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax); - phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etax); + phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etay); } else { - phi_out[0][i*NY+j] = SWATER_MIN_HEIGHT; +// phi_out[0][i*NY+j] = SWATER_MIN_HEIGHT; + phi_out[0][i*NY+j] = eta; phi_out[1][i*NY+j] = 0.0; phi_out[2][i*NY+j] = 0.0; } @@ -1507,7 +1584,11 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax); phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etay); } - + if (ADD_CORIOLIS_FORCE) + { + phi_out[1][i*NY+j] += intstep*SWATER_CORIOLIS_FORCE*v*wsphere[i*NY+j].ctheta; + phi_out[2][i*NY+j] -= intstep*SWATER_CORIOLIS_FORCE*u*wsphere[i*NY+j].reg_cottheta; + } if (VISCOSITY > 0.0) { phi_out[1][i*NY+j] += intstep*VISCOSITY*delta_u[i*NY+j]; @@ -1518,18 +1599,32 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY]) phi_out[1][i*NY+j] -= intstep*DISSIPATION*u; phi_out[2][i*NY+j] -= intstep*DISSIPATION*v; } + if ((DISSIPATION_EXT > 0.0)&&(!wsphere[i*NY+j].indomain)) + { + phi_out[1][i*NY+j] -= intstep*DISSIPATION_EXT*u; + phi_out[2][i*NY+j] -= intstep*DISSIPATION_EXT*v; + } if (ADD_FORCE_FIELD) { phi_out[1][i*NY+j] += intstep*gfield[i*NY+j]; phi_out[2][i*NY+j] += intstep*gfield[NX*NY+i*NY+j]; } - if (VARIABLE_DEPTH) + if (ADD_MOON_FORCING) { - phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].depth*(ux + vy); - phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].gradx*u; - phi_out[0][i*NY+j] -= intstep*rde[i*NY+j].grady*v; + phi_out[0][i*NY+j] += intstep*wsphere[i*NY+j].force; + } + if ((VARIABLE_DEPTH)&&(j > DSMOOTH)&&(j < NY-DSMOOTH)) + { + phi_out[0][i*NY+j] -= intstep*SWATER_VARDEPTH_FACTOR*rde[i*NY+j].depth*(ux + vy); + phi_out[0][i*NY+j] -= intstep*SWATER_VARDEPTH_FACTOR*rde[i*NY+j].gradx*u; + phi_out[0][i*NY+j] -= intstep*SWATER_VARDEPTH_FACTOR*rde[i*NY+j].grady*v; } + if ((j <= DSMOOTH)||(j >= NY-DSMOOTH)) + { +// printf("Pole reset\n"); + phi_out[0][i*NY+j] = SWATER_MIN_HEIGHT; + } break; } } @@ -1605,7 +1700,15 @@ void update_tracer_table(double tracers[2*N_TRACERS*NSTEPS], t_rde rde[NX*NY], i /* update tracer information in rde */ { int tracer, t, t1, maxtime, i, j, n, ij[2], length = 50, cell, oldcell; - double x, y; + double x, y, x1, y1, dist; + static double maxlength; + static int first = 1; + + if (first) + { + maxlength = pow(0.01*(XMAX - XMIN), 2.0); + first = 0; + } #pragma omp parallel for private(cell) for (cell=0; cell 0)) + dist = (x-x1)*(x-x1) + (y-y1)*(y-y1); + if ((cell != oldcell)&&(t > 0)&&(dist < maxlength)) rde[cell].prev_cell = oldcell; oldcell = cell; + x1 = x; + y1 = y; } } } @@ -1654,6 +1762,7 @@ void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], t_ { int tracer, i, j, n, t, ij[2], iplus, jplus, prev_cell, new_cell; double x, y, xy[2], vx, vy; + static int n_respawn = 0; step = TRACERS_STEP; @@ -1717,7 +1826,14 @@ void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], t_ if (y < YMIN) y += (YMAX - YMIN); if (time+1 < NSTEPS) - { + { + if ((RESPAWN_TRACERS)&&((double)rand()/RAND_MAX < RESPAWN_PROBABILTY)) + { + x = XMIN + 0.05 + (XMAX - XMIN - 0.1)*rand()/RAND_MAX; + y = YMIN + 0.05 + (YMAX - YMIN - 0.1)*rand()/RAND_MAX; + n_respawn++; + printf("Respawning tracer %i, %i respawns\n", tracer, n_respawn); + } tracers[(time+1)*2*N_TRACERS + 2*tracer] = x; tracers[(time+1)*2*N_TRACERS + 2*tracer + 1] = y; } @@ -2024,7 +2140,7 @@ void animation() if (VARIABLE_DEPTH) { // water_depth = (t_swater_depth *)malloc(NX*NY*sizeof(t_swater_depth)); - max_depth = initialize_water_depth(rde, wsphere); + max_depth = initialize_water_depth(rde, wsphere, wsphere_hr); printf("Max depth = %.3lg\n", max_depth); } @@ -2119,8 +2235,15 @@ void animation() // init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in); - init_swater_laminar_flow_sphere(0, -0.75, 0.0075, 6, 1, 0.0025, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); - + +// init_swater_laminar_flow_sphere(0, 0.25, 0.00075, 4, 4, 0.0025, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); + +// init_tidal_state(1, 0.0015, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); + +// init_linear_blob_sphere(0, 1.3*PI, 0.65*PI, 0.0, 0.0, 0.3, 0.1, 0.1, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); + + init_expanding_blob_sphere(0, (220.0/180.0)*PI, (140.0/180.0)*PI, 1.0, 0.075, 0.1, SWATER_MIN_HEIGHT, phi, xy_in, wsphere); + // add_gaussian_wave(-1.6, -0.5, 0.015, 0.25, SWATER_MIN_HEIGHT, phi, xy_in); if (SMOOTHBLOCKS) for (k=0; k DPI) + { + segment[i].angle1 -= DPI; + segment[i].angle2 -= DPI; + } + while (segment[i].angle2 < 0.0) + { + segment[i].angle1 += DPI; + segment[i].angle2 += DPI; + } + } +} + + +int add_shovel_to_belt(double x, double y, double angle, double size, double width, t_segment segment[NMAXSEGMENTS], t_belt belt) +/* add showver to conveyor belt */ +{ + int i, n, iminus, iplus; + double sq2, angle1, angle2, cg, sg, t; + + if (nsegments + 6 > NMAXSEGMENTS) + { + printf("NMAXSEGMENTS too small"); + return(0); + } + + n = nsegments; +// sq2 = sqrt(2.0)/2.0; + cg = cos(PI/6.0); + sg = sin(PI/6.0); + t = size + width*(cg - 1.0)/sg; + + segment[n].x1 = 0.0; + segment[n].y1 = 0.0; + + segment[n+1].x1 = 0.0; + segment[n+1].y1 = size; + + segment[n+2].x1 = size*sg; + segment[n+2].y1 = size*(1.0+cg); + + segment[n+3].x1 = size*sg + width*cg; + segment[n+3].y1 = size*(1.0+cg) - width*sg; + + segment[n+4].x1 = width; + segment[n+4].y1 = size*(1.0 + cg) - width*sg - t*cg; + +// segment[n+2].x1 = size*sq2; +// segment[n+2].y1 = size*(1.0+sq2); +// +// segment[n+3].x1 = size*sq2 + width*sq2; +// segment[n+3].y1 = size*(1.0+sq2) - width*sq2; +// +// segment[n+4].x1 = width; +// segment[n+4].y1 = size + width*(1.0 - sqrt(2.0)); + + segment[n+5].x1 = width; + segment[n+5].y1 = 0.0; + + for (i=0; i<6; i++) + { + segment[n+i].x2 = segment[n+i+1].x1; + segment[n+i].y2 = segment[n+i+1].y1; + } + segment[n+6].x2 = segment[n].x1; + segment[n+6].y2 = segment[n].y1; + + for (i=n; i n + 5) iplus = n; + angle1 = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID; + angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID; + if (angle2 < angle1) angle2 += DPI; + segment[i].angle1 = angle1; + segment[i].angle2 = angle2; + } + + for (i=n; i n + 13) iplus = n; - angle1 = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID; - angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID; - if (angle2 < angle1) angle2 += DPI; - segment[i].angle1 = angle1; - segment[i].angle2 = angle2; - -// printf("i = %i, iplus = %i, iminus = %i, angle1 = %.0f, angle2 = %.0f\n", i, iplus, iminus, angle*360.0/DPI, angle2*360.0/DPI); - } + { + iminus = i-1; + iplus = i+1; + if (iminus < n) iminus = n + 13; + if (iplus > n + 13) iplus = n; + angle1 = argument(segment[iplus].x1 - segment[i].x1, segment[iplus].y1 - segment[i].y1) + PID; + angle2 = argument(segment[i].x1 - segment[iminus].x1, segment[i].y1 - segment[iminus].y1) + PID; + if (angle2 < angle1) angle2 += DPI; + segment[i].angle1 = angle1; + segment[i].angle2 = angle2; + } nsegments += 14; @@ -2039,8 +2374,67 @@ int add_conveyor_belt(double x1, double y1, double x2, double y2, double width, conveyor_belt[nbelts].angle = angle; conveyor_belt[nbelts].tx = tx; conveyor_belt[nbelts].ty = ty; + conveyor_belt[nbelts].nshovels = nshovels; nbelts++; + if (nshovels > NMAXSHOVELS) + { + printf("NMAXSHOVELS too small"); + conveyor_belt[nbelts].nshovels = 0; + return(0); + } + + if (nshovels > 0) + { + beltlength = 2.0*length + DPI*width; + shovel_dist = beltlength/(double)nshovels; + shovel_pos = 0.0; + shovel = 0; + + while (shovel_pos < length) + { + x = x1 - width*ty + shovel_pos*tx; + y = y1 + width*tx + shovel_pos*ty; + conveyor_belt[nbelts-1].shovel_segment[shovel] = nsegments; + add_shovel_to_belt(x, y, angle, width, 0.3*width, segment, conveyor_belt[nbelts-1]); + conveyor_belt[nbelts-1].shovel_pos[shovel] = shovel_pos; + shovel_pos += shovel_dist; + shovel++; + } + while (shovel_pos < length + width*PI) + { + beta = (shovel_pos - length)/width; + x = x2 + width*cos(angle + PID - beta); + y = y2 + width*sin(angle + PID - beta); + conveyor_belt[nbelts-1].shovel_segment[shovel] = nsegments; + add_shovel_to_belt(x, y, angle - beta, width, 0.3*width, segment, conveyor_belt[nbelts-1]); + conveyor_belt[nbelts-1].shovel_pos[shovel] = shovel_pos; + shovel_pos += shovel_dist; + shovel++; + } + while (shovel_pos < 2.0*length + width*PI) + { + x = x2 + width*ty - (shovel_pos - length - width*PI)*tx; + y = y2 - width*tx - (shovel_pos - length - width*PI)*ty; + conveyor_belt[nbelts-1].shovel_segment[shovel] = nsegments; + add_shovel_to_belt(x, y, angle - PI, width, 0.3*width, segment, conveyor_belt[nbelts-1]); + conveyor_belt[nbelts-1].shovel_pos[shovel] = shovel_pos; + shovel_pos += shovel_dist; + shovel++; + } + while (shovel_pos < beltlength) + { + beta = (shovel_pos - 2.0*length - width*PI)/width; + x = x1 + width*cos(angle - PID - beta); + y = y1 + width*sin(angle - PID - beta); + conveyor_belt[nbelts-1].shovel_segment[shovel] = nsegments; + add_shovel_to_belt(x, y, angle - PI - beta, width, 0.3*width, segment, conveyor_belt[nbelts-1]); + conveyor_belt[nbelts-1].shovel_pos[shovel] = shovel_pos; + shovel_pos += shovel_dist; + shovel++; + } + } + return(1); } @@ -2051,8 +2445,12 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[N double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, yy, dy, ca, sa, padding; double lfoot[NTREES][2], rfoot[NTREES][2]; - /* set default to no conveyor */ - for (i=0; i DPI) - { - segment[i].angle1 -= DPI; - segment[i].angle2 -= DPI; - } - while (segment[i].angle2 < 0.0) - { - segment[i].angle1 += DPI; - segment[i].angle2 += DPI; - } - } -} - /* Computation of interaction force */ double lennard_jones_force(double r, t_particle particle) @@ -4188,7 +4596,7 @@ void polygon_interaction(t_particle part_i, t_particle part_k, double ca, double /* computes force and torque of particle k on particle i */ { int s, k, s0, s1; - double x0, y0, x, y, f, fx, fy, cb, sb, cab, sab, r, r0, z2, phi, ckt, skt, d; + double x0, y0, x, y, f, fx, fy, cb, sb, cab, sab, r, r0, z2, phi, ckt, skt, d, fx1, fy1, dmax2; double torque, mu_i, mu_k, r3, width, gamma_beta, rot, cr, sr, xx[NPOLY], yy[NPOLY], dd[NPOLY], z, z0, z02, d1, dmin; static double theta, twotheta, ctheta, stheta, cornerx[NPOLY+1], cornery[NPOLY+1], nx[NPOLY], ny[NPOLY]; static int first = 1; @@ -4220,6 +4628,7 @@ void polygon_interaction(t_particle part_i, t_particle part_k, double ca, double torque = 0.0; mu_i = part_i.radius; mu_k = part_k.radius; + dmax2 = 4.0*mu_i*mu_i; width = 0.1*mu_i; mu_i -= width; gamma_beta = part_k.angle - part_i.angle; @@ -4250,7 +4659,7 @@ void polygon_interaction(t_particle part_i, t_particle part_k, double ca, double } /* compute force from vertices */ - for (s = 0; s width) f = width; - fx -= f*nx[k]; - fy -= f*ny[k]; - torque += x*fy - y*fx; + fx1 = f*nx[k]; + fy1 = f*ny[k]; + fx -= fx1; + fy -= fy1; + torque += x*fy1 - y*fx1; +// fx -= f*nx[k]; +// fy -= f*ny[k]; +// torque += x*fy - y*fx; } else for (s1=0; s1 width) f = width; - fx -= f*(x - mu_i*cornerx[s1])/d1; - fy -= f*(y - mu_i*cornery[s1])/d1; + fx1 -= f*(x - mu_i*cornerx[s1])/d1; + fy1 -= f*(y - mu_i*cornery[s1])/d1; } - torque += x*fy - y*fx; + fx += fx1; + fy += fy1; + torque += x*fy1 - y*fx1; } } } @@ -6046,19 +6464,27 @@ int add_particle(double x, double y, double vx, double vy, double mass, short in particle[i].energy = 0.0; particle[i].emean = 0.0; particle[i].dirmean = 0.0; - - if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); + + if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(RANDOM_RADIUS_MIN + RANDOM_RADIUS_RANGE*((double)rand()/RAND_MAX)); particle[i].mass_inv = 1.0/mass; if (particle[i].type == 0) particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; else particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; + + if ((RANDOM_RADIUS)&&(ADAPT_MASS_TO_RADIUS > 0.0)) + particle[i].mass_inv *= 1.0/(1.0 + pow(particle[i].radius/MU, ADAPT_MASS_TO_RADIUS)); + if ((RANDOM_RADIUS)&&(ADAPT_DAMPING_TO_RADIUS > 0.0)) + particle[i].damping = 1.0 + ADAPT_DAMPING_FACTOR*pow(particle[i].radius/MU, ADAPT_DAMPING_TO_RADIUS); + particle[i].vx = vx; particle[i].vy = vy; particle[i].energy = (particle[i].vx*particle[i].vx + particle[i].vy*particle[i].vy)*particle[i].mass_inv; particle[i].angle = DPI*(double)rand()/RAND_MAX; - particle[i].omega = 0.0; + particle[i].omega = OMEGA_INITIAL*gaussian(); + +// printf("Particle[%i].omega = %.4lg\n", i, particle[i].omega); // if (particle[i].type == 1) // { @@ -6407,6 +6833,16 @@ void compute_particle_colors(t_particle particle, t_cluster cluster[NMAXCIRCLES] *width = BOUNDARY_WIDTH; break; } + case (P_RADIUS): + { +// hue = atan(5.0*(particle.radius/MU - 0.75))/PID; + hue = (particle.radius/MU - RANDOM_RADIUS_MIN)/RANDOM_RADIUS_RANGE; +// hue = 0.5*(hue + 1.0); + hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue; + *radius = particle.radius; + *width = BOUNDARY_WIDTH; + break; + } } switch (plot) { @@ -7240,24 +7676,27 @@ 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) +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) /* draw tracer particle trajectory */ { - int i, j, time; - double x1, x2, y1, y2, rgb[3], lum; + int i, j, time, p, width; + double x1, x2, y1, y2, rgb[3], rgbx[3], rgby[3], radius, lum; // blank(); glLineWidth(TRAJECTORY_WIDTH); printf("drawing trajectory\n"); - for (j=0; j 0)) color_background(particle, bg_color, hashgrid); - else blank(); + else if (!TRACER_PARTICLE) blank(); glColor3f(1.0, 1.0, 1.0); @@ -7783,6 +8223,14 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES] draw_line(x - r, y, x + r, y); draw_line(x, y - r, x, y + r); } + if (ROTATE_OBSTACLES) + { + ca = cos(obstacle[i].angle); + sa = sin(obstacle[i].angle); +// printf("Obstacle %i angle = %.5lg\n", i, obstacle[i].angle); + draw_line(x - r*ca, y - r*sa, x + r*ca, y + r*sa); + draw_line(x + r*sa, y - r*ca, x - r*sa, y + r*ca); + } } } if (ADD_FIXED_SEGMENTS) @@ -8785,8 +9233,16 @@ void print_particles_speeds(t_particle particle[NMAXCIRCLES]) double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall, double krepel) { - int i, k; - double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin, charge, d2, speed; + int i, k, corner; + double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, y1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin, charge, d2, speed, pangle, distance1, sangle; + static int first = 1; + static double seqangle; + + if (first) + { + seqangle = PI/(double)NPOLY; + first = 0; + } /* compute force from fixed circular obstacles */ if (ADD_FIXED_OBSTACLES) for (i=0; i 0.0)&&(proj < 1.0)) { // distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; - r = 1.5*particle[j].radius; + /* case of interacting polygons */ +// if (POLYGON_INTERACTION) +// { +// // distance = segment[i].nx*x + segment[i].ny*y - segment[i].c; +// pangle = DPI/(double)NPOLY; +// dx = 0.0; +// dy = 0.0; +// for (corner=0; corner 0.0)&&(pvect < pvmin)) pvect = pvmin; else if ((pvect < 0.0)&&(pvect > -pvmin)) pvect = -pvmin; // particle[j].torque += KTORQUE_BOUNDARY*pvect; - particle[j].torque += KTORQUE_BOUNDARY*pvect*(1.5*r - vabs(distance)); + particle[j].torque += KTORQUE_BOUNDARY*pvect*(SEGMENT_FORCE_EQR*r - vabs(distance)); } } @@ -8883,7 +9380,7 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl /* added 24/9/22 */ else if (angle > segment[i].angle2) angle -= DPI; - r = 1.5*particle[j].radius; + r = SEGMENT_FORCE_EQR*particle[j].radius; if ((distance < r)&&(angle > segment[i].angle1)&&(angle < segment[i].angle2)) { @@ -8917,7 +9414,7 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl case (BC_RECTANGLE): { /* add harmonic force outside rectangular box */ - padding = MU + 0.01; + padding = particle[j].radius + 0.01; xmin = xleft + padding; xmax = xright - padding; ymin = BCYMIN + padding; @@ -9740,6 +10237,7 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) particle[i].thermostat = 1; particle[i].close_to_boundary = 0; particle[i].emean = 0.0; + particle[i].damping = 1.0; particle[i].dirmean = 0.0; particle[i].paired = 0; particle[i].collision = 0; @@ -9751,7 +10249,7 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) // if (y >= YMAX) y -= particle[i].radius; // if (y <= YMIN) y += particle[i].radius; - if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(0.75 + 0.5*((double)rand()/RAND_MAX)); + if (RANDOM_RADIUS) particle[i].radius = particle[i].radius*(RANDOM_RADIUS_MIN + RANDOM_RADIUS_RANGE*((double)rand()/RAND_MAX)); if (particle[i].type == 0) { @@ -9760,6 +10258,10 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) particle[i].spin_range = SPIN_RANGE; particle[i].spin_freq = SPIN_INTER_FREQUENCY; particle[i].mass_inv = 1.0/PARTICLE_MASS; + if ((RANDOM_RADIUS)&&(ADAPT_MASS_TO_RADIUS > 0.0)) + particle[i].mass_inv *= 1.0/(1.0 + pow(particle[i].radius/MU, ADAPT_MASS_TO_RADIUS)); + if ((RANDOM_RADIUS)&&(ADAPT_DAMPING_TO_RADIUS > 0.0)) + particle[i].damping = 1.0 + ADAPT_DAMPING_FACTOR*pow(particle[i].radius/MU, ADAPT_DAMPING_TO_RADIUS); particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; particle[i].charge = CHARGE; } @@ -9851,6 +10353,7 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) particle[i].thermostat = 0; particle[i].energy = 0.0; particle[i].emean = 0.0; + particle[i].damping = 1.0; particle[i].dirmean = 0.0; particle[i].mass_inv = 1.0/PARTICLE_MASS; particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT; @@ -9914,7 +10417,7 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES]) } /* change type of tracer particle */ - if (TRACER_PARTICLE) for (j=0; j 0.0) - add_particle(x, y, 5.0*V_INITIAL*(double)rand()/RAND_MAX, -10.0*V_INITIAL, PARTICLE_MASS, type, particle); - else - add_particle(x, y, 5.0*V_INITIAL*(double)rand()/RAND_MAX, 10.0*V_INITIAL, PARTICLE_MASS, type, particle); +// if (y > 0.0) +// add_particle(x, y, 5.0*V_INITIAL*(double)rand()/RAND_MAX, -10.0*V_INITIAL, PARTICLE_MASS, type, particle); +// else +// add_particle(x, y, 5.0*V_INITIAL*(double)rand()/RAND_MAX, 10.0*V_INITIAL, PARTICLE_MASS, type, particle); particle[ncircles - 1].eq_dist = EQUILIBRIUM_DIST; particle[ncircles - 1].thermostat = 1; px[ncircles - 1] = particle[ncircles - 1].vx; py[ncircles - 1] = particle[ncircles - 1].vy; + if ((TRACER_PARTICLE)&&(n_tracers < N_TRACER_PARTICLES)) + { + tracer_n[n_tracers] = ncircles - 1; + n_tracers++; + printf("%i tracer particles\n", n_tracers); + } + // init_molecule_data(particle, molecule); return (ncircles); @@ -10403,6 +10920,17 @@ int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin, t return(nthermo); } + +int partial_bfield(double x, double y) +/* */ +{ + switch (BFIELD_REGION) { + case (BF_CONST): return(1); + case (BF_SQUARE): return(vabs(x) < YMAX); + } +} + + double compute_mean_energy(t_particle particle[NMAXCIRCLES]) { int i, nactive = 0; @@ -14671,6 +15199,100 @@ void compute_cluster_force(t_cluster cluster[NMAXCIRCLES], t_particle particle[N } } + +void update_conveyor_belts(t_segment segment[NMAXSEGMENTS], t_belt belt[NMAXBELTS]) +/* move shovels of conveyor belt */ +{ + int b, shovel, firstseg, seg; + double position, length, width, angle, newpos, deltapos, beltlength, shift, beta; + + for (b = 0; b < nbelts; b++) + { + length = belt[b].length; + width = belt[b].width; + angle = belt[b].angle; + beltlength = 2.0*length + DPI*width; + for (shovel = 0; shovel < belt[b].nshovels; shovel++) + { + position = belt[b].shovel_pos[shovel]; + deltapos = belt[b].speed*DT_PARTICLE; + newpos = position + deltapos; + firstseg = belt[b].shovel_segment[shovel]; + if (newpos < length) + { + shift = deltapos; + for (seg = firstseg; seg < firstseg+6; seg++) + translate_one_segment(segment, seg, shift*belt[b].tx, shift*belt[b].ty); + } + else if (newpos < length + deltapos) + { + shift = length - position; + beta = (newpos - length)/width; + for (seg = firstseg; seg < firstseg+6; seg++) + { + translate_one_segment(segment, seg, shift*belt[b].tx, shift*belt[b].ty); + rotate_one_segment(segment, seg, -beta, belt[b].x2, belt[b].y2); + } + } + else if (newpos < length + PI*width) + { + beta = deltapos/width; + for (seg = firstseg; seg < firstseg+6; seg++) + { + rotate_one_segment(segment, seg, -beta, belt[b].x2, belt[b].y2); + } + } + else if (newpos < length + PI*width + deltapos) + { + beta = (length + PI*width - position)/width; + shift = newpos - length - PI*width; + for (seg = firstseg; seg < firstseg+6; seg++) + { + rotate_one_segment(segment, seg, -beta, belt[b].x2, belt[b].y2); + translate_one_segment(segment, seg, -shift*belt[b].tx, -shift*belt[b].ty); + } + } + else if (newpos < 2.0*length + PI*width) + { + shift = deltapos; + for (seg = firstseg; seg < firstseg+6; seg++) + translate_one_segment(segment, seg, -shift*belt[b].tx, -shift*belt[b].ty); + } + else if (newpos < 2.0*length + PI*width + deltapos) + { + shift = 2.0*length + PI*width - position; + beta = (newpos - 2.0*length - PI*width)/width; + for (seg = firstseg; seg < firstseg+6; seg++) + { + translate_one_segment(segment, seg, -shift*belt[b].tx, -shift*belt[b].ty); + rotate_one_segment(segment, seg, -beta, belt[b].x2, belt[b].y2); + } + } + else if (newpos < beltlength) + { + beta = deltapos/width; + for (seg = firstseg; seg < firstseg+6; seg++) + { + rotate_one_segment(segment, seg, -beta, belt[b].x1, belt[b].y1); + } + } + else + { + beta = (beltlength - position)/width; + shift = newpos - beltlength; + for (seg = firstseg; seg < firstseg+6; seg++) + { + rotate_one_segment(segment, seg, -beta, belt[b].x1, belt[b].y1); + translate_one_segment(segment, seg, shift*belt[b].tx, shift*belt[b].ty); + } + } + + if (newpos > beltlength) newpos -= beltlength; + belt[b].shovel_pos[shovel] = newpos; + } + } +} + void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_position, int traj_length, int wall, double pressure[N_PRESSURES], double pleft, double pright, @@ -14681,10 +15303,10 @@ void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_positio t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], t_group_data *group_speeds, t_group_segments *segment_group, - t_belt *conveyor_belt) + t_belt *conveyor_belt, int *tracer_n) /* draw a movie frame */ { - if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length); + if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length, particle, cluster, tracer_n, plot); draw_particles(particle, cluster, plot, params.beta, collisions, ncollisions, bg_color, hashgrid, params); draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, conveyor_belt, wall); if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, refresh); diff --git a/sub_rde.c b/sub_rde.c index 0c98fb5..03c54f7 100644 --- a/sub_rde.c +++ b/sub_rde.c @@ -1421,6 +1421,142 @@ double init_linear_blob_sphere(int add, double phi0, double theta0, double vx, d } } for (j=0; j NX) ii -= NX; + /* NEW */ + if (FLOODING) + draw = ((wsphere[i*NY+j].indomain)||((*rde[i*NY+j].p_zfield[movie] >= wsphere[i*NY+j].altitude + FLOODING_VSHIFT))); + else draw = wsphere[i*NY+j].indomain; - glVertex2i(HRES*ii, HRES*j); - glVertex2i(HRES*(ii+1), HRES*j); - glVertex2i(HRES*(ii+1), HRES*(j+1)); - glVertex2i(HRES*ii, HRES*(j+1)); + if (draw) + { + for (k=0; k<3; k++) rgb[k] = rde[i*NY+j].rgb[k]; + glColor3f(rgb[0], rgb[1], rgb[2]); + + ii = NX-i-1+ishift; + if (ii > NX) ii -= NX; + + glVertex2i(HRES*ii, HRES*j); + glVertex2i(HRES*(ii+1), HRES*j); + glVertex2i(HRES*(ii+1), HRES*(j+1)); + glVertex2i(HRES*ii, HRES*(j+1)); + } } glEnd (); @@ -6006,7 +6207,17 @@ void draw_wave_2d_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere ws for (i=0; i= wsphere[i1*NY+j1].altitude + FLOODING_VSHIFT))); + else draw = wsphere_hr[i*HRES*NY+j].indomain; +// if (FLOODING) draw = ((wsphere_hr[i*HRES*NY+j].indomain)||((*rde[i1*NY+j1].p_zfield[movie] >= wsphere[i1*NY+j1].altitude + FLOODING_VSHIFT))); +// else draw = wsphere_hr[i*HRES*NY+j].indomain; +// draw_hr[p*HRES+q] = (draw)&&(wsphere_hr[(HRES*i+p)*HRES*NY+HRES*j+q].indomain); +// if (!wsphere_hr[i*HRES*NY+j].indomain) + if (!draw) { ca = wsphere_hr[i*HRES*NY+j].cos_angle; ca = (ca + 1.0)*0.4 + 0.2; @@ -6028,6 +6239,116 @@ void draw_wave_2d_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere ws if (DRAW_BILLIARD) draw_billiard(0, 1.0); } +void draw_wave_2d_rde_ij(int i, int j, int movie, short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY], t_wave_sphere wsphere_hr[HRES*HRES*NX*NY], int fade, double fade_value) +{ + int p, q, k, ii, i1, j1, iplus1, jplus1; + double ca, rgb[3], r_hr[HRES*HRES], g_hr[HRES*HRES], b_hr[HRES*HRES]; + short int draw, draw_hr[HRES*HRES]; + static int ishift, first = 1; + + if (first) + { + ishift = (int)((double)NX*PHISHIFT/360.0); + first = 0; + } + + if (FLOODING) + draw = ((wsphere[i*NY+j].indomain)||((*rde[i*NY+j].p_zfield[movie] >= wsphere[i*NY+j].altitude + FLOODING_VSHIFT_2D))); + else draw = wsphere[i*NY+j].indomain; + + for (p=0; p NX) ii -= NX; + + if (draw) + { + for (k=0; k<3; k++) rgb[k] = rde[i*NY+j].rgb[k]; + glColor3f(rgb[0], rgb[1], rgb[2]); + + glVertex2i(HRES*ii, HRES*j); + glVertex2i(HRES*(ii+1), HRES*j); + glVertex2i(HRES*(ii+1), HRES*(j+1)); + glVertex2i(HRES*ii, HRES*(j+1)); + } + + /* draw continents in higher resolution */ + if (RDE_PLANET) for (p=0; p NX) ii -= NX; + else if (ii < 0) ii += NX; + glColor3f(fade_value, fade_value, fade_value); + glBegin(GL_LINE_STRIP); + glVertex2i(ii*HRES, 0); + glVertex2i(ii*HRES, NY*HRES); + glEnd(); + } +} + void draw_wave_sphere_rde_ij(int i, int iplus, int j, int jplus, int jcolor, int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY], t_wave_sphere wsphere_hr[HRES*HRES*NX*NY], int zplot, int cplot, int palette, int fade, double fade_value) /* draw wave at simulation grid point (i,j) */ { @@ -6039,10 +6360,13 @@ void draw_wave_sphere_rde_ij(int i, int iplus, int j, int jplus, int jcolor, int draw_bc = (xy_in[i*NY+j])&&(xy_in[iplus*NY+j])&&(xy_in[i*NY+jplus])&&(xy_in[iplus*NY+jplus]); // draw = wsphere[i*NY+j].indomain; - draw = 1; + if (FLOODING) + draw = ((wsphere[i*NY+j].indomain)||((*rde[i*NY+j].p_zfield[movie] >= wsphere[i*NY+j].altitude + FLOODING_VSHIFT))); + else draw = wsphere[i*NY+j].indomain; +// draw = 1; for (p=0; p XMAX) x1 += XMIN - XMAX; y1 = tracers[2*(t-1)*N_TRACERS + 2*tracer + 1]; - x2 = -tracers[2*t*N_TRACERS + 2*tracer]; + x2 = -tracers[2*t*N_TRACERS + 2*tracer] + xshift; + if (x2 > XMAX) x2 += XMIN - XMAX; y2 = tracers[2*t*N_TRACERS + 2*tracer + 1]; lum = 1.0 - 0.75*(double)(time - t)/(double)length; @@ -6740,7 +7075,7 @@ void draw_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int glColor3f(lum, lum, lum); - if (module2(x2 - x1, y2 - y1) < 0.2) draw_line_hres(x1, y1, x2, y2); + if (module2(x2 - x1, y2 - y1) < maxlength) draw_line_hres(x1, y1, x2, y2); // printf("time = %i, tracer = %i, coord = %i, x1 = %.2lg, y1 = %.2lg, x2 = %.2lg, y2 = %.2lg\n", t, tracer,2*t*N_TRACERS + 2*tracer, x1, y1, x2, y2); } diff --git a/sub_sphere.c b/sub_sphere.c index cbb2606..4e5ea51 100644 --- a/sub_sphere.c +++ b/sub_sphere.c @@ -480,6 +480,11 @@ void init_dem(t_wave_sphere wsphere[NX*NY], int dem_number) // else wsphere[i*NY+j].radius = 1.0; wsphere[i*NY+j].radius_dem = 1.0 + RSCALE_DEM*(wsphere[i*NY+j].altitude - vshift); + /* set domain of evolution for simulations with flooding */ + for (i=0; i x2)&&(vabs(y) < WALL_WIDTH)) return(1); + return(0); + } + case (D_POLYCIRCLES): + { + if (module2(x, y) < LAMBDA) return(1); + for (k=0; k 0.0)&&(x1 < LAMBDA + 2.0*MU)&&(vabs(y1) < WALL_WIDTH)) return(1); + } + return(0); + } + case (D_POLYCIRCLES_ANGLED): + { + if (module2(x, y) < LAMBDA) return(1); + alpha = -0.5*PID; + cb = cos(alpha); + sb = sin(alpha); + for (k=0; k 0.5*LAMBDA)&&(x1 < LAMBDA+MU)&&(b < c1)&&(b > c2)) return(1); + } return(0); } case (D_MENGER): @@ -4029,7 +4074,7 @@ void draw_rc_hyp() void draw_billiard(int fade, double fade_value) /* draws the billiard boundary */ { - double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, d, r1, r2, ymax, height, xmin, xmax, ca, sa, xshift, x5, x6, f, fp, xratio, w; + double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, alpha2, dphi, omega, z, l, width, a, b, c, d, r1, r2, ymax, height, xmin, xmax, ca, sa, xshift, x5, x6, f, fp, xratio, w; int i, j, k, k1, k2, mr2, ntiles; static int first = 1, nsides; static double h, hh, sqr3, ll, salpha, arcangle; @@ -5843,6 +5888,57 @@ void draw_billiard(int fade, double fade_value) /* draws the billiard bound /* TODO */ break; } + case (D_TWOCIRCLES): + { + x1 = 0.25*(LAMBDA + 5.0*MU); + x2 = -0.75*(LAMBDA + 2.0*MU); + alpha = asin(WALL_WIDTH/LAMBDA); + alpha2 = asin(WALL_WIDTH/MU); + draw_circle_arc(x1, 0.0, LAMBDA, -PI + alpha, DPI - 2.0*alpha, NSEG); + draw_line(x1 - LAMBDA*cos(alpha), WALL_WIDTH, x2 + MU*cos(alpha2), WALL_WIDTH); + draw_circle_arc(x2, 0.0, MU, alpha2, DPI - 2.0*alpha2, NSEG); + draw_line(x2 + MU*cos(alpha2), -WALL_WIDTH, x1 - LAMBDA*cos(alpha), -WALL_WIDTH); + break; + } + case (D_POLYCIRCLES): + { + alpha = asin(WALL_WIDTH/LAMBDA); + alpha2 = asin(WALL_WIDTH/MU); + for (k=0; k nmaxpixels) + { + printf("bathymetric data file too large, increase nmaxpixels in read_negative_dem_values()\n"); + exit(0); + } + + /* shift due to min/max latitudes of image */ + sshift = 0 + DPOLE; + nshift = 0 + DPOLE; + + /* choice of zero meridian */ + ishift = (int)(nx*ZERO_MERIDIAN/360.0); + + /* read rgb values */ + for (j=0; j hcont)) hmin = rgbtot; + if (rgbtot > hmax) hmax = rgbtot; + int_height_values[3*(j*nx+i)] = rgbtot; + } + + printf("hmin = %i, hmax = %i\n", hmin, hmax); + + /* remove remaining black continents */ + for (i=0; i<3*nx*ny; i++) if (int_height_values[i] < hcont) int_height_values[i] = hmax; + + cratio = 1.0/(double)(hmax-hmin); + rx = (double)nx/(double)NX; + ry = (double)ny/(double)(NY - sshift - nshift); + + /* build underwater height table */ + for (i=0; i nx-1) ii -= nx; + jj = (int)(ry*(double)(NY-nshift - j)); + if (jj > ny-1) jj = ny-1; + if (jj < 0) jj = 0; + + if (wsphere[i*NY+j].indomain) + { + /* set height to zero if color is black (due to black patches in bathymetric map) */ + if (int_height_values[3*(jj*nx+ii)] < hcont) height = 0.0; + else height = -1.0 + (double)(int_height_values[3*(jj*nx+ii)])*cratio; + if (height > 0.0) height = 0.0; + height_values_tmp[i*NY+j] = height; + wsphere[i*NY+j].altitude = height; + +// if (int_height_values[3*(jj*nx+ii)] > hcont) +// { +// wsphere[i*NY+j].r = 0.9*(double)rgb_values[3*(jj*nx+ii)]*cratio; +// wsphere[i*NY+j].g = 0.9*(double)rgb_values[3*(jj*nx+ii)+1]*cratio; +// wsphere[i*NY+j].b = 0.9*(double)rgb_values[3*(jj*nx+ii)+2]*cratio; +// } +// else +// { +// wsphere[i*NY+j].r = 0.29; +// wsphere[i*NY+j].g = 0.29; +// wsphere[i*NY+j].b = 0.29; +// } + } + else + { + height_values_tmp[i*NY+j] = 0.0; + height_values_tmp2[i*NY+j] = 0.0; + } + } + + /* smoothen values at low depth */ + if (SMOOTH_DEM) for (k=1; k= -0.25)) + { + height_values_tmp2[i*NY+j] = height_values_tmp[i*NY+j] + 0.1*(height_values_tmp[(i+1)*NY+j] + height_values_tmp[(i-1)*NY+j] + height_values_tmp[i*NY+j+1] + height_values_tmp[i*NY+j-1] - 4.0*height_values_tmp[i*NY+j]); + + height_values_tmp[i*NY+j] = height_values_tmp2[i*NY+j] + 0.1*(height_values_tmp2[(i+1)*NY+j] + height_values_tmp2[(i-1)*NY+j] + height_values_tmp2[i*NY+j+1] + height_values_tmp2[i*NY+j-1] - 4.0*height_values_tmp2[i*NY+j]); + } + } + + if (SMOOTH_DEM) for (i=1; i= -0.25)) + { + wsphere[i*NY+j].altitude = height_values_tmp[i*NY+j]; + } + + for (i=0; i= 0.0) wsphere[i*rny+j].radius_dem = 1.0 + RSCALE_DEM*altitude; +// else wsphere[i*rny+j].radius_dem = 1.0; +// } + + /* set domain in which wave is evolved */ for (i=0; i