Add files via upload

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 */

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;
 

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 */

lennardjones.c

@@ -34,10 +34,10 @@
 #include <sys/types.h>
 #include <tiffio.h>     /* Sam Leffler's libtiff library. */
 #include <omp.h>
-#include <time.h>
+#include <time.h>   
 
 #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 INITXMIN -2.1
-#define INITXMAX 1.8	/* x interval for initial condition */
+#define INITXMAX -2.0	/* x interval for initial condition */
-#define INITYMIN 0.9
+#define INITYMIN -0.7
-#define INITYMAX 2.2	/* y interval for initial condition */
+#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 ADDXMIN -2.1
-#define ADDXMAX 1.5	/* x interval for adding particles */
+#define ADDXMAX -2.0	/* x interval for adding particles */
-#define ADDYMIN -1.0
+#define ADDYMIN -0.6
-#define ADDYMAX 1.0	/* y interval for adding particles */
+#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_SPEED1 10.0     /* speed of first conveyor belt */
-#define BELT_SPEED2 -10.0   /* speed of second 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 4.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 SPIN_INTER_FREQUENCY_B 4.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 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 NGRIDX 1           /* number of grid point for grid of disks */
-#define NGRIDY 36           /* 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 NSTEPS 4800      /* number of frames of movie */
-#define NVID 200          /* number of iterations between images displayed on screen */
+#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 23
-#define PLOT_B 14        /* plot type for second movie */
+#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 V_INITIAL 800.0        /* initial velocity range */
-#define OMEGA_INITIAL 50.0        /* initial angular 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 NBH_DIST_FACTOR 4.0       /* radius in which to count neighbours */
-#define GRAVITY 4000.0            /* gravity acting on all particles */
+#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 1.0e5         /* force constant in angular dynamics */
-#define KTORQUE_BOUNDARY 5.0e5  /* constant in torque from the boundary */
+#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_TIME 0       /* time at which to add first particle */
-#define ADD_PERIOD 2       /* time interval between adding further particles */
+#define ADD_PERIOD 15       /* time interval between adding further particles */
-#define N_ADD_PARTICLES 8  /* number of particles to add */
+#define N_ADD_PARTICLES 1  /* number of particles to add */
-#define FINAL_NOADD_PERIOD 4700  /* final period where no particles are added */
+#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 TRACER_PARTICLE 1   /* set to 1 to have a tracer particle */
-#define N_TRACER_PARTICLES 3    /* number of tracer particles */
+#define N_TRACER_PARTICLES 500    /* number of tracer particles */
-#define TRAJECTORY_LENGTH 8000   /* length of recorded trajectory */
+#define TRAJECTORY_LENGTH 4800    /* length of recorded trajectory */
-#define TRACER_PARTICLE_MASS 4.0    /* relative mass of tracer particle */
+#define TRACER_LUM_FACTOR 5.0     /* controls luminosity decrease of trajectories with time */
-#define TRAJECTORY_WIDTH 3      /* width of tracer particle trajectory */
+#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 HASHX 40     /* size of hashgrid in x direction */
-#define HASHY 88      /* size of hashgrid in y 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<ncircles; j++) if (particle[j].active) 
     {
+        if (ADAPT_DAMPING_TO_RADIUS > 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);
+            px[j] *= exp(- DT_PARTICLE*damping1);
-            py[j] *= exp(- DT_PARTICLE*damping);
+            py[j] *= exp(- DT_PARTICLE*damping1);
-            pangle[j] *= exp(- DT_PARTICLE*DAMPING_ROT);
+            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;
+                    bfield = params.bfield;
-                    particle[j].fy -= params.bfield*particle[j].charge*particle[j].vx*particle[j].mass_inv;
+                    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<N_TRACER_PARTICLES; j++)
+            for (j=0; j<n_tracers; j++)
             {
                 trajectory[j*TRAJECTORY_LENGTH + traj_position].xc = particle[tracer_n[j]].xc;
                 trajectory[j*TRAJECTORY_LENGTH + traj_position].yc = particle[tracer_n[j]].yc;
@@ -2126,7 +2178,7 @@ void animation()
             change_type_proportion(particle, params.prop);
         }
  
-        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);
 
@@ -2137,13 +2189,17 @@ void animation()
             if ((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_DNA_ENZYME)||(RD_REACTION == CHEM_DNA_ENZYME_REPAIR)))
             {
                 for (k=0; k<N_ADD_PARTICLES; k++)
-                    nadd_particle = add_particles(particle, px, py, nadd_particle, 7, molecule);
+                    nadd_particle = add_particles(particle, px, py, nadd_particle, 7, molecule, tracer_n);
             }
             else for (k=0; k<N_ADD_PARTICLES; k++)
-                nadd_particle = add_particles(particle, px, py, nadd_particle, 0, molecule);
+                nadd_particle = add_particles(particle, px, py, nadd_particle, 0, molecule, tracer_n);
 //             params.nactive = nadd_particle;
             params.nactive = 0;
-            for (j=0; j<ncircles; j++) if (particle[j].active) params.nactive++;
+            for (j=0; j<ncircles; j++) if (particle[j].active) 
+            {
+                params.nactive++;
+                pangle[j] = particle[j].omega;
+            }
         }
         
         /* change particle radius */
@@ -2180,7 +2236,7 @@ void animation()
         
         draw_frame(i, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, 
         wall, pressure, pleft, pright, particle_numbers, 1, params, particle, cluster, 
-        collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt);
+        collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n);
         
 	if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
 
@@ -2212,7 +2268,7 @@ void animation()
             {
                 draw_frame(i, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, 
                 wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, 
-                collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt);
+                collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n);
                 glutSwapBuffers();
                 save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter);
                 counter++;
@@ -2253,7 +2309,7 @@ void animation()
             blank();
             draw_frame(NSTEPS, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, 
             wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, 
-            collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt);
+            collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n);
         }
         if (DOUBLE_MOVIE) for (i=0; i<MID_FRAMES; i++) 
         {
@@ -2265,7 +2321,7 @@ void animation()
         {
             draw_frame(NSTEPS, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, 
             wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, 
-            collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt);
+            collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n);
             if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
         }
         if ((TIME_LAPSE)&&(!DOUBLE_MOVIE)) 

rde.c

@@ -40,7 +40,7 @@
 #include <time.h>
 
 #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 G_FIELD 0.025    /* gravity/Coriolis force */
-#define BC_FIELD 1.0e-5     /* constant in repulsive field from obstacles */
+#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 SWATER_MIN_HEIGHT 0.05     /* min height of initial condition for shallow water equation */
-#define DEPTH_FACTOR 0.075       /* proportion of min height in variable depth */
+#define DEPTH_FACTOR 0.00002       /* proportion of min height in variable depth */
-#define TANH_FACTOR 1.0         /* steepness of 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 NSTEPS 1400       /* number of frames of movie */
-#define NVID 85           /* number of iterations between images displayed on screen */
+#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 11       /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 12      /* 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 VSCALE_WATER_HEIGHT 15.0     /* vertical scaling of water height */
-#define ADD_HEIGHT_CONSTANT -0.025   /* constant added to wave 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 light[3] = {-0.816496581, 0.40824829, 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 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 ADD_NEGATIVE_DEM 1      /* add DEM with bathymetric data */
-#define RSCALE_DEM 0.1          /* scaling factor of radial component for DEM */
+#define RSCALE_DEM 0.075        /* scaling factor of radial component for DEM */
-#define SMOOTH_DEM 0            /* set to 1 to smoothen DEM (to make altitude less constant) */
+#define SMOOTH_DEM 1            /* set to 1 to smoothen DEM (to make altitude less constant) */
-#define DEM_SMOOTH_STEPS 1      /* number of smoothening steps */
+#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 RSCALE 0.15              /* scaling factor of radial component */
-#define RSHIFT -0.01             /* shift in 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<NX-1; i++){
-                for (j=1; j<NY-1; j++){
+                for (j=1; j<NY-1; j++)
-                    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;
+                    if (wsphere[i*NY+j].altitude > 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; i++){
         for (j=0; j<NY; j++){
@@ -1157,6 +1197,30 @@ double initialize_water_depth(t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
     return(max);
 }
 
+void compute_forcing_schedule(int t, t_wave_sphere wsphere[NX*NY])
+/* compute periodic forcing */
+{
+    int i, j;
+    double phase;
+    
+    phase = DPI*(double)t/(double)FORCING_PERIOD + ZERO_MERIDIAN*PI/180.0 + FORCING_SHIFT*PID;
+    
+    for (i=0; i<NX; i++) for (j=0.0; j<NY; j++)
+    {
+        wsphere[i*NY+j].force = FORCING_AMP*wsphere[i*NY+j].stheta*sin(2.0*(wsphere[i*NY+j].phi - phase));
+        wsphere[i*NY+j].force += FORCING_CONST_AMP*wsphere[i*NY+j].stheta*cos(2.0*(wsphere[i*NY+j].phi));
+    }
+    
+    moon_position = (int)((double)NX*(phase/DPI + 0.625));
+    while (moon_position < 0) moon_position += NX;
+    while (moon_position >= 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*wsphere[i*NY+j].force;
-                        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;
+                    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<NX*NY; cell++) 
@@ -1621,6 +1724,8 @@ void update_tracer_table(double tracers[2*N_TRACERS*NSTEPS], t_rde rde[NX*NY], i
     #pragma omp parallel for private(tracer)
     for (tracer = 0; tracer < N_TRACERS; tracer++)
     {
+        x1 = XMIN;
+        y1 = YMIN;
         for (t = 0; t < maxtime; t++)
         {
             t1 = time - t;
@@ -1641,10 +1746,13 @@ void update_tracer_table(double tracers[2*N_TRACERS*NSTEPS], t_rde rde[NX*NY], i
             }
 //             else printf("More than %i tracer points per cell\n", NMAX_TRACER_PTS);
             
-            if ((cell != oldcell)&&(t > 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<NFIELDS; k++) block_poles(phi[k]);
@@ -2141,6 +2264,8 @@ void animation()
         tracers[2*i+1] = YMIN + 0.05 + (YMAX - YMIN - 0.1)*rand()/RAND_MAX;
     }
     
+    if (ADD_MOON_FORCING) compute_forcing_schedule(0, wsphere);
+    
     blank();
     glColor3f(0.0, 0.0, 0.0);
     
@@ -2179,6 +2304,7 @@ void animation()
         if (CHANGE_RPSLZB) rpslzb = rpslzb_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);
         
         if (ROTATE_VIEW) 
         {

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<DPOLE; j++)
+        {
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*sin(theta0);
+            pscal += wsphere[i*NY+j].ctheta*cos(theta0);
+            dist2 = acos(pscal);
+            dist2 *= dist2;
+            module = amp*exp(-dist2/varx);
+            
+            phi[0][i*NY+j] = min + a*module*wsphere[i*NY+DPOLE].stheta;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+        for (j=NY-DPOLE; j<NY; j++)
+        {
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*sin(theta0);
+            pscal += wsphere[i*NY+j].ctheta*cos(theta0);
+            dist2 = acos(pscal);
+            dist2 *= dist2;
+            module = amp*exp(-dist2/varx);
+            
+            phi[0][i*NY+j] = min + a*module*wsphere[i*NY+NY-DPOLE].stheta;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+    }
+}
+
+void init_expanding_blob_sphere(int add, double phi0, double theta0, double v, double amp, double radius, double min, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
+/* initialise gaussian wave height for shallow water equation */
+{
+    int i, j;
+    double xy[2], var, a, pscal, dist2, module, x0, y0, z0, x1, y1, z1, x2, y2, z2, r, vphi, vtheta;
+    
+    var = radius*radius; 
+    a = amp/(sqrt(DPI*var));
+    
+    printf("[init_gaussian_wave_sphere]: a = %.3lg, var = %.3lg\n", a, var); 
+    
+    if (add) for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*sin(theta0);
+            pscal += wsphere[i*NY+j].ctheta*cos(theta0);
+            dist2 = acos(pscal);
+            dist2 *= dist2;
+            module = amp*exp(-dist2/var);
+            
+            /* compute radial vector between X0 and X */
+            x0 = cos(phi0)*sin(theta0);
+            y0 = sin(phi0)*sin(theta0);
+            z0 = -cos(theta0);
+            
+            x1 = y0*wsphere[i*NY+j].z - z0*wsphere[i*NY+j].y;
+            y1 = z0*wsphere[i*NY+j].x - x0*wsphere[i*NY+j].z;
+            z1 = x0*wsphere[i*NY+j].y - y0*wsphere[i*NY+j].x;
+            
+            x2 = y1*wsphere[i*NY+j].z - z1*wsphere[i*NY+j].y;
+            y2 = z1*wsphere[i*NY+j].x - x1*wsphere[i*NY+j].z;
+            z2 = x1*wsphere[i*NY+j].y - y1*wsphere[i*NY+j].x;
+            
+            r = x2*x2 + y2*y2 + z2*z2;
+            r = 1.0/sqrt(r);
+            x2 *= r;
+            y2 *= r;
+            z2 *= r;
+            
+            vphi = -x2*sin(wsphere[i*NY+j].phi) + y2*cos(wsphere[i*NY+j].phi);
+            vtheta = (x2*cos(wsphere[i*NY+j].phi) + y2*sin(wsphere[i*NY+j].phi))*cos(wsphere[i*NY+j].theta) + z2*sin(wsphere[i*NY+j].theta);
+            
+            
+//             printf("[init_gaussian_wave_sphere]: dist2 = %.3lg, module = %.3lg\n", dist2, module);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] += a*module*wsphere[i*NY+j].stheta;
+                phi[1][i*NY+j] += v*module*vphi;
+                phi[2][i*NY+j] += v*module*vtheta;
+            }
+        }
+    }
+    else for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*sin(theta0);
+            pscal += wsphere[i*NY+j].ctheta*cos(theta0);
+            dist2 = acos(pscal);
+            dist2 *= dist2;
+            module = amp*exp(-dist2/var);
+            
+            /* compute radial vector between X0 and X */
+            x0 = cos(phi0)*sin(theta0);
+            y0 = sin(phi0)*sin(theta0);
+            z0 = -cos(theta0);
+            
+            x1 = y0*wsphere[i*NY+j].z - z0*wsphere[i*NY+j].y;
+            y1 = z0*wsphere[i*NY+j].x - x0*wsphere[i*NY+j].z;
+            z1 = x0*wsphere[i*NY+j].y - y0*wsphere[i*NY+j].x;
+            
+            x2 = y1*wsphere[i*NY+j].z - z1*wsphere[i*NY+j].y;
+            y2 = z1*wsphere[i*NY+j].x - x1*wsphere[i*NY+j].z;
+            z2 = x1*wsphere[i*NY+j].y - y1*wsphere[i*NY+j].x;
+            
+            r = x2*x2 + y2*y2 + z2*z2;
+            r = 1.0/sqrt(r);
+            x2 *= r;
+            y2 *= r;
+            z2 *= r;
+            
+            vphi = -x2*sin(wsphere[i*NY+j].phi) + y2*cos(wsphere[i*NY+j].phi);
+            vtheta = (x2*cos(wsphere[i*NY+j].phi) + y2*sin(wsphere[i*NY+j].phi))*cos(wsphere[i*NY+j].theta) + z2*sin(wsphere[i*NY+j].theta);
+            
+//             printf("[init_gaussian_wave_sphere]: dist2 = %.3lg, module = %.3lg\n", dist2, module);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] = min + a*module*wsphere[i*NY+j].stheta;
+                phi[1][i*NY+j] = v*module*vphi;
+                phi[2][i*NY+j] = v*module*vtheta;
+//                 phi[1][i*NY+j] = v*module*sin(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta;
+//                 phi[2][i*NY+j] = v*module*sin(wsphere[i*NY+j].theta - theta0)*wsphere[i*NY+j].stheta;
+            }
+            else
+            {
+                phi[0][i*NY+j] = min;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+        }
+        for (j=0; j<DPOLE; j++)
         {
             phi[0][i*NY+j] = min;
             phi[1][i*NY+j] = 0.0;
@@ -1617,6 +1753,63 @@ double init_swater_laminar_flow_sphere(int add, double vx, double amp, int nx, i
     }
 }
 
+void init_tidal_state(int add, double amp, double min, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
+/* initial state for tide simulations */
+{
+    int i, j;
+    double phase, xy[2];
+    
+    phase = ZERO_MERIDIAN*PI/180.0 + 0.5*PID;
+        
+    if (add) for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] += amp*wsphere[i*NY+j].stheta*sin(2.0*(wsphere[i*NY+j].phi - phase));
+            }
+        }
+    }
+    else for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] = min + amp*wsphere[i*NY+j].stheta*sin(2.0*(wsphere[i*NY+j].phi - phase));;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+            else
+            {
+                phi[0][i*NY+j] = min;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+        }
+        for (j=0; j<DPOLE; j++)
+        {
+            phi[0][i*NY+j] = min;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+        for (j=NY-DPOLE; j<NY; j++)
+        {
+            phi[0][i*NY+j] = min;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+    }
+    
+}
+
 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) */
 {
@@ -1909,7 +2102,7 @@ void initialize_bcfield(double bc_field[NX*NY], double bc_field2[NX*NY], t_recta
                 {
                     /* set dummy values */
                     height = wsphere[i*NY+j+1].altitude;
-                    if (height == 0.0) bc_field[i*NY+j] = 1.0;
+                    if (height <= 0.0) bc_field[i*NY+j] = 1.0;
                     else 
                     {   
                         f = tanh(BC_STIFFNESS*height);
@@ -5965,9 +6158,9 @@ void compute_cfield_rde(short int xy_in[NX*NY], int cplot, int palette, t_rde rd
 
 
 
-void draw_wave_2d_rde(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)
+void draw_wave_2d_rde_old(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 movie, int fade, double fade_value)
 {
-    int i, j, k, ii;
+    int i, j, k, ii, i1, j1;
     short int draw;
     double ca, rgb[3];
     static int ishift, first = 1;
@@ -5986,16 +6179,24 @@ void draw_wave_2d_rde(short int xy_in[NX*NY], t_rde rde[NX*NY], t_wave_sphere ws
     for (i=0; i<NX; i++)
         for (j=0; j<NY; j++)
         {
-            for (k=0; k<3; k++) rgb[k] = rde[i*NY+j].rgb[k];
+            /* NEW */
-            glColor3f(rgb[0], rgb[1], rgb[2]);
+            if (FLOODING) 
-                                            
+                draw = ((wsphere[i*NY+j].indomain)||((*rde[i*NY+j].p_zfield[movie] >= wsphere[i*NY+j].altitude + FLOODING_VSHIFT)));
-            ii = NX-i-1+ishift;
+            else draw = wsphere[i*NY+j].indomain;
-            if (ii > NX) ii -= NX;
             
-            glVertex2i(HRES*ii, HRES*j);
+            if (draw)
-            glVertex2i(HRES*(ii+1), HRES*j);
+            {
-            glVertex2i(HRES*(ii+1), HRES*(j+1));
+                for (k=0; k<3; k++) rgb[k] = rde[i*NY+j].rgb[k];
-            glVertex2i(HRES*ii, HRES*(j+1));
+                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<HRES*NX; i++)
             for (j=0; j<HRES*NY; j++)
             {
-                if (!wsphere_hr[i*HRES*NY+j].indomain)
+                i1 = i/HRES;
+                j1 = j/HRES;
+//                 printf("i = %i, j = %j, i1 = %i, j1 = %i\n", i, j, i1, j1);
+                /* NEW */
+                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;
+//                 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<HRES; p++)
+        for (q=0; q<HRES; q++)
+            draw_hr[p*HRES+q] = (draw)||(wsphere_hr[(HRES*i+p)*HRES*NY+HRES*j+q].indomain);
+    
+    for (p=0; p<HRES; p++)
+        for (q=0; q<HRES; q++)
+        {
+            i1 = HRES*i + p;
+            j1 = HRES*j + q;
+            ca = wsphere_hr[i1*HRES*NY+j1].cos_angle;
+            ca = (ca + 1.0)*0.4 + 0.2;
+            if (fade) ca *= fade_value;
+            if (RDE_PLANET)
+            {
+                r_hr[p*HRES + q] = wsphere_hr[i1*HRES*NY+j1].r*ca;
+                g_hr[p*HRES + q] = wsphere_hr[i1*HRES*NY+j1].g*ca;
+                b_hr[p*HRES + q] = wsphere_hr[i1*HRES*NY+j1].b*ca;
+            }
+            else
+            {
+                r_hr[p*HRES + q] = COLOR_OUT_R*ca;
+                g_hr[p*HRES + q] = COLOR_OUT_G*ca;
+                b_hr[p*HRES + q] = COLOR_OUT_B*ca;
+            }
+        }
+    
+    ii = NX-i-1+ishift;
+    if (ii > 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<HRES; p++)
+    {
+        for (q=0; q<HRES; q++) if (!draw_hr[p*HRES+q])
+        {
+            i1 = HRES*ii + p;
+            j1 = HRES*j + q;
+            iplus1 = HRES*(ii+1) + p;
+            jplus1 = HRES*(j+1) + q;
+                
+            glColor3f(r_hr[p*HRES + q], g_hr[p*HRES + q], b_hr[p*HRES + q]);
+            
+            glVertex2i(i1, j1);
+            glVertex2i(iplus1, j1);
+            glVertex2i(iplus1, jplus1);
+            glVertex2i(i1, jplus1);
+        }
+    } 
+}
+
+void draw_wave_2d_rde(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 movie, int fade, double fade_value)
+{
+    int i, j, ii;
+    static int ishift, first;
+    
+    if (first)
+    {
+        ishift = (int)((double)NX*PHISHIFT/360.0);
+        first = 0;
+    }
+    
+    /* draw the field */
+    glBegin(GL_QUADS);
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+            draw_wave_2d_rde_ij(i, j, movie, xy_in, rde, wsphere, wsphere_hr, fade, fade_value);
+    glEnd ();
+    
+    if (DRAW_BILLIARD) draw_billiard(0, 1.0);
+    
+    if (DRAW_MOON_POSITION)
+    {
+        ii = NX-moon_position-1+ishift;
+        if (ii > 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<HRES; p++)
         for (q=0; q<HRES; q++)
-            draw_hr[p*HRES+q] = wsphere_hr[(HRES*i+p)*HRES*NY+HRES*j+q].indomain;
+            draw_hr[p*HRES+q] = (draw)&&(wsphere_hr[(HRES*i+p)*HRES*NY+HRES*j+q].indomain);
     for (k=0; k<3; k++) rgb[k] = rde[i*NY+jcolor].rgb[k];
         glColor3f(rgb[0], rgb[1], rgb[2]);
     
@@ -6070,7 +6394,7 @@ void draw_wave_sphere_rde_ij(int i, int iplus, int j, int jplus, int jcolor, int
 
     if (draw_bc)
     {
-//         if (draw)
+        if (draw)
         {
             glBegin(GL_TRIANGLE_FAN);
             glColor3f(rgb[0], rgb[1], rgb[2]);
@@ -6709,7 +7033,7 @@ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rd
             draw_periodicised_wave_3d(movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
         else draw_wave_3d_rde(movie, phi, xy_in, rde, potential, zplot, cplot, palette, fade, fade_value);
     }
-    else draw_wave_2d_rde(xy_in, rde, wsphere, wsphere_hr, fade, fade_value);
+    else draw_wave_2d_rde(xy_in, rde, wsphere, wsphere_hr, movie, fade, fade_value);
 }
 
 void draw_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int time, int fade, double fade_value)
@@ -6717,6 +7041,15 @@ void draw_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int
 {
     int tracer, t, t1, length = 50, ij[2];
     double x1, y1, x2, y2, lum;
+    static double xshift, maxlength;
+    static int first = 1;
+    
+    if (first)
+    {
+        xshift = PHISHIFT*(XMAX - XMIN)/360.0;
+        maxlength = 0.01*(XMAX - XMIN);
+        first = 0;
+    }
     
     glColor3f(1.0, 1.0, 1.0);
     glLineWidth(1);
@@ -6729,10 +7062,12 @@ void draw_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], int
     for (t = t1 + 1; t < time; t++) 
         for (tracer = 0; tracer < N_TRACERS; tracer++)
         {
-            x1 = -tracers[2*(t-1)*N_TRACERS + 2*tracer];
+            x1 = -tracers[2*(t-1)*N_TRACERS + 2*tracer] + xshift;
+            if (x1 > 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);
         }

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<NX; i++)
+        for (j=0; j<NY; j++)
+            wsphere[i*NY+j].evolve_wave = (wsphere[i*NY+j].altitude < vshift + 0.1); 
+    
     /* compute light angle */  
     dx = 2.0*(XMAX - XMIN)/(double)NX;
     dy = 2.0*(YMAX - YMIN)/(double)NY;

sub_wave.c

@@ -2700,7 +2700,7 @@ int init_xyin_from_image(short int * xy_in[NX])
 int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_circle *circles)
 /* returns 1 if (x,y) represents a point in the billiard */
 {
-    double l2, r1, r2, r2mu, omega, b, c, d, angle, z, x1, y1, x2, y2, y3, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height, ca, sa, l, ht, xshift, zz[9][2], x5, x6, f, fp, h1;
+    double l2, r1, r2, r2mu, omega, b, c, d, angle, z, x1, y1, x2, y2, y3, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height, ca, sa, l, ht, xshift, zz[9][2], x5, x6, f, fp, h1, cb, sb, c1, c2;
     int i, j, k, k1, k2, condition = 0, m;
     static int first = 1, nsides;
     static double h, hh, ra, rb, ll, salpha;
@@ -3720,8 +3720,6 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
                 r = module2(x - x1, y - y1);
                 if (r < MU) return(1);
             }
-            
-            
             return(0);
         }
         case (D_MAGNETRON_CATHODE):
@@ -3741,8 +3739,55 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
                 r = module2(x - x1, y - y1);
                 if (r < MU) return(1);
             }
-            
+            return(0);
-            
+        }
+        case (D_TWOCIRCLES):
+        {
+            x1 = 0.25*(LAMBDA + 5.0*MU);
+            x2 = -0.75*(LAMBDA + 2.0*MU);
+            if (module2(x - x1, y) < LAMBDA) return(1);
+            if (module2(x - x2, y) < MU) return(1);
+            if ((x < x1)&&(x > x2)&&(vabs(y) < WALL_WIDTH)) return(1);
+            return(0);
+        }
+        case (D_POLYCIRCLES):
+        {
+            if (module2(x, y) < LAMBDA) return(1);
+            for (k=0; k<NPOLY; k++)
+            {
+                angle = APOLY*PID + (double)k*DPI/(double)NPOLY;
+                ca = cos(angle);
+                sa = sin(angle);
+                x1 = (LAMBDA + 2.0*MU)*ca;
+                y1 = (LAMBDA + 2.0*MU)*sa;
+                if (module2(x - x1, y - y1) < MU) return(1);
+                x1 = x*ca + y*sa;
+                y1 = -x*sa + y*ca;
+                if ((x1 > 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<NPOLY; k++)
+            {
+                angle = APOLY*PID + (double)k*DPI/(double)NPOLY;
+                ca = cos(angle);
+                sa = sin(angle);
+                x1 = x*ca + y*sa;
+                y1 = -x*sa + y*ca;
+                x2 = LAMBDA + 2.0*MU*cb;
+                y2 = 2.0*MU*sb;
+                if (module2(x1 - x2, y1 - y2) < MU) return(1);
+                c1 = -sb*LAMBDA + cb*WALL_WIDTH;
+                c2 = -sb*LAMBDA - cb*WALL_WIDTH;
+                b = -sb*x1 + cb*y1;
+                if ((x1 > 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<NPOLY; k++)
+            {
+                phi = APOLY*PID + (double)k*DPI/(double)NPOLY;
+                x1 = (LAMBDA + 2.0*MU)*cos(phi);
+                y1 = (LAMBDA + 2.0*MU)*sin(phi);
+                draw_line(LAMBDA*cos(-alpha + phi), LAMBDA*sin(-alpha + phi), x1 - MU*cos(alpha2 + phi), y1 - MU*sin(alpha2 + phi));
+                draw_circle_arc(x1, y1, MU, -PI + alpha2 + phi, DPI - 2.0*alpha2, NSEG); 
+                draw_line(x1 - MU*cos(-alpha2 + phi), y1 - MU*sin(-alpha2 + phi), LAMBDA*cos(alpha + phi), LAMBDA*sin(alpha + phi)); 
+                draw_circle_arc(0.0, 0.0, LAMBDA, alpha + phi, DPI/(double)NPOLY - 2.0*alpha, NSEG); 
+            }
+            break;
+        }
+        case (D_POLYCIRCLES_ANGLED):
+        {
+            alpha = atan(WALL_WIDTH/LAMBDA);
+            alpha2 = asin(WALL_WIDTH/(MU*sqrt(2.0)));
+            for (k=0; k<NPOLY; k++)
+            {
+                phi = APOLY*PID + (double)k*DPI/(double)NPOLY;
+                x1 = LAMBDA*cos(phi);
+                y1 = LAMBDA*sin(phi);
+                x2 = x1 + 2.0*MU*cos(phi - 0.5*PID);
+                y2 = y1 + 2.0*MU*sin(phi - 0.5*PID);
+                draw_line(LAMBDA*cos(phi-alpha), LAMBDA*sin(phi-alpha), 
+                          x2 - MU*cos(phi -0.5*PID + alpha2), 
+                          y2 - MU*sin(phi -0.5*PID + alpha2));
+                draw_circle_arc(x2, y2, MU, PI + phi - 0.5*PID + alpha2, DPI - 2.0*alpha2, NSEG); 
+                draw_line(x2 - MU*cos(phi -0.5*PID - alpha2), 
+                          y2 - MU*sin(phi -0.5*PID - alpha2), 
+                          LAMBDA*cos(phi+alpha), LAMBDA*sin(phi+alpha)); 
+                draw_circle_arc(0.0, 0.0, LAMBDA, phi + alpha, DPI/(double)NPOLY - 2.0*alpha, NSEG);
+            }
+            
+            break;
+        }
         case (D_MENGER):
         {
             glLineWidth(3);

sub_wave_3d_rde.c

@@ -70,6 +70,9 @@ void init_wave_sphere_rde(t_wave_sphere *wsphere, int res)
             wsphere[i*res*NY+j].y2d = xy[1]; 
             
             wsphere[i*res*NY+j].cos_angle_sphere = wsphere[i*res*NY+j].x*light[0] + wsphere[i*res*NY+j].y*light[1] + wsphere[i*res*NY+j].z*light[2];
+            
+            /* default value, to be changed by init_dem */
+            wsphere[i*res*NY+j].evolve_wave = 1;
         }
         
         /* cotangent, taking care of not dividing by zero */
@@ -94,12 +97,152 @@ void init_wave_sphere_rde(t_wave_sphere *wsphere, int res)
     }
 }
 
+
+void read_negative_dem_values_rde(double *height_values, t_wave_sphere *wsphere)
+/* init bathymetric data */
+{
+    int i, j, k, ii, jj, nx, ny, maxrgb, nmaxpixels = 6480000, hmin, hmax, ishift, nshift, sshift, rgbval, scan, rgbtot;
+    int *rgb_values, *int_height_values;
+    int hcont = 50, rgbdiff;
+    double cratio, rx, ry, height;
+    double *height_values_tmp, *height_values_tmp2;
+    FILE *image_file;
+    
+    printf("Reading bathymetric data\n");
+    
+    rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int));
+    int_height_values = (int *)malloc(3*nmaxpixels*sizeof(int));
+    height_values_tmp = (double *)malloc(NX*NY*sizeof(double));
+    height_values_tmp2 = (double *)malloc(NX*NY*sizeof(double));
+    
+//     image_file = fopen("bathymetry_gebco_3600x1800_color.ppm", "r");
+    image_file = fopen("bathymetry_gebco_2560_1280_mod2_color.ppm", "r");
+    scan = fscanf(image_file,"%i %i\n", &nx, &ny);
+    scan = fscanf(image_file,"%i\n", &maxrgb);
+    
+    for (i=0; i<NX*NY; i++) 
+    {
+        height_values_tmp[i] = 0.0;
+        height_values_tmp2[i] = 0.0;
+    }
+    
+    hmin = maxrgb;
+    hmax = 0;
+    
+    if (nx*ny > 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<ny; j++)
+        for (i=0; i<nx; i++)
+        {
+            rgbtot = 0;
+            for (k=0; k<3; k++)
+            {
+                scan = fscanf(image_file,"%i\n", &rgbval);
+                rgb_values[3*(j*nx+i)+k] = rgbval;
+                rgbtot += rgbval;
+            }
+            if ((rgbtot < hmin)&&(rgbtot > 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; i++)
+        for (j=0; j<NY; j++) 
+        {
+            ii = (int)(rx*(double)(NX-1 - i)) + nx/2 + ishift;
+            if (ii > 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<DEM_SMOOTH_STEPS; k++)
+    {
+        printf("Smoothing step %i\n", k);
+        for (i=1; i<NX-1; i++)
+            for (j=1; j<NY-1; j++) 
+                if ((wsphere[i*NY+j].indomain)&&(height_values[i*NY+j] >= -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<NX-1; i++)
+        for (j=1; j<NY-1; j++) 
+            if ((wsphere[i*NY+j].indomain)&&(height_values[i*NY+j] >= -0.25))
+            {
+                wsphere[i*NY+j].altitude = height_values_tmp[i*NY+j];
+            }
+    
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++) 
+            height_values[i*NY+j] = height_values_tmp[i*NY+j];
+    
+    free(rgb_values);
+    free(int_height_values);
+    free(height_values_tmp);
+    free(height_values_tmp2);
+}
+
 void init_dem_rde(t_wave_sphere *wsphere, int dem_number, int res)
 /* init heights from digital elevation map */
 {
     int i, j, ii, jj, k, nx, ny, maxrgb, nmaxpixels = 4915200, scan, rgbval, diff, sshift, nshift, hmin, hmax, ishift, hsum, rnx, rny;
     int *rgb_values;
-    double cratio, rx, ry, cy, dx, dy, pscal, norm, vscale1, vscale2, gradx, grady, deltar, deltai[3], deltaj[3], dphi, dtheta, n[3], hsea, hmean, vshift;
+    double cratio, rx, ry, cy, dx, dy, pscal, norm, vscale1, vscale2, gradx, grady, deltar, deltai[3], deltaj[3], dphi, dtheta, n[3], hsea, hmean, vshift, altitude;
     double *height_values, *height_values_tmp;
     FILE *image_file;
     
@@ -328,13 +471,23 @@ void init_dem_rde(t_wave_sphere *wsphere, int dem_number, int res)
                 wsphere[i*rny+j].altitude = height_values[i*rny+j];
             }
                 
-//     if ((ADD_NEGATIVE_DEM)&&(dem_number == DEM_EARTH)) 
+    if ((ADD_NEGATIVE_DEM)&&(dem_number == DEM_EARTH)) 
-//         read_negative_dem_values(height_values, wsphere);
+        read_negative_dem_values_rde(height_values, wsphere);
             
     /* set radius */
+//     for (i=0; i<rnx; i++)
+//         for (j=0; j<rny; j++)
+//         {
+//             altitude = wsphere[i*rny+j].altitude - vshift;
+//             if (!wsphere[i*rny+j].indomain) wsphere[i*rny+j].radius_dem = 1.0 + RSCALE_DEM*altitude;
+// //             if (altitude >= 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<rnx; i++)
         for (j=0; j<rny; j++)
-                wsphere[i*rny+j].radius_dem = 1.0 + RSCALE_DEM*(wsphere[i*rny+j].altitude - vshift);
+            wsphere[i*rny+j].evolve_wave = (wsphere[i*rny+j].altitude < vshift + 0.01); 
     
     /* compute light angle */  
     dx = 2.0*(XMAX - XMIN)/(double)rnx;
@@ -599,7 +752,10 @@ void init_earth_map_rde(t_wave_sphere *wsphere, int res)
     for (i=0; i<res*NX; i++)
         for (j=0; j<res*NY; j++)
         {
-            wsphere[i*res*NY+j].radius_dem = 1.0 + RSCALE_DEM*(wsphere[i*res*NY+j].altitude - vshift);
+//             wsphere[i*res*NY+j].indomain = (wsphere[i*res*NY+j].altitude < vshift + 1.0e-7);
+            if (!wsphere[i*res*NY+j].indomain)
+                wsphere[i*res*NY+j].radius_dem = 1.0 + RSCALE_DEM*(wsphere[i*res*NY+j].altitude - vshift);
+            else wsphere[i*res*NY+j].radius_dem = 1.0;
 //             printf("Radius_dem at (%i,%i) = %.3lg\n", i, j, wsphere[i*NY+j].radius_dem);
         }
 }

wave_billiard.c

@@ -62,8 +62,8 @@
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
-#define YMIN -1.197916667
+#define YMIN -0.997916667
-#define YMAX 1.197916667	/* y interval for 9/16 aspect ratio */
+#define YMAX 1.397916667	/* y interval for 9/16 aspect ratio */
 
 #define HIGHRES 1       /* set to 1 if resolution of grid is double that of displayed image */
 #define HRES 1          /* dummy, only used by rde.c */
@@ -72,7 +72,7 @@
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 6        /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 792        /* 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 */
 
@@ -85,11 +85,11 @@
 #define PDISC_FACTOR 3.5    /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 0.3	    /* parameter controlling the dimensions of domain */
+#define LAMBDA 0.85	    /* parameter controlling the dimensions of domain */
-#define MU 1.0              /* parameter controlling the dimensions of domain */
+#define MU 0.18             /* parameter controlling the dimensions of domain */
 #define MU_B 1.0            /* parameter controlling the dimensions of domain */
-#define NPOLY 6             /* number of sides of polygon */
+#define NPOLY 3            /* number of sides of polygon */
-#define APOLY -0.666666666666          /* angle by which to turn polygon, in units of Pi/2 */ 
+#define APOLY -0.17           /* 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 */
@@ -97,7 +97,7 @@
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
 #define NGRIDX 60           /* number of grid point for grid of disks */
 #define NGRIDY 25           /* number of grid point for grid of disks */
-#define WALL_WIDTH 0.05      /* width of wall separating lenses */
+#define WALL_WIDTH 0.05     /* width of wall separating lenses */
 #define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 
 #define X_SHOOTER -0.2
@@ -116,7 +116,7 @@
 
 /* Physical parameters of wave equation */
 
-#define TWOSPEEDS 1         /* set to 1 to replace hardcore boundary by medium with different speed */
+#define TWOSPEEDS 0         /* set to 1 to replace hardcore boundary by medium with different speed */
 #define OSCILLATE_LEFT 0     /* set to 1 to add oscilating boundary condition on the left */
 #define OSCILLATE_TOPBOT 0   /* set to 1 to enforce a planar wave on top and bottom boundary */
 #define OSCILLATION_SCHEDULE 61  /* oscillation schedule, see list in global_pdes.c */
@@ -146,9 +146,10 @@
 /* For similar wave forms, COURANT^2*GAMMA should be kept constant */
 
 #define ADD_OSCILLATING_SOURCE 1        /* set to 1 to add an oscillating wave source */
-#define OSCILLATING_SOURCE_PERIOD 3     /* period of oscillating source */
+#define OSCILLATING_SOURCE_PERIOD 16    /* 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 3                     /* number of sources, for option draw_sources */
+#define ALTERNATE_SOURCE_PHASES 0       /* set to 1 to alternate initial phases of sources */
 
 #define ADD_WAVE_PACKET_SOURCES 0       /* set to 1 to add several sources emitting wave packets */
 #define WAVE_PACKET_SOURCE_TYPE 3       /* type of wave packet sources */
@@ -163,9 +164,8 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 1000      /* number of frames of movie */
+#define NSTEPS 2300         /* number of frames of movie */
-// #define NSTEPS 2000         /* number of frames of movie */
+#define NVID 8              /* number of iterations between images displayed on screen */
-#define NVID 10             /* number of iterations between images displayed on screen */
 #define NSEG 1000           /* number of segments of boundary */
 #define INITIAL_TIME 0      /* time after which to start saving frames */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
@@ -182,7 +182,7 @@
 
 /* Parameters of initial condition */
 
-#define INITIAL_AMP 0.75             /* amplitude of initial condition */
+#define INITIAL_AMP 1.5            /* amplitude of initial condition */
 #define INITIAL_VARIANCE 0.00001    /* variance of initial condition */
 #define INITIAL_WAVELENGTH  0.025   /* wavelength of initial condition */
 
@@ -194,8 +194,8 @@
 
 /* Color schemes */
 
-#define COLOR_PALETTE 11     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 13     /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 13   /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 17   /* Color palette, see list in global_pdes.c  */
 
 #define BLACK 1          /* background */
 
@@ -205,20 +205,18 @@
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
 #define PHASE_FACTOR 1.0       /* factor in computation of phase in color scheme P_3D_PHASE */
 #define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
-// #define COLOR_PHASE_SHIFT 0.0   /* phase shift of color scheme, in units of Pi */
 #define ATTENUATION 0.0   /* exponential attenuation coefficient of contrast with time */
 #define VSHIFT_AMPLITUDE 0.0   /* additional shift for wave amplitude */
 #define VSCALE_AMPLITUDE 0.5   /* additional scaling factor for wave amplitude */
-#define E_SCALE 50.0       /* scaling factor for energy representation */
+#define E_SCALE 10.0       /* scaling factor for energy representation */
 #define LOG_SCALE 0.75     /* scaling factor for energy log representation */
 #define LOG_SHIFT 0.75     /* shift of colors on log scale */
 #define FLUX_SCALE 250.0    /* scaling factor for energy flux represtnation */
-#define AVRG_E_FACTOR 0.75   /* controls time window size in P_AVERAGE_ENERGY scheme */
+#define AVRG_E_FACTOR 0.8   /* controls time window size in P_AVERAGE_ENERGY scheme */
 #define RESCALE_COLOR_IN_CENTER 0   /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
 #define FADE_IN_OBSTACLE 1      /* set to 1 to fade color inside obstacles */
-#define SHADE_2D 1       /* set to 1 to add pseudo-3d shading effect */ 
+#define SHADE_2D 0       /* set to 1 to add pseudo-3d shading effect */ 
-// #define SHADE_SCALE_2D 0.25  /* lower value increases sensitivity of shading */
+#define SHADE_SCALE_2D 0.25  /* lower value increases sensitivity of shading */
-#define SHADE_SCALE_2D 0.05  /* lower value increases sensitivity of shading */
 
 #define COLORHUE 260     /* initial hue of water color for scheme C_LUM */
 #define COLORDRIFT 0.0   /* how much the color hue drifts during the whole simulation */
@@ -227,9 +225,9 @@
 #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 1     /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 1.4     /* scale of color scheme bar */
+#define COLORBAR_RANGE 2.0      /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 1.2   /* scale of color scheme bar for 2nd part */
+#define COLORBAR_RANGE_B 0.2    /* scale of color scheme bar for 2nd part */
 #define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
 #define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
 #define CIRC_COLORBAR_B 0       /* set to 1 to draw circular color scheme */
@@ -599,7 +597,7 @@ void animation()
     double time, scale, ratio, startleft[2], startright[2], sign[N_SOURCES], r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle = 0.0, x1, ior_angle = 0.0, omega, phase_shift, vshift, dsource, finv, source_amp; 
     double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *average_energy[NX], *color_scale[NX], *total_flux, *tcc_table[NX], *tgamma_table[NX], *fade_table;
     short int *xy_in[NX];
-    int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q, first_source = 1, imin, imax, ij[2], source, source_period, source_shift;
+    int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q, first_source = 1, imin, imax, ij[2], source, source_period, source_shift[N_SOURCES];
 //     static int image_counter = 0;
     int image_counter = 0;
     long int wave_value;
@@ -636,7 +634,8 @@ void animation()
         init_wave_packets(packet, WAVE_PACKET_RADIUS);
     }
     
-    for (i=0; i<N_SOURCES; i++) sign[i] = 1.0;
+    if (ALTERNATE_SOURCE_PHASES) for (i=0; i<N_SOURCES; i++) sign[i] = pow(-1.0,(double)i);
+    else for (i=0; i<N_SOURCES; i++) sign[i] = 1.0;
     
     /* initialise positions and radii of circles */
     if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) ncircles = init_circle_config(circles);
@@ -792,8 +791,6 @@ void animation()
 
     glutSwapBuffers();
 
-
-
     sleep(SLEEP1);
 
     for (i=0; i<=INITIAL_TIME + NSTEPS; i++)
@@ -830,13 +827,17 @@ void animation()
         if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, fade, fade_value);
         
         /* add oscillating waves */
-        wave_source_x[0] = XMIN + 2.0*(XMAX - XMIN)*(double)i/(double)NSTEPS;
+        for (source = 0; source < N_SOURCES; source++)
-        wave_source_y[0] = 0.3;
-        if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
         {
-            if (ALTERNATE_OSCILLATING_SOURCE) sign[0] = -sign[0];
+            angle = APOLY*PID + (double)source*DPI/(double)NPOLY;
-            for (source = 0; source < N_SOURCES; source++)
+            wave_source_x[source] = LAMBDA*cos(angle) + 2.0*MU*cos(angle - 0.5*PID);
-                add_circular_wave(-sign[0]*INITIAL_AMP, wave_source_x[source], wave_source_y[source], phi, psi, xy_in);
+            wave_source_y[source] = LAMBDA*sin(angle) + 2.0*MU*sin(angle - 0.5*PID);
+            source_shift[source] = source*OSCILLATING_SOURCE_PERIOD/N_SOURCES;
+            if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == source_shift[source]))
+            {
+                if (ALTERNATE_OSCILLATING_SOURCE) sign[source] = -sign[source];
+                add_circular_wave(-sign[source]*INITIAL_AMP, wave_source_x[source], wave_source_y[source], phi, psi, xy_in);
+            }
         }
 
         

wave_sphere.c

@@ -50,14 +50,10 @@
 #define WINHEIGHT 	1150  /* window height */
 #define NX 2560          /* number of grid points on x axis */
 #define NY 1280           /* number of grid points on y axis */
-// #define NX 2048          /* number of grid points on x axis */
-// #define NY 1024          /* number of grid points on y axis */
-// #define NX 1024          /* number of grid points on x axis */
-// #define NY 512           /* number of grid points on y axis */
 
-#define DPOLE 30         /* safety distance to poles */
+#define DPOLE 15         /* safety distance to poles */
 #define SMOOTHPOLE 0.1     /* smoothing coefficient at poles */
-#define ZERO_MERIDIAN 190.0     /* choice of zero meridian (will be at left/right boundary of 2d plot) */
+#define ZERO_MERIDIAN 0.0     /* choice of zero meridian (will be at left/right boundary of 2d plot) */
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
@@ -70,20 +66,18 @@
 #define JULIA_ROT -20.0       /* rotation of Julia set, in degrees */
 #define JULIA_RE 0.5    
 #define JULIA_IM 0.462    /* parameters for Julia sets */
-// #define JULIA_IM 0.45    /* parameters for Julia sets */
 
 /* Choice of the billiard table */
 
-// #define B_DOMAIN 999         /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 81         /* choice of domain shape, see list in global_pdes.c */
-#define B_DOMAIN 84         /* choice of domain shape, see list in global_pdes.c */
 
-#define CIRCLE_PATTERN 33   /* pattern of circles or polygons, see list in global_pdes.c */
+#define CIRCLE_PATTERN 31   /* pattern of circles or polygons, see list in global_pdes.c */
 
 #define COMPARISON 0        /* set to 1 to compare two different patterns */
 #define B_DOMAIN_B 20       /* second domain shape, for comparisons */
 #define CIRCLE_PATTERN_B 0  /* second pattern of circles or polygons */
 
-#define VARIABLE_IOR 1      /* set to 1 for a variable index of refraction */
+#define VARIABLE_IOR 0      /* set to 1 for a variable index of refraction */
 #define IOR 20              /* choice of index of refraction, see list in global_pdes.c */
 #define IOR_TOTAL_TURNS 1.0 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
 #define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
@@ -94,7 +88,7 @@
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 0.75	    /* parameter controlling the dimensions of domain */
-#define MU 0.1              /* parameter controlling the dimensions of domain */
+#define MU 0.2              /* parameter controlling the dimensions of domain */
 #define MU_B 1.0            /* parameter controlling the dimensions of domain */
 #define NPOLY 6             /* number of sides of polygon */
 #define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
@@ -125,7 +119,7 @@
 
 /* Physical parameters of wave equation */
 
-#define TWOSPEEDS 1          /* set to 1 to replace hardcore boundary by medium with different speed */
+#define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
 #define OSCILLATE_LEFT 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 3  /* oscillation schedule, see list in global_pdes.c */
@@ -134,9 +128,8 @@
 #define AMPLITUDE 0.8     /* amplitude of periodic excitation */ 
 #define ACHIRP 0.2        /* acceleration coefficient in chirp */
 #define DAMPING 0.0       /* damping of periodic excitation */
-#define COURANT 0.08       /* Courant number */
+#define COURANT 0.1       /* Courant number */
-#define COURANTB 0.005     /* Courant number in medium B */
+#define COURANTB 0.005    /* Courant number in medium B */
-// #define GAMMA 2.0e-6         /* damping factor in wave equation */
 #define GAMMA 0.0         /* damping factor in wave equation */
 #define GAMMAB 0.0        /* damping factor in wave equation */
 #define GAMMA_SIDES 1.0e-4      /* damping factor on boundary */
@@ -150,9 +143,9 @@
 /* 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 ADD_OSCILLATING_SOURCE 1        /* set to 1 to add an oscillating wave source */
 #define OSCILLATING_SOURCE_PERIOD 25    /* period of oscillating source */
-#define ALTERNATE_OSCILLATING_SOURCE 0  /* set to 1 to alternate sign of oscillating source */
+#define ALTERNATE_OSCILLATING_SOURCE 1  /* set to 1 to alternate sign of oscillating source */
 
 #define ADD_WAVE_PACKET_SOURCES 0       /* set to 1 to add several sources emitting wave packets */
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
@@ -161,7 +154,7 @@
 
 #define ADD_FORCING 0       /* set to 1 to add periodic forcing */
 #define FORCING_AMP 0.0     /* amplitude of periodic forcing */
-#define FORCING_CONST_AMP 1.0-10    /* amplitude of periodic forcing */
+#define FORCING_CONST_AMP 1.0e-10    /* amplitude of periodic forcing */
 #define FORCING_PERIOD 2400  /* period of forcing */
 
 #define DRIFT_WAVE 0        /* add drift of wave to the right (experimental) */
@@ -178,10 +171,9 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 1500       /* number of frames of movie */
+// #define NSTEPS 2800       /* number of frames of movie */
-// #define NSTEPS 500         /* number of frames of movie */
+#define NSTEPS 1800         /* number of frames of movie */
-// #define NVID 6            /* number of iterations between images displayed on screen */
+#define NVID 4            /* number of iterations between images displayed on screen */
-#define NVID 10            /* number of iterations between images displayed on screen */
 #define NSEG 1000          /* number of segments of boundary */
 #define INITIAL_TIME 0      /* time after which to start saving frames */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
@@ -197,9 +189,9 @@
 
 /* Parameters of initial condition */
 
-#define INITIAL_AMP 0.03           /* amplitude of initial condition */
+#define INITIAL_AMP 2.5            /* amplitude of initial condition */
-#define INITIAL_VARIANCE 0.00005    /* variance of initial condition */
+#define INITIAL_VARIANCE 0.003    /* variance of initial condition */
-#define INITIAL_WAVELENGTH  0.002  /* wavelength of initial condition */
+#define INITIAL_WAVELENGTH  0.1   /* wavelength of initial condition */
 
 /* Plot type, see list in global_pdes.c  */
 
@@ -215,8 +207,8 @@
 #define SHADE_3D 0              /* set to 1 to change luminosity according to normal vector */
 #define SHADE_2D 1              /* set to 1 to change luminosity according to normal vector to plane */
 #define SHADE_WAVE 1            /* set to 1 to have luminosity depend on wave height */
-#define NON_DIRICHLET_BC 0      /* set to 1 to draw only facets in domain, if field is not zero on boundary */
+#define NON_DIRICHLET_BC 1      /* set to 1 to draw only facets in domain, if field is not zero on boundary */
-#define FLOOR_ZCOORD 1          /* set to 1 to draw only facets with z not too negative */
+#define FLOOR_ZCOORD 0          /* set to 1 to draw only facets with z not too negative */
 #define DRAW_BILLIARD 0         /* set to 1 to draw boundary */
 #define DRAW_BILLIARD_FRONT 0   /* set to 1 to draw front of boundary after drawing wave */
 #define DRAW_CONSTRUCTION_LINES 0   /* set to 1 to draw construction lines of certain domains */
@@ -234,14 +226,13 @@
 #define REPRESENTATION_3D 1     /* choice of 3D representation */ 
 #define PLOT_2D 1               /* switch to 2D representation, equirectangular projection */
 #define PHISHIFT 0.0            /* shift of phi in 2D plot (in degrees) */
-#define FLOODING 1              /* set to 1 to draw waves above altitude (for Earth representations) */
+#define FLOODING 0              /* set to 1 to draw waves above altitude (for Earth representations) */
 
 #define REP_AXO_3D 0        /* linear projection (axonometry) */
 #define REP_PROJ_3D 1       /* projection on plane orthogonal to observer line of sight */
 
 #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 360.0   /* total angle of rotation during simulation */
-// #define ROTATE_ANGLE -50.0   /* total angle of rotation during simulation */
 #define ROTATE_VIEW_WHILE_FADE 1    /* set to 1 to keep rotating viewpoint during fade */
 
 #define VIEWPOINT_TRAJ 1    /* type of viewpoint trajectory */
@@ -249,7 +240,7 @@
 
 /* Color schemes */
 
-#define COLOR_PALETTE 11      /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 14      /* Color palette, see list in global_pdes.c  */
 #define COLOR_PALETTE_B 16     /* Color palette, see list in global_pdes.c  */
 
 #define BLACK 1          /* background */
@@ -261,7 +252,7 @@
 
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
-#define VSCALE_AMPLITUDE 50.0  /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define VSCALE_AMPLITUDE 0.5  /* additional scaling factor for color scheme P_3D_AMPLITUDE */
 #define VSHIFT_AMPLITUDE 0.0   /* additional shift for wave amplitude */
 #define VSCALE_ENERGY 5.0     /* additional scaling factor for color scheme P_3D_ENERGY */
 #define PHASE_FACTOR 20.0      /* factor in computation of phase in color scheme P_3D_PHASE */
@@ -298,7 +289,7 @@
 #define MAZE_WIDTH 0.02     /* half width of maze walls */
 
 #define DRAW_COLOR_SCHEME 0       /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 3.0      /* scale of color scheme bar */
+#define COLORBAR_RANGE 1.5      /* scale of color scheme bar */
 #define COLORBAR_RANGE_B 6.0    /* 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 */
@@ -327,7 +318,7 @@
 #define INITIAL_SHIFT 20.0          /* time shift of initial wave packet (in oscillation periods) */
 #define WAVE_PACKET_SHIFT 200.0     /* time shift between wave packets (in oscillation periods) */
 #define FADE_IN_OBSTACLE 0      /* set to 1 to fade color inside obstacles */
-#define N_SOURCES 1                     /* number of sources, for option draw_sources */
+#define N_SOURCES 20                     /* number of sources, for option draw_sources */
 #define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
 /* end of constants only used by sub_wave and sub_maze */
 
@@ -344,28 +335,27 @@ double light[3] = {-0.816496581, -0.40824829, 0.40824829};      /* vector of "li
 double observer[3] = {-5.0, -5.0, 3.5};    /* location of observer for REP_PROJ_3D representation */ 
 int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
 
-#define ADD_DEM 1               /* add DEM (digital elevation model) */
+#define ADD_DEM 0               /* add DEM (digital elevation model) */
-#define ADD_NEGATIVE_DEM 1      /* add DEM with bathymetric data */
+#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 5            /* set to 1 to smoothen DEM (to make altitude less constant) */
 #define DEM_SMOOTH_STEPS 1      /* 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 0.0      /* sea level for flooded planet */
+#define PLANET_SEALEVEL -10000.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 TRANSPARENT_WAVE 0      /* set to 1 for waves to be "transparent" */
 
-#define RSCALE 1.0             /* scaling factor of radial component */
+#define RSCALE 1.0              /* scaling factor of radial component */
-#define RMAX 10.0               /* max value of radial component */
+#define RMAX 1.1               /* max value of radial component */
-#define Z_SCALING_FACTOR 0.8    /* overall scaling factor of z axis for REP_PROJ_3D representation */
+#define RMIN 0.9               /* min value of radial component */
-#define RMIN 0.5               /* min value of radial component */
 #define Z_SCALING_FACTOR 0.8   /* overall scaling factor of z axis for REP_PROJ_3D representation */
 #define XY_SCALING_FACTOR 2.2   /* 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 COS_VISIBLE -0.3        /* limit on cosine of normal to shown facets */
+#define COS_VISIBLE -0.02        /* limit on cosine of normal to shown facets */
 
 #include "global_pdes.c"        /* constants and global variables */
 #include "global_3d.c"          /* constants and global variables */
@@ -719,7 +709,7 @@ void viewpoint_schedule(int i)
 
 void animation()
 {
-    double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle = 0.0, lambda1, y, x1, sign1, omega, phase_shift, theta, amp; 
+    double time, scale, ratio, startleft[2], startright[2], sign = -1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, angle = 0.0, lambda1, y, x1, sign1, omega, phase_shift, theta, amp; 
     double *phi, *psi, *tmp, *color_scale, *tc, *tcc, *tgamma;
 //     double *total_energy;
     short int *xy_in;
@@ -841,11 +831,11 @@ void animation()
     
 //     init_moving_tidal_wave_sphere(0.0, 0.05, 0.025, phi, psi, xy_in, wsphere);
 
-    init_moving_tidal_wave_sphere(0.0, 2.5, 1.25, phi, psi, xy_in, wsphere);
+//     init_moving_tidal_wave_sphere(0.0, 2.5, 1.25, phi, psi, xy_in, wsphere);
     
 //     init_circular_wave_sphere(0.0, -0.5, phi, psi, xy_in, wsphere);
     
-//      init_wave_flat_sphere(phi, psi, xy_in, wsphere);
+     init_wave_flat_sphere(phi, psi, xy_in, wsphere);
     
 //     init_wave_flat_sphere(phi, psi, xy_in, wsphere);
 //     for (j=0; j<6; j++)
@@ -854,6 +844,31 @@ void animation()
 //     add_circular_wave_sphere(1.0, 1.0 - PI/9.0 + DPI/3.0, 0.15, phi, psi, xy_in, wsphere);
 //     add_circular_wave_sphere(1.0, 1.0 - PI/9.0 + 2.0*DPI/3.0, 0.15, phi, psi, xy_in, wsphere);
     
+    for (j=0; j<5; j++)
+    {
+        a = circ_sphere[1].x + circ_sphere[2].x;
+        b = circ_sphere[1].y + circ_sphere[2].y;
+        c = 1.0 + circ_sphere[1].z + circ_sphere[2].z;
+        theta = acos(c/sqrt(a*a + b*b + c*c));
+        
+        wave_source_x[4*j] = (double)j*DPI/5.0;
+        wave_source_y[4*j] = PID - theta;
+        wave_source_x[4*j+1] = (double)j*DPI/5.0 + PI/5.0;
+        wave_source_y[4*j+1] = theta - PID;
+        
+        a = circ_sphere[1].x + circ_sphere[2].x + circ_sphere[6].x;
+        b = circ_sphere[1].y + circ_sphere[2].y + circ_sphere[6].y;
+        c = circ_sphere[1].z + circ_sphere[2].z + circ_sphere[6].z;
+        theta = acos(c/sqrt(a*a + b*b + c*c));
+        
+        wave_source_x[4*j+2] = (double)j*DPI/5.0 + PI/5.0;
+        wave_source_y[4*j+2] = theta - PID;
+        wave_source_x[4*j+3] = (double)j*DPI/5.0;
+        wave_source_y[4*j+3] = PID - theta;
+        
+    }
+    
+    
 //     printf("Wave initialized\n");
 
     /* initialize table of wave speeds/dissipation */
@@ -958,12 +973,11 @@ void animation()
         if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, fade, fade_value); 
         
         /* add oscillating waves */
-//         if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == OSCILLATING_SOURCE_PERIOD - 1))
         if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
         {
             if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
-            for (j=0; j<6; j++)
+            for (j=0; j<N_SOURCES; j++)
-                add_circular_wave_sphere(sign, (double)j*PI/3.0, 0.925*PID, phi, psi, xy_in, wsphere);
+                add_circular_wave_sphere(sign,wave_source_x[j], wave_source_y[j], phi, psi, xy_in, wsphere);
         }
         if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0);