Add files via upload

2026-04-11 16:12:23 +02:00
1 parent 77e9c16fd0
commit 96f1e1760f
@@ -16,6 +16,8 @@
 #define E_KURAMOTO_SPHERE 91    /* Kuramoto on the sphere, with cubic grid */
 #define E_KELLER_SEGEL 10   /* Keller-Segel model */
 #define E_KELLER_SEGEL_SPHERE 11    /* Keller-Segel model on the sphere, cubic grid */
 #define E_GRAY_SCOTT 12     /* Gray-Scott model */
 #define E_GRAY_SCOTT_SPHERE 121  /* Gray-Scott model on the sphere, cubic grid */
 /* Choice of potential */
@@ -102,7 +104,7 @@
 #define RDE_PLANET (((ADAPT_STATE_TO_BC)&&((OBSTACLE_GEOMETRY == D_SPHERE_EARTH)||(OBSTACLE_GEOMETRY == D_SPHERE_MARS)||(OBSTACLE_GEOMETRY == D_SPHERE_VENUS))))
 // #define RDE_PLANET (((ADAPT_STATE_TO_BC)&&((OBSTACLE_GEOMETRY == D_SPHERE_EARTH)||(OBSTACLE_GEOMETRY == D_SPHERE_MARS)||(OBSTACLE_GEOMETRY == D_SPHERE_VENUS)))||((RDE_EQUATION == E_SHALLOW_WATER)&&(SWATER_DEPTH == SH_EARTH)))
-#define NMAXCIRC_SPHERE 100     /* max number of circles on sphere */
+#define NMAXCIRC_SPHERE 1000     /* max number of circles on sphere */
 #define NMAX_TRACER_PTS 20      /* max number of tracer points recorded per cell */
 #define NMAX_SPHERE_NEIGHB 10    /* max number of neighbours in Poisson simulation grid */
@@ -157,6 +159,7 @@ typedef struct
    double depth;               /* water depth */
    double cos_depth_angle;     /* cos of angle of depth profile */
    double gradx, grady;        /* gradient of water depth */
    double time;                /* time when threshold is reached */
    short int tracer;           /* has value 1 if cell contains a tracer */
    short int n_tracer_pts;     /* number of recorded tracer points per cell */
    double tracerx[NMAX_TRACER_PTS], tracery[NMAX_TRACER_PTS];    /* coordinates of tracer */
@@ -195,6 +198,7 @@ typedef struct
    int convert_grid;           /* convert field from simulation grid to longitude-latitude */
    short int edge;             /* has value 1 on edges of cubic simulation grid */
    double cos_grid, sin_grid;  /* cosine and sine of grid angle */
    double laplace_weight;      /* weight for adjusted Laplacian */
 } t_wave_sphere;
@@ -15,7 +15,7 @@
 #define NMAXCIRCLES 200000       /* total number of circles/polygons (must be at least NCX*NCY for square grid) */
 #define MAXNEIGH 20         /* max number of neighbours kept in memory */
-#define NMAXOBSTACLES 5000  /* max number of obstacles */
+#define NMAXOBSTACLES 10000  /* max number of obstacles */
 #define NMAXSEGMENTS 1200   /* max number of repelling segments */
 #define NMAXGROUPS 100      /* max number of groups of segments */
 #define NMAXCOLLISIONS 200000   /* max number of collisions */
@@ -26,7 +26,7 @@
 #define NMAXSHOVELS 50      /* max number of shovels */
 #define NMAX_TRIANGLES_PER_OBSTACLE 10  /* max number of triangles per obstacle */
 #define NMAX_TRIANGLES_PER_FACET 4  /* max number of triangles per facet between obstacles */
-#define NMAX_ABSORBERS 10   /* max number of absorbing discs */
+#define NMAX_ABSORBERS 100   /* max number of absorbing discs */
 #define C_SQUARE 0          /* square grid of circles */
 #define C_HEX 1             /* hexagonal/triangular grid of circles */
@@ -88,6 +88,7 @@
 #define S_RECTANGLE 0       /* segments forming a rectangle */
 #define S_CYLINDRICAL 100   /* two lines at top and bottom */
 #define S_BOTTOM 110        /* one line at bottom */
 #define S_VERTICAL 111      /* vertical double line */
 #define S_CUP 1             /* segments forming a cup (for increasing gravity) */
 #define S_HOURGLASS 2       /* segments forming an hour glass */
 #define S_PENTA 3           /* segments forming a pentagon with 3 angles of 120° and 2 right angles */
@@ -347,6 +348,12 @@
 #define IP_X 0      /* color depends on x coordinate of initial position */
 #define IP_Y 1      /* color depends on y coordinate of initial position */
 /* Position-dependent initial conditions (for POSITION_DEPENDENT_TYPE) */
 #define PDIC_TWO 1        /* two layers */
 #define PDIC_THREE 2      /* three layers */
 #define PDIC_MASSY 3      /* mass depends on y coordinate*/
 /* Space dependence of electric field */
 #define EF_CONST 0  /* constant magnetic field */
@@ -396,6 +403,7 @@
 #define BG_POS_OBSTACLES 8  /* background color depends on displacement of obstacles */
 #define BG_CURRENTX 9       /* background color depends on x-component of current */
 #define BG_DIRECTION 10     /* background color depends on particle orientation */
 #define BG_POTENTIAL 12     /* potential added to sphere */
 /* Obstacle color schemes */
@@ -422,6 +430,18 @@
 #define AP_CUBE 2       /* absorbers on the vertices of a cube */
 #define AP_DODECA 3     /* absorbers on the vertices of a dodecahedron */
 #define AP_ICOSA 4      /* absorbers on the vertices of an icosahedron */
 #define AP_PDISC 9      /* Poisson disc process */
 #define AP_POLES 10     /* absorbers at the North and South pole */
 /* Type of sphere potential */
 #define SPP_WELLS 0      /* potential wells located at points in table wells[] */
 #define SPP_CUBE 1       /* potential having shape of a cube */
 /* Type of sphere gravity */
 #define SG_POLAR 1       /* gravity moving from pole to pole */
 #define SG_EQUATOR 2     /* gravity moving along equator */
 /* Color schemes */
@@ -706,6 +726,7 @@ typedef struct
    double fboundary;           /* boundary force */
    double pressure;            /* pressure */
    double gravity;             /* gravity */
    double g_angle;             /* angle of gravity (for sphere) */
    double radius;              /* particle radius */
    double angle;               /* orientation of obstacle */
    double omega;               /* angular speed of obstacle */
@@ -714,6 +735,7 @@ typedef struct
    double prop;                /* proportion of types */
    double thermo_radius;       /* radius of thermostat region */
    double current;             /* electrical current */
    double omega_sphere;        /* angular frequency of rotating sphere */
 } t_lj_parameters;
@@ -725,12 +747,16 @@ typedef struct
    double reg_cottheta;        /* regularized cotangent of theta */
    double x, y, z;             /* x, y, z coordinates of point on sphere */
    double radius;              /* radius with wave height */
    double initial_radius;      /* radius before particles are added */
    double r, g, b;             /* RGB values for image */
    double cos_angle;           /* cosine of light angle */
    double cos_angle_sphere;    /* cosine of light angle for perfect sphere */
    double cos_angle_pot;       /* cosine of light angle for 2D potential */
    double potential;           /* value of potential */
    double nablax, nablay;      /* gradient of potential */
    short int locked;           /* has value 1 if color of pixel is fixed */
 } t_lj_sphere;
-int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0, n_otriangles = 0, n_ofacets = 0, nabsorbers = 0;
+int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0, n_tracers = 0, n_otriangles = 0, n_ofacets = 0, nabsorbers = 0, nwells = 0;
 FILE *lj_log;
@@ -146,6 +146,7 @@
 #define D_DISC_WAVEGUIDE_ELLIPSE 982     /* an ellipse with a shifted attached waveguide */
 #define D_WAVEGUIDE_ELLIPSE_OBLIQUE 983     /* an ellipse with an angled attached waveguide */
 #define D_CIRCLE_PAIRS 99       /* pairs of circles connected by channels */
 #define D_CIRCLE_POLYGON 991    /* circles connected to polygons by channels */
 /* for wave_sphere.c */
@@ -208,6 +209,8 @@
 #define C_SPH_ICOSA 31      /* icosahedron (on sphere) */
 #define C_SPH_OCTA 32       /* octahedron (on sphere) */
 #define C_SPH_CUBE 33       /* cube (on sphere) */
 #define C_SPH_CUBE_B 34     /* rotated cube (on sphere) */
 #define C_SPH_POISSON 35    /* Poisson process (on sphere) */
 #define C_ONE 97            /* one single circle, as for Sinai */
 #define C_TWO 98            /* two concentric circles of different type */
@@ -415,6 +418,9 @@
 #define Z_SWATER_SPEED 71            /* speed */
 #define Z_SWATER_DIRECTION_SPEED 74  /* hue for direction of velocity, luminosity for speed */
 /* for Gray-Scott model */
 #define Z_TIME_SLICE 80           /* first time field 0 reaches threshold */
 /* special boundary conditions for Euler equation */
 #define BCE_TOPBOTTOM 1         /* special flow at top and bottom */
 #define BCE_TOPBOTTOMLEFT 2     /* special flow at top, bottom and left side */
@@ -58,18 +58,18 @@
 #define YMIN 0.0
 #define YMAX 3.141592654	/* y interval for 9/16 aspect ratio */
-#define INITXMIN 2.95
+#define INITXMIN 0.0
-#define INITXMAX 3.4	/* x interval for initial condition */
+#define INITXMAX 6.28	/* x interval for initial condition */
-#define INITYMIN 1.37
+#define INITYMIN 0.0
-#define INITYMAX 1.77	/* y interval for initial condition */
+#define INITYMAX 3.14	/* y interval for initial condition */
 #define THERMOXMIN -1.25
 #define THERMOXMAX 1.25	/* x interval for initial condition */
 #define THERMOYMIN 0.0
 #define THERMOYMAX 0.75	/* y interval for initial condition */
-#define ADDXMIN 0.0
+#define ADDXMIN 0.1
-#define ADDXMAX 0.1	/* x interval for adding particles */
+#define ADDXMAX 0.2	/* x interval for adding particles */
 #define ADDYMIN 1.57
 #define ADDYMAX 1.57	/* y interval for adding particles */
 #define ADDRMIN 2.0 
@@ -85,7 +85,7 @@
 #define OBSYMIN -1.125
 #define OBSYMAX 1.125     /* x interval for motion of obstacle */
-#define CIRCLE_PATTERN 20  /* pattern of circles, see list in global_ljones.c */
+#define CIRCLE_PATTERN 81  /* pattern of circles, see list in global_ljones.c */
 #define ADD_INITIAL_PARTICLES 0 /* set to 1 to add a second type of particles */
 #define CIRCLE_PATTERN_B 0  /* pattern of circles for additional particles */
@@ -108,7 +108,7 @@
 #define COUPLE_MINLENGTH 0.5        /* length at which bonds decouple */
 #define ADD_FIXED_SEGMENTS 0    /* set to 1 to add fixed segments as obstacles */
-#define SEGMENT_PATTERN 100     /* pattern of repelling segments, see list in global_ljones.c */
+#define SEGMENT_PATTERN 111     /* pattern of repelling segments, see list in global_ljones.c */
 #define ROCKET_SHAPE 3        /* shape of rocket combustion chamber, see list in global_ljones.c */
 #define ROCKET_SHAPE_B 3      /* shape of second rocket */
 #define NOZZLE_SHAPE 6        /* shape of nozzle, see list in global_ljones.c */
@@ -126,21 +126,21 @@
 #define CENTER_PY 0         /* set to 1 to center vertical momentum */
 #define CENTER_PANGLE 0     /* set to 1 to center angular momentum */
-#define INTERACTION 1       /* particle interaction, see list in global_ljones.c */
+#define INTERACTION 1        /* particle interaction, see list in global_ljones.c */
-#define INTERACTION_B 1     /* particle interaction for second type of particle, see list in global_ljones.c */
+#define INTERACTION_B 1      /* particle interaction for second type of particle, see list in global_ljones.c */
 #define SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */
 #define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */
 #define MOL_ANGLE_FACTOR 1.0    /* rotation angle for P_MOL_ANGLE color scheme */
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 100        /* number of points for Poisson C_RAND_POISSON arrangement */
-#define PDISC_DISTANCE 5.0  /* minimal distance in Poisson disc process, controls density of particles */
+#define PDISC_DISTANCE 1.8  /* minimal distance in Poisson disc process, controls density of particles */
 #define PDISC_CANDIDATES 100 /* number of candidates in construction of Poisson disc process */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 #define LAMBDA 0.2	    /* parameter controlling the dimensions of domain */
-#define MU 0.02 	    /* parameter controlling radius of particles */
+#define MU 0.018 	    /* parameter controlling radius of particles */
-#define MU_B 0.02           /* parameter controlling radius of particles of second type */
+#define MU_B 0.018          /* parameter controlling radius of particles of second type */
 #define MU_ADD 0.022        /* parameter controlling radius of added particles */
 #define MU_ADD_B 0.022      /* parameter controlling radius of added particles */
 #define NPOLY 3             /* number of sides of polygon */
@@ -171,7 +171,7 @@
 /* Parameters for length and speed of simulation */
-#define NSTEPS 3800      /* number of frames of movie */
+#define NSTEPS 1750      /* number of frames of movie */
 #define NVID 100          /* number of iterations between images displayed on screen */
 #define NSEG 25          /* number of segments of boundary of circles */
 #define INITIAL_TIME 0     /* time after which to start saving frames */
@@ -193,15 +193,17 @@
 /* Plot type, see list in global_ljones.c  */
-#define PLOT 13
+#define PLOT 14
 #define PLOT_B 13        /* plot type for second movie */
 /* Background color depending on particle properties */
 #define COLOR_BACKGROUND 0  /* set to 1 to color background */
-#define BG_COLOR 1         /* type of background coloring, see list in global_ljones.c */
+#define BG_COLOR 2         /* type of background coloring, see list in global_ljones.c */
 #define BG_COLOR_B 3        /* type of background coloring, see list in global_ljones.c */
 #define OBSTACLE_COLOR 0    /* type of obstacle, see OC_ in global_ljones.c */
 #define SHADE_BG_COLOR_2D 1 /* set to 1 to shade BG color, for option BG_POTENTIAL */
 #define SHADE_SCALE_BG_2D 0.1   /* controls 2D shading */
 #define DRAW_BONDS 0    /* set to 1 to draw bonds between neighbours */
 #define COLOR_BONDS 1   /* set to 1 to color bonds according to length */
@@ -212,7 +214,7 @@
 #define ALTITUDE_LINES 0    /* set to 1 to add horizontal lines to show altitude */
 #define COLOR_SEG_GROUPS 0  /* set to 1 to collor segment groups differently */
 #define N_PARTICLE_COLORS 5   /* number of colors for P_NUMBER color scheme */
-#define INITIAL_POS_TYPE 0     /* type of initial position dependence */
+#define INITIAL_POS_TYPE 1     /* type of initial position dependence */
 #define ERATIO 0.995          /* ratio for time-averaging in P_EMEAN color scheme */
 #define DRATIO 0.999          /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */
 #define OBSTACLE_AREA_SHADE_FACTOR 75.0     /* controls sensitivity of triangle shade for option FILL_OBSTACLE_TRIANGLES */
@@ -222,7 +224,7 @@
 #define COLOR_PALETTE 10             /* Color palette, see list in global_ljones.c  */
 #define COLOR_PALETTE_EKIN 10        /* Color palette for kinetic energy */
-#define COLOR_PALETTE_ANGLE 0        /* Color palette for angle representation */
+#define COLOR_PALETTE_ANGLE 17       /* Color palette for angle representation */
 #define COLOR_PALETTE_DIRECTION 0    /* Color palette for direction representation */
 #define COLOR_PALETTE_INITIAL_POS 10 /* Color palette for initial position representation */
 #define COLOR_PALETTE_DIFFNEIGH 10   /* Color palette for different neighbours representation */
@@ -262,17 +264,17 @@
 #define PRINT_SEGMENTS_FORCE 0      /* set to 1 to print force on segments */
 #define PRINT_NPARTICLES 0          /* print number of active particles */
 #define PRINT_TYPE_PROP 0           /* print type proportion */
-#define PRINT_NABSORBED 1           /* print number of absorbed particles */
+#define PRINT_NABSORBED 0           /* print number of absorbed particles */
 #define FORCE_FACTOR 0.1            /* factor controlling length of force vector */
 /* particle properties */
-#define ENERGY_HUE_MIN 330.0        /* color of original particle */
+#define ENERGY_HUE_MIN 350.0        /* color of original particle */
 #define ENERGY_HUE_MAX 50.0         /* color of saturated particle */
 #define PARTICLE_HUE_MIN 359.0      /* color of original particle */
 #define PARTICLE_HUE_MAX 0.0        /* color of saturated particle */
 #define PARTICLE_EMIN 0.0           /* energy of particle with coolest color */
-#define PARTICLE_EMAX 250000.0        /* energy of particle with hottest color */
+#define PARTICLE_EMAX 4000.0        /* energy of particle with hottest color */
 #define PARTICLE_DMIN 200.0         /* energy of particle with largest local density */
 #define PARTICLE_DMAX 500.0         /* energy of particle with largest local density */
 #define SEGMENT_HUE_MIN 275.0       /* color of original segment */
@@ -280,7 +282,7 @@
 #define OBSTACLE_EMAX 1000000.0         /* energy of obstacle with hottest color */
 #define OBSTACLE_VMAX 4.0           /* speed of obstacle with largest luminosity */
 #define HUE_TYPE0 320.0      /* hue of particles of type 0 */
-#define HUE_TYPE1 60.0       /* hue of particles of type 1 */
+#define HUE_TYPE1 40.0       /* hue of particles of type 1 */
 #define HUE_TYPE2 100.0      /* hue of particles of type 2 */
 #define HUE_TYPE3 140.0      /* hue of particles of type 3 */
 #define HUE_TYPE4 180.0       /* hue of particles of type 4 */
@@ -291,6 +293,7 @@
 #define BG_LOG_EKIN_SHIFT 1.0    /* constant in BG_LOG_EKIN background color scheme */
 #define BG_FORCE_SLOPE 1.0e-6    /* constant in BG_FORCE backgound color scheme */
 #define BG_CHARGE_SLOPE 1.0     /* constant in BG_CHARGE backgound color scheme (default: 0.5) */
 #define BG_POTENTIAL_SLOPE 50.0  /* constant in BG_POTENTIAL background color scheme */
 #define CHARGE_HUE_RANGE 0.5    /* range of charge colors */
 #define PARTICLE_LMAX 1.5e4     /* angular momentum particle with brightest color */
@@ -300,14 +303,14 @@
 #define ADAPT_MASS_TO_RADIUS 0   /* set to positive value to for mass prop to power of radius */
 #define ADAPT_DAMPING_TO_RADIUS 0.0   /* set to positive value to for friction prop to power of radius */
 #define ADAPT_DAMPING_FACTOR 0.0    /* factor by which damping is adapted to radius */
-#define DT_PARTICLE 1.0e-6    /* time step for particle displacement */
+#define DT_PARTICLE 2.0e-6    /* time step for particle displacement */
 #define KREPEL 40.0           /* constant in repelling force between particles */
-#define EQUILIBRIUM_DIST 2.5    /* Lennard-Jones equilibrium distance */
+#define EQUILIBRIUM_DIST 3.5    /* Lennard-Jones equilibrium distance */
-#define EQUILIBRIUM_DIST_B 2.5  /* Lennard-Jones equilibrium distance for second type of particle */
+#define EQUILIBRIUM_DIST_B 4.5  /* Lennard-Jones equilibrium distance for second type of particle */
 #define SEGMENT_FORCE_EQR 1.0   /* equilibrium distance factor for force from segments (default 1.5) */
 #define REPEL_RADIUS 25.0    /* radius in which repelling force acts (in units of particle radius) */
-#define DAMPING 0.0          /* damping coefficient of particles */
+#define DAMPING 25.0          /* damping coefficient of particles */
-#define INITIAL_DAMPING 1000.0  /* damping coefficient of particles during initial phase */
+#define INITIAL_DAMPING 0.0  /* damping coefficient of particles during initial phase */
 #define DAMPING_ROT 0.0      /* damping coefficient for rotation of particles */
 #define DAMPING_PAIRS 0.0    /* damping between paired particles */
 #define PARTICLE_MASS 2.0    /* mass of particle of radius MU */
@@ -316,39 +319,43 @@
 #define PARTICLE_ADD_MASS_B 1.0   /* mass of added particles */
 #define PARTICLE_INERTIA_MOMENT 0.1     /* moment of inertia of particle */
 #define PARTICLE_INERTIA_MOMENT_B 0.1     /* moment of inertia of second type of particle */
-#define V_INITIAL 0.0        /* initial velocity range */
+#define V_INITIAL 50.0        /* initial velocity range */
-#define V_INITIAL_ADD 5500.0        /* initial velocity range for added particles */
+#define V_INITIAL_ADD 4500.0        /* initial velocity range for added particles */
 #define OMEGA_INITIAL 100.0        /* initial angular velocity range */
 #define VICSEK_VMIN 1.0    /* minimal speed of particles in Vicsek model */
 #define VICSEK_VMAX 40.0    /* minimal speed of particles in Vicsek model */
 #define COULOMB_LJ_FACTOR 1.0   /* relative intensity of LJ interaction in I_COULOMB_LJ interaction (default: 0.01) */
 #define KCOULOMB_FACTOR 500.0  /* relative intensity of Coulomb interaction in I_COULOMB_LJ (default: 100.0) */
 #define COULOMB_ALWAYS_REPEL 1  /* set to 1 to always repel with I_COULOMB_IMAGINARY */
 #define OBSTACLE_DAMPING 0.0   /* damping of oscillating obstacles */
 #define V_INITIAL_TYPE 0    /* type of initial speed distribution (see VI_ in global_ljones.c) */
-#define THERMOSTAT 0        /* set to 1 to switch on thermostat */
+#define THERMOSTAT 1        /* set to 1 to switch on thermostat */
 #define VARY_THERMOSTAT 0   /* set to 1 for time-dependent thermostat schedule */
 #define SIGMA 5.0           /* noise intensity in thermostat */
-#define BETA 0.00005           /* initial inverse temperature */
+#define BETA 0.002          /* initial inverse temperature */
 #define MU_XI 0.005           /* friction constant in thermostat */
 #define KSPRING_BOUNDARY 2.0e11    /* confining harmonic potential outside simulation region */
 #define KSPRING_OBSTACLE 5.0e11    /* harmonic potential of obstacles */
-#define NBH_DIST_FACTOR 6.0       /* radius in which to count neighbours */
+#define NBH_DIST_FACTOR 6.0        /* radius in which to count neighbours */
 #define BOND_DIST_FACTOR 6.0       /* radius in which to draw bonds */
 #define GRAVITY 0.0            /* gravity acting on all particles */
 #define GRAVITY_X 0.0          /* horizontal gravity acting on all particles */
 #define CIRCULAR_GRAVITY 0     /* set to 1 to have gravity directed to center */
 #define SPHERE_GRAVITY 0       /* set to 1 to have gravity at constant angle wrt sphere */
 #define INCREASE_GRAVITY 0     /* set to 1 to increase gravity during the simulation */
 #define GRAVITY_SCHEDULE 1     /* type of gravity schedule, see list in global_ljones.c */
-#define GRAVITY_FACTOR 10.0     /* factor by which to increase gravity */
+#define GRAVITY_FACTOR 2000.0     /* factor by which to increase gravity */
-#define GRAVITY_INITIAL_TIME 250    /* time at start of simulation with constant gravity */
+#define GRAVITY_INITIAL_TIME 100    /* time at start of simulation with constant gravity */
-#define GRAVITY_RESTORE_TIME 500    /* time at end of simulation with gravity restored to initial value */
+#define GRAVITY_RESTORE_TIME 750    /* time at end of simulation with gravity restored to initial value */
 #define GRAVITY_INITIAL_ANGLE 0.0   /* initial angle for SPHERE_GRAVITY */
 #define GRAVITY_DELTA_ANGLE 1440.0   /* increase of angle for SPHERE_GRAVITY */
 #define KSPRING_VICSEK 0.2   /* spring constant for I_VICSEK_SPEED interaction */
 #define VICSEK_REPULSION 10.0    /* repulsion between particles in Vicsek model */
 #define ADD_EFIELD 0      /* set to 1 to add an electric field */
-#define EFIELD 150000.0      /* value of electric field */
+#define EFIELD 20000.0    /* value of electric field */
 #define EFIELD_Y 0.0      /* value of electric field */
 #define ADD_BFIELD 0      /* set to 1 to add a magnetic field */
 #define BFIELD 20000.0       /* value of magnetic field */
@@ -359,7 +366,7 @@
 #define INCREASE_E 0      /* set to 1 to increase electric field */
 #define OSCILLATE_E 0     /* set to 1 for oscillating electric field */
 #define E_PERIOD 1000      /* period of oscillating electric field */
-#define EFIELD_FACTOR 180000.0    /* factor by which to increase electric field */
+#define EFIELD_FACTOR 1000.0    /* factor by which to increase electric field */
 #define INCREASE_B 0     /* set to 1 to increase magnetic field */
 #define BFIELD_FACTOR 1000.0    /* factor by which to increase magnetic field */
 #define CHARGE_OBSTACLES 1      /* set to 1 for obstacles to be charged */
@@ -382,6 +389,11 @@
 #define WIND_FORCE 1.35e6    /* force of wind */
 #define WIND_YMIN -0.6      /* min altitude of region with wind */
 #define ROTATE_SPHERE 1     /* set to 1 to add Coriolis and centripetal force */
 #define OMEGA_SPHERE 4.0    /* angular frequency of rotating sphere */
 #define CHANGE_OMEGA_SPHERE 1   /* set to 1 to change sphere rotation frequency */
 #define OMEGA_SPHERE_FACTOR 200.0    /* change factor of sphere rotation frequency */
 #define ROTATION 0          /* set to 1 to include rotation of particles */
 #define COUPLE_ANGLE_TO_THERMOSTAT 0    /* set to 1 to couple angular degrees of freedom to thermostat */
 #define DIMENSION_FACTOR 0.25  /* scaling factor taking into account number of degrees of freedom */  
@@ -398,15 +410,15 @@
 #define QUADRUPOLE_RATIO 0.6  /* anisotropy in quadrupole potential */ 
 #define INCREASE_BETA 0  /* set to 1 to increase BETA during simulation */
-#define BETA_SCHEDULE 0    /* type of temperature schedule, see TS_* in global_ljones */
+#define BETA_SCHEDULE 3    /* type of temperature schedule, see TS_* in global_ljones */
-#define BETA_FACTOR 50000.0    /* factor by which to change BETA during simulation */
+#define BETA_FACTOR 0.002    /* factor by which to change BETA during simulation */
 #define TS_SLOPE 8.5          /* controls speed of change of BETA for TS_TANH schedule (default 1.0) */
 #define N_TOSCILLATIONS 1.0   /* number of temperature oscillations in BETA schedule */
 #define NO_OSCILLATION 0        /* set to 1 to have exponential BETA change only */
 #define INITIAL_CONSTANT_PHASE 200 /* initial phase in which temperature is constant */
 #define MIDDLE_CONSTANT_PHASE 0   /* middle phase in which temperature is constant */
 #define FINAL_DECREASE_PHASE 1    /* final phase in which temperature decreases */ 
-#define FINAL_CONSTANT_PHASE 200     /* final phase in which temperature is constant */
+#define FINAL_CONSTANT_PHASE 400     /* final phase in which temperature is constant */
 #define DECREASE_CONTAINER_SIZE 0   /* set to 1 to decrease size of container */
 #define SMOOTH_CONTAINER_DECREASE 1 /* set to 1 to decrease size smoothly at each simulation step */
@@ -429,10 +441,10 @@
 #define N_T_AVERAGE 1        /* size of temperature averaging window */
 #define MAX_PRESSURE 3.0e10  /* pressure shown in "hottest" color */
 #define PARTIAL_THERMO_COUPLING 0   /* set to 1 to couple only some particles to thermostat */
-#define PARTIAL_THERMO_REGION 9     /* region for partial thermostat coupling (see list in global_ljones.c) */
+#define PARTIAL_THERMO_REGION 2     /* region for partial thermostat coupling (see list in global_ljones.c) */
 #define PARTIAL_THERMO_SHIFT 0.2    /* distance from obstacle at the right of which particles are coupled to thermostat */
-#define PARTIAL_THERMO_WIDTH 1.0    /* vertical size of partial thermostat coupling */
+#define PARTIAL_THERMO_WIDTH 0.3    /* vertical size of partial thermostat coupling */
-#define PARTIAL_THERMO_HEIGHT 0.0   /* vertical size of partial thermostat coupling */
+#define PARTIAL_THERMO_HEIGHT 0.2   /* vertical size of partial thermostat coupling */
 #define PARTIAL_THERMO_RIN 0.5      /* initial radius of region without coupling */
 #define PARTIAL_THERMO_RFIN 1.3     /* final radius of region without coupling */
@@ -443,20 +455,20 @@
 #define MASS_PART_BOTTOM 10000.0 /* mass of particles at bottom */
 #define NPART_BOTTOM 100     /* number of particles at the bottom */
-#define ADD_PARTICLES 1   /* set to 1 to add particles */
+#define ADD_PARTICLES 0   /* set to 1 to add particles */
 #define ADD_REGION 0      /* shape of add regions, cf ADD_* in global_ljones */
 #define ADD_TIME 20        /* time at which to add first particle */
 #define ADD_PERIOD 10000      /* time interval between adding further particles */
 #define ADD_TYPE 1         /* type of added particles */
 #define N_ADD_PARTICLES 1  /* number of particles to add */
-#define FINAL_NOADD_PERIOD 1800  /* final period where no particles are added */
+#define FINAL_NOADD_PERIOD 100  /* final period where no particles are added */
 #define SAFETY_FACTOR 10.0  /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */
 #define ADD_ALTERNATE_CHARGE 0   /* set to 1 to randomly select sign of added charge */
 #define TIME_DEPENDENT_ADD_CHARGE 0     /* set to 1 to have added charge depend on time */
 #define ALTERNATE_CHARGE_PROPORTION 0.5    /* proportion of particles of opposite charge */
-#define TRACER_PARTICLE 1   /* set to 1 to have a tracer particle */
+#define TRACER_PARTICLE 0   /* set to 1 to have a tracer particle */
-#define N_TRACER_PARTICLES 1000    /* number of tracer particles */
+#define N_TRACER_PARTICLES 5500    /* number of tracer particles */
 #define TRACER_STEPS 5           /* number of tracer steps recorded between images */
 #define TRAJECTORY_LENGTH 7000    /* length of recorded trajectory */
 #define TRAJECTORY_DRAW_LENGTH 1000 /* length of drawn trajectory */
@@ -506,10 +518,12 @@
 #define TRACK_SEGMENT_GROUPS 0      /* set to 1 for view to track group of segments */
 #define TRACK_X_PADDING 2.0         /* distance from x boundary where tracking starts */
-#define POSITION_DEPENDENT_TYPE 0   /* set to 1 to make particle type depend on initial position */
+#define POSITION_DEPENDENT_TYPE 0   /* set to PDIC_* to make particle type depend on initial position */
-#define POSITION_Y_DEPENDENCE 0     /* set to 1 for the separation between particles to be horizontal */
+#define POSITION_Y_DEPENDENCE 1     /* set to 1 for the separation between particles to be horizontal */
 #define POSITION_DEP_SIGN -1.0      /* sign in position dependence condition */
-#define POSITION_DEP_X -0.625         /* threshold value for position-dependent type */
+#define POSITION_DEP_X 1.5         /* threshold value for position-dependent type */
 #define POSITION_DEP_Y 1.5         /* threshold value for position-dependent type */
 #define POSITION_DEP_MASS_RATIO 5.0    /* position-dependent mass factor */
 #define PRINT_ENTROPY 0     /* set to 1 to compute entropy */
 #define SPECIAL_IC 0              /* set to 1 for choosing special initial condition RD_INITIAL_COND */
@@ -602,11 +616,21 @@
 #define PARTICLE_MASS_D 1.0  /* mass or partner particle */
 #define CHARGE_D -1.0         /* charge of partner particle */
-#define ADD_ABSORBERS 1     /* set to 1 to add absorbing discs */
+#define ADD_ABSORBERS 0     /* set to 1 to add absorbing discs */
-#define ABSORBER_PATTERN 4  /* pattern of absorbers, see AP_* in global_ljones */
+#define ABSORBER_PATTERN 3  /* pattern of absorbers, see AP_* in global_ljones */
 #define ABSORBER_X 0.0
 #define ABSORBER_Y 0.0      /* coordinates of first absorber */
-#define ABSORBER_R 0.16     /* radius of absorber */
+#define ABSORBER_R 0.015     /* radius of absorber */
 #define ABSORBER_PDIST 0.4  /* parameter of Poisson disc process */
 #define ADD_POTENTIAL_SPHERE 0  /* add potential for gradient force field on sphere */
 #define DRAW_POTENTIAL_SPHERE 1 /* draw sphere radius depending on potential */
 #define SPHERE_POTENTIAL 1      /* type of sphere potential */
 #define SPHERE_POT_PATTERN 3    /* pattern of local minma of SPP_WELLS sphere potential */
 #define POT_SPHERE_AMP 1.0      /* amplitude in definition of potential on sphere */
 #define POT_SPHERE_RADIUS 0.1   /* radius in definition of potential on sphere */
 #define POT_SPHERE_SMOOTH 0.5   /* smoothing of potential well */
 #define POT_SPHERE_STRENGTH 1.0e3    /* coefficient of gradient force */
 #define NXMAZE 16      /* width of maze */
 #define NYMAZE 16      /* height of maze */
@@ -620,8 +644,8 @@
 #define FLOOR_OMEGA 0      /* set to 1 to limit particle momentum to PMAX */
 #define PMAX 1000.0        /* maximal force */
-#define HASHX 60      /* size of hashgrid in x direction */
+#define HASHX 120     /* size of hashgrid in x direction */
-#define HASHY 30      /* size of hashgrid in y direction */
+#define HASHY 60      /* size of hashgrid in y direction */
 #define HASHMAX 100   /* maximal number of particles per hashgrid cell */
 #define HASHGRID_PADDING 0.1    /* padding of hashgrid outside simulation window */
@@ -634,27 +658,28 @@
 /* constants related to evolution on a sphere */
 #define SPHERE 1        /* set to 1 to compute evolution in spherical geometry */
-#define SIN_THETA_REG 0.01  /* regularization of sin(theta) for motion on sphere */
+#define SIN_THETA_REG 0.01   /* regularization of sin(theta) for motion on sphere */
-#define POLAR_PADDING 0.05   /* region around poles that belong to the same hashcell */
+#define POLAR_PADDING 0.01   /* region around poles that belong to the same hashcell */
 #define DRAW_SPHERE 1    /* set to 1 to draw 3D sphere */
 #define DRAW_ELLIPSES_ON_SPHERE 1   /* set to 1 to draw ellipses instead of circles on sphere in 2D */
 #define NX_SPHERE 3000
 #define NY_SPHERE 1500   /* number of points on sphere */
 #define Z_SCALING_FACTOR 0.75   /* overall scaling factor of z axis for REP_PROJ_3D representation */
-#define XY_SCALING_FACTOR 1.9  /* overall scaling factor for on-screen (x,y) coordinates after projection */
+#define XY_SCALING_FACTOR 2.0  /* overall scaling factor for on-screen (x,y) coordinates after projection */
 #define FLIPX -1.0             /* set to -1 to flip left/right */
 #define ZMAX_FACTOR 1.0        /* max value of z coordinate for REP_PROJ_3D representation */
 #define XSHIFT_3D -0.0          /* overall x shift for REP_PROJ_3D representation */
 #define YSHIFT_3D -0.0          /* overall y shift for REP_PROJ_3D representation */
 #define COS_VISIBLE -0.35       /* limit on cosine of normal to shown facets */
 #define RSCALE_POTENTIAL 1.0    /* radial scaling of potential */
 #define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
-#define ROTATE_ANGLE 750.0   /* total angle of rotation during simulation */
+#define ROTATE_ANGLE 360.0   /* total angle of rotation during simulation */
 #define VIEWPOINT_TRAJ 1    /* type of viewpoint trajectory */
 #define MAX_LATITUDE 45.0   /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
 double light[3] = {-0.40824829, 0.816496581, 0.40824829};      /* vector of "light" direction for P_3D_ANGLE color scheme */
-double observer[3] = {-3.0, 0.0, 1.5};    /* location of observer for REP_PROJ_3D representation */ 
+double observer[3] = {-3.0, 2.0, 2.5};    /* location of observer for REP_PROJ_3D representation */ 
 #define NO_WRAP_BC ((BOUNDARY_COND != BC_PERIODIC)&&(BOUNDARY_COND != BC_PERIODIC_CIRCLE)&&(BOUNDARY_COND != BC_PERIODIC_TRIANGLE)&&(BOUNDARY_COND != BC_KLEIN)&&(BOUNDARY_COND != BC_PERIODIC_FUNNEL)&&(BOUNDARY_COND != BC_BOY)&&(BOUNDARY_COND != BC_GENUS_TWO))
 #define PERIODIC_BC ((BOUNDARY_COND == BC_PERIODIC)||(BOUNDARY_COND == BC_PERIODIC_CIRCLE)||(BOUNDARY_COND == BC_PERIODIC_FUNNEL)||(BOUNDARY_COND == BC_PERIODIC_TRIANGLE))
@@ -992,6 +1017,18 @@ double gravity_schedule(int i, int j)
    }
 }
 double g_angle_schedule(int i, int j)
 {
    double time;
    if (i < INITIAL_TIME + GRAVITY_INITIAL_TIME) return(0.0);
    time = ((double)(i - INITIAL_TIME - GRAVITY_INITIAL_TIME) + (double)j/(double)NVID)/(double)(NSTEPS - GRAVITY_INITIAL_TIME);
    return(GRAVITY_INITIAL_ANGLE*PI/180.0 + time*GRAVITY_DELTA_ANGLE*PI/180.0);
 //     return(PI*time);
 }
 double rotation_angle(double phase)
 {
    /* case of rotating hourglass */
@@ -1851,13 +1888,11 @@ void evolve_obstacles(t_obstacle obstacle[NMAXOBSTACLES])
 void animation()
 {
    double time, scale, diss, rgb[3], dissip, gradient[2], x, y, dx, dy, dt, xleft, xright, 
-            a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, angle, theta, sum, alpha, bfield, track_x0, track_y0, efield, efieldy, dist;
+            a, b, length, fx, fy, force[2], totalenergy = 0.0, pos[2], prop, vx, xi = 0.0, torque, torque_ij, pleft = 0.0, pright = 0.0, entropy[2], speed_ratio, xmin, xmax, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t, angle, theta, sum, alpha, bfield, track_x0, track_y0, efield, efieldy, dist, dphi_inv, dtheta_inv, sincos;
    double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds, *currents;
    double *cqx, *cqy, *cpx, *cpy, *cqangle, *cpangle;
-    int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, cl, 
+    int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2,       
-        min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, 
+        total_neighbours = 0, cl, min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0, group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt, obs, reset_obs;
        tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0,
        group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt, obs, reset_obs;
    int *particle_numbers;
    static int imin, imax;
    static short int first = 1;
@@ -1876,7 +1911,7 @@ void animation()
    t_otriangle *otriangle;
    t_ofacet *ofacet;
    t_lj_sphere *wsphere;
-    t_absorber *absorber;
+    t_absorber *absorber, *wells;
    char message[100];
    ratioc = 1.0 - ratio;
@@ -1888,16 +1923,19 @@ void animation()
    params.xmaxcontainer = BCXMAX; 
    params.fboundary = 0.0;
    params.gravity = GRAVITY;
    params.g_angle = 0.0;
    params.radius = MU;
    params.nabsorbed = 0;
    particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle));    /* particles */ 
    /* sphere data */
-    if (DRAW_SPHERE) 
+    if (SPHERE) 
    {
        wsphere = (t_lj_sphere *)malloc(NX_SPHERE*NY_SPHERE*sizeof(t_lj_sphere));    
        init_lj_sphere(wsphere);
        dphi_inv = (double)NX_SPHERE/DPI;
        dtheta_inv = (double)NY_SPHERE/PI;        
    }
    if (CLUSTER_PARTICLES) 
@@ -1949,10 +1987,20 @@ void animation()
        for (i=0; i<2*NMAXCOLLISIONS; i++) collisions[i].time = 0;
    }
    if (ADD_POTENTIAL_SPHERE)
    {
        if (SPHERE_POTENTIAL == SPP_WELLS)
        {
            wells = (t_absorber *)malloc(2*NMAX_ABSORBERS*sizeof(t_absorber));
            nwells = init_absorbers(wells, SPHERE_POT_PATTERN, POT_SPHERE_RADIUS);
        }
        init_potential_sphere(wsphere, wells);
    }
    if (ADD_ABSORBERS)
    {
        absorber = (t_absorber *)malloc(2*NMAX_ABSORBERS*sizeof(t_absorber));
-        nabsorbers = init_absorbers(absorber); 
+        nabsorbers = init_absorbers(absorber, ABSORBER_PATTERN, ABSORBER_R); 
        if (DRAW_SPHERE) draw_absorbers_sphere(absorber, wsphere); 
    }
@@ -1965,7 +2013,11 @@ void animation()
    lj_log = fopen("lj_logfile.txt", "w");
    if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle, otriangle, ofacet);
-    if (ADD_FIXED_SEGMENTS) init_segment_config(segment, conveyor_belt);
+    if (ADD_FIXED_SEGMENTS) 
    {
        init_segment_config(segment, conveyor_belt);
        if (DRAW_SPHERE) draw_segments_sphere(segment, wsphere);
    }
    if ((MOVE_SEGMENT_GROUPS)&&(ADD_FIXED_SEGMENTS)) 
    {
@@ -2050,6 +2102,8 @@ void animation()
        else params.efield = EFIELD;
        if (INCREASE_B) params.bfield = bfield_schedule(i);
        else params.bfield = BFIELD;
        if (CHANGE_OMEGA_SPHERE) params.omega_sphere = omega_sphere_schedule(i);
        else params.omega_sphere = OMEGA_SPHERE;
        if ((PARTIAL_THERMO_COUPLING)&&(PARTIAL_THERMO_REGION == TH_RING_EXPAND))
            params.thermo_radius = PARTIAL_THERMO_RIN + (double)i/(double)NSTEPS*(PARTIAL_THERMO_RFIN - PARTIAL_THERMO_RIN);
        if (DECREASE_CONTAINER_SIZE) 
@@ -2137,6 +2191,7 @@ void animation()
            if (INCREASE_GRAVITY) params.gravity = gravity_schedule(i,n); 
            if (SPHERE_GRAVITY > 0) params.g_angle = g_angle_schedule(i,n);
            if ((BOUNDARY_COND == BC_RECTANGLE_WALL)&&(i < INITIAL_TIME + WALL_TIME)) wall = 1;
            else wall = 0;
@@ -2244,7 +2299,18 @@ void animation()
 //                 }
                /* add gravity */
-                if (INCREASE_GRAVITY) 
+                if (SPHERE_GRAVITY == SG_POLAR)
                {
                    particle[j].fx += params.gravity*sin(params.g_angle)*sin(particle[j].xc)/particle[j].mass_inv;
                    particle[j].fy -= params.gravity*sin(params.g_angle)*cos(particle[j].xc)*cos(particle[j].yc)/particle[j].mass_inv;
                    particle[j].fy -= params.gravity*cos(params.g_angle)*sin(particle[j].yc)/particle[j].mass_inv;
                }
                else if (SPHERE_GRAVITY == SG_EQUATOR)
                {
                    particle[j].fx += params.gravity*sin(params.g_angle - particle[j].xc)/particle[j].mass_inv;
                    particle[j].fy += params.gravity*cos(params.g_angle - particle[j].xc)*cos(particle[j].yc)/particle[j].mass_inv;
                }
                else if (INCREASE_GRAVITY) 
                {
                    if (CIRCULAR_GRAVITY)
                    {
@@ -2268,7 +2334,11 @@ void animation()
                if (ADD_EFIELD) 
                {
                    if ((INCREASE_E)||(OSCILLATE_E)) 
-                        particle[j].fx += params.efield*particle[j].charge;
+                    {
                        if (SPHERE) /* zero E field at poles */
                            particle[j].fx += params.efield*particle[j].charge*sin(particle[j].yc);
                        else particle[j].fx += params.efield*particle[j].charge;
                    }
                    else 
                    {
                        efield = EFIELD;
@@ -2278,7 +2348,9 @@ void animation()
                            efield *= (double)partial_efield(particle[j].xc, particle[j].yc); 
                            efieldy *= (double)partial_efield(particle[j].xc, particle[j].yc); 
                        }
-                        particle[j].fx += efield*particle[j].charge;
+                        if (SPHERE) /* zero E field at poles */
                            particle[j].fx += efield*particle[j].charge*sin(particle[j].yc);
                        else particle[j].fx += efield*particle[j].charge;                            
                        particle[j].fy += efieldy*particle[j].charge;
                    }
                }
@@ -2299,6 +2371,26 @@ void animation()
                    particle[j].fx += WIND_FORCE*particle[j].radius*particle[j].mass_inv;
                }
                if ((SPHERE)&&(ROTATE_SPHERE))
                {
                    /* centrifugal force */ 
                    sincos = cos(particle[j].yc)*sin(particle[j].yc);
                    particle[j].fy += params.omega_sphere*params.omega_sphere*sincos;
                    /* Coriolis force */
                    particle[j].fx += 2.0*params.omega_sphere*cos(particle[j].yc)*particle[j].vy;
                    particle[j].fy += 2.0*params.omega_sphere*sincos*particle[j].vx;
                }
                /* add potential force on sphere */
                if ((SPHERE)&&(ADD_POTENTIAL_SPHERE))
                {
                    i0 = (int)(particle[j].xc*dphi_inv);
                    j0 = (int)(particle[j].yc*dtheta_inv);
 //                     printf("(fx, fy) = (%.3lg, %.3lg)\n", POT_SPHERE_STRENGTH*wsphere[i0+NY_SPHERE+j0].nablax, POT_SPHERE_STRENGTH*wsphere[i0+NY_SPHERE+j0].nablay);
                    particle[j].fx -= POT_SPHERE_STRENGTH*wsphere[i0*NY_SPHERE+j0].nablax;
                    particle[j].fy -= POT_SPHERE_STRENGTH*wsphere[i0*NY_SPHERE+j0].nablay;
                }
                if (FLOOR_FORCE)
                {
                    if (particle[j].fx > FMAX) particle[j].fx = FMAX;
@@ -2814,13 +2906,15 @@ void animation()
    if (ADD_ABSORBERS) free(absorber); 
    if ((ADD_POTENTIAL_SPHERE)&&(SPHERE_POTENTIAL == SPP_WELLS)) free(wells); 
    if (SAVE_TIME_SERIES) 
    {
        fclose(lj_time_series);
        fclose(lj_final_position);
    }
-    if (DRAW_SPHERE) free(wsphere); 
+    if (SPHERE) free(wsphere); 
    fclose(lj_log);
 }
@@ -111,6 +111,7 @@
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
 #define COLOR_RANGE 1.0    /* max range of color (default: 1.0) */
 #define HUE_CLOSED 350.0   /* color hue of closed cells */
 #define HUE_OPEN 200.0     /* color hue of open (dry) cells */
@@ -40,7 +40,7 @@
 #include <time.h>
 #define MOVIE 0         /* set to 1 to generate movie */
-#define DOUBLE_MOVIE 0  /* set to 1 to produce movies for wave height and energy simultaneously */
+#define DOUBLE_MOVIE 1  /* set to 1 to produce movies for wave height and energy simultaneously */
 #define SAVE_MEMORY 1           /* set to 1 to save memory when writing tiff images */
 #define NO_EXTRA_BUFFER_SWAP 1    /* some OS require one less buffer swap when recording images */
@@ -48,8 +48,8 @@
 #define WINWIDTH 	1920  /* window width */
 #define WINHEIGHT 	1150  /* window height */
-#define NX 870           /* number of grid points on x axis */
+#define NX 1600           /* number of grid points on x axis */
-#define NY 520           /* number of grid points on y axis */
+#define NY 960           /* number of grid points on y axis */
 #define HRES 1           /* factor for high resolution plots */
 #define XMIN -2.0
@@ -59,7 +59,7 @@
 /* Choice of simulated equation */
-#define RDE_EQUATION 11  /* choice of reaction term, see list in global_3d.c */
+#define RDE_EQUATION 121  /* choice of reaction term, see list in global_3d.c */
 #define NFIELDS 3       /* number of fields in reaction-diffusion equation */
 #define NLAPLACIANS 0   /* number of fields for which to compute Laplacian */
@@ -141,7 +141,7 @@
 /* Physical parameters of wave equation */
-#define DT 0.000002
+#define DT 0.000007
 #define VISCOSITY 0.02
 #define POISSON_STIFFNESS 1.0   /* stiffness of Poisson equation solver for incompressible Euler */
@@ -161,6 +161,8 @@
 #define K_KURAMOTO 20.0 /* constant in Kuramoto model */
 #define NU_KURAMOTO 0.0 /* viscosity in Kuramoto model */
 #define OMEGA_KURAMOTO 0.02  /* frequency in Kuramoto model */
 /* Keller-Segel model */
 #define A_KS 0.15        /* coupling constant in Keller-Segel model */
 #define C_KS 3.5         /* coupling constant in Keller-Segel model */
 #define D_KSU 1.0        /* organism diffusion in Keller-Segel model */
@@ -168,6 +170,17 @@
 #define KS_RSCALE 0.07   /* scaling factor for reaction term of Keller-Segel model */
 #define KS_UMAX 100.0    /* max value for u component of Keller-Segel model */
 #define KS_CHIMAX_FACTOR 30.0   /* limiter on the value of chi/D_KSU */
 /* Gray-Scott model */
 #define A_GS 0.029       /* coupling constant in Gray-Scott model */
 #define B_GS 0.06        /* coupling constant in Gray-Scott model */
 #define D_GSU 1.0        /* u diffusion in Gray-Scott model */
 #define D_GSV 2.0        /* v diffusion in Gray-Scott model */
 #define GS_RSCALE 0.1    /* scaling factor for reaction term of Gray-Scott model */
 #define GS_DSCALE 0.1    /* scaling factor for diffusion term of Gray-Scott model */
 #define PDISC_DISTANCE 0.5  /* constant for Poisson disc initial condition */
 #define V_MAZE 0.4      /* potential in walls of maze */
 #define BZQ 0.0008      /* parameter in BZ equation */
 #define BZF 1.2         /* parameter in BZ equation */
@@ -217,6 +230,16 @@
 #define RPSLZB_INITIAL_TIME 0   /* initial time during which rpslzb remains constant */
 #define RPSLZB_FINAL_TIME 500   /* final time during which rpslzb remains constant */
 #define CHANGE_A_GS 0           /* set to 1 to change first parameter in Gray-Scott model */
 #define A_GS_CHANGE -0.02       /* amount by which to decrease a_gs parameter */
 #define A_GS_INITIAL_TIME 0     /* initial time during which a_gs remains constant */
 #define A_GS_FINAL_TIME 100     /* final time during which a_gs remains constant */
 #define CHANGE_B_GS 0           /* set to 1 to change second parameter in Gray-Scott model */
 #define B_GS_CHANGE -0.01       /* amount by which to decrease b_gs parameter */
 #define B_GS_INITIAL_TIME 0     /* initial time during which b_gs remains constant */
 #define B_GS_FINAL_TIME 100     /* final time during which b_gs remains constant */
 #define CHANGE_FLOW_SPEED 0     /* set to 1 to change speed of laminar flow */
 #define IN_OUT_FLOW_BC 0          /* type of in-flow/out-flow boundary conditions for Euler equation, 0 for no b.c. */
 #define IN_OUT_BC_FACTOR 0.001    /* factor of convex combination between old and new flow */
@@ -251,7 +274,7 @@
 /* Parameters for length and speed of simulation */
 #define NSTEPS 1400       /* number of frames of movie */
-#define NVID 12           /* number of iterations between images displayed on screen */
+#define NVID 30           /* number of iterations between images displayed on screen */
 #define ACCELERATION_FACTOR 1.0 /* factor by which to increase NVID in course of simulation */
 #define DT_ACCELERATION_FACTOR 1.0 /* factor by which to increase time step in course of simulation  */
 #define MAX_DT 0.024     /* maximal value of integration step */
@@ -273,6 +296,8 @@
 #define PLOT_SPHERE 1   /* draws fields on a sphere */
 #define SIMULATION_GRID 1   /* type of simulation grid, for certain equations */
 #define OPTIMIZE_GRID 1     /* set to 1 to optimize grid by evolving grid points */
 #define OPTIMIZE_STEPS 200  /* number of grid optimization steps */
 #define WEIGHT_GRID 1       /* set to 1 to compute correcting weights for Laplacian */
 #define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
 #define ROTATE_ANGLE 360.0   /* total angle of rotation during simulation */
@@ -287,12 +312,12 @@
 /* Plot type - color scheme */
 #define CPLOT 0
-#define CPLOT_B 0
+#define CPLOT_B 80
 /* Plot type - height of 3D plot */
-#define ZPLOT 0    /* z coordinate in 3D plot */
+#define ZPLOT 80    /* z coordinate in 3D plot */
-#define ZPLOT_B 0    /* z coordinate in second 3D plot */
+#define ZPLOT_B 80   /* z coordinate in second 3D plot */
 #define AMPLITUDE_HIGH_RES 1    /* set to 1 to increase resolution of P_3D_AMPLITUDE plot */
 #define NON_DIRICHLET_BC 0      /* set to 1 to draw only facets in domain, if field is not zero on boundary */
@@ -308,12 +333,15 @@
 #define FLOODING 0              /* set to 1 for drawing water when higher than continents */
 #define FLOODING_VSHIFT -10.0    /* controls when wave is considered higher than land */
 #define FLOODING_VSHIFT_2D 0.56    /* controls when wave is considered higher than land */
 #define TIME_SLICE_THRESHOLD 0.1    /* threshold of field 0 for Z_TIME_SLICE */
 #define PLOT_SCALE_ENERGY 0.05      /* vertical scaling in energy plot */
 #define PRINT_TIME 0        /* set to 1 to print running time */
 #define PRINT_VISCOSITY 0   /* set to 1 to print viscosity */
 #define PRINT_RPSLZB 0      /* set to 1 to print rpslzb parameter */
 #define PRINT_A_GS 0        /* set to 1 to print a_gs parameter */
 #define PRINT_B_GS 0        /* set to 1 to print b_gs parameter */
 #define PRINT_PROBABILITIES 0   /* set to 1 to print probabilities (for Ehrenfest urn configuration) */
 #define PRINT_NOISE 0       /* set to 1 to print noise intensity */
 #define PRINT_FLOW_SPEED 0      /* set to 1 to print speed of flow */
@@ -337,8 +365,8 @@
 /* Color schemes, see list in global_pdes.c  */
-#define COLOR_PALETTE 15        /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 17        /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 11      /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 14      /* Color palette, see list in global_pdes.c  */
 #define BLACK 1          /* black background */
 #define COLOR_OUT_R 1.0    /* color outside domain */
@@ -352,8 +380,9 @@
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
 #define COLOR_RANGE 1.0    /* max range of color (default: 1.0) */
-#define VSHIFT_AMPLITUDE -0.75   /* additional shift for wave amplitude */
+#define VSHIFT_AMPLITUDE -0.4    /* additional shift for wave amplitude */
-#define VSCALE_AMPLITUDE 1.5      /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define VSCALE_AMPLITUDE 2.5      /* additional scaling factor for color scheme P_3D_AMPLITUDE */
 #define VSCALE_TIMESLICE 1.6      /* additional scaling factor for color scheme Z_TIMESLICE */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
 #define CURL_SCALE 1.25   /* scaling factor for curl representation */
 #define RESCALE_COLOR_IN_CENTER 0   /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
@@ -395,8 +424,8 @@
 #define MAZE_XSHIFT 0.0     /* horizontal shift of maze */
 #define MAZE_WIDTH 0.04     /* half width of maze walls */
-#define DRAW_COLOR_SCHEME 1     /* set to 1 to plot the color scheme */
+#define DRAW_COLOR_SCHEME 0     /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 2.5      /* scale of color scheme bar */
+#define COLORBAR_RANGE 1.5      /* scale of color scheme bar */
 #define COLORBAR_RANGE_B 2.5    /* scale of color scheme bar for 2nd part */
 #define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
 #define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
@@ -478,16 +507,16 @@ int reset_view = 0;         /* switch to reset 3D view parameters (for option RO
 #define VENUS_NODATA_FACTOR 0.5     /* altitude to assign to DEM points without data (fraction of mean altitude) */
 #define Z_SCALING_FACTOR 0.75   /* overall scaling factor of z axis for REP_PROJ_3D representation */
-#define XY_SCALING_FACTOR 1.9  /* overall scaling factor for on-screen (x,y) coordinates after projection */
+#define XY_SCALING_FACTOR 1.7  /* overall scaling factor for on-screen (x,y) coordinates after projection */
 #define ZMAX_FACTOR 1.0        /* max value of z coordinate for REP_PROJ_3D representation */
-#define XSHIFT_3D -0.0          /* overall x shift for REP_PROJ_3D representation */
+#define 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 YSHIFT_3D 0.0          /* overall y shift for REP_PROJ_3D representation */
 #define BORDER_PADDING 0       /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */
-#define DRAW_ARROW 1           /* set to 1 to draw arrow above sphere */
+#define DRAW_ARROW 0           /* set to 1 to draw arrow above sphere */
 #define ZMAX 3.0             /* maximal value of z coordinate */
 #define ZMIN -3.0            /* maximal value of z coordinate */
-#define RSCALE 0.02             /* scaling factor of radial component */
+#define RSCALE 0.25              /* scaling factor of radial component */
 #define RSCALE_B 1.00           /* experimental, additional radial scaling factor in ij_to_sphere */
 #define RSHIFT 0.0              /* shift in radial component */
 #define RMAX 4.0                /* max value of radial component */
@@ -496,13 +525,12 @@ int reset_view = 0;         /* switch to reset 3D view parameters (for option RO
 /* For debugging purposes only */
 #define FLOOR 1         /* set to 1 to limit wave amplitude to VMAX */
-#define VMAX 100.0        /* max value of wave amplitude */
+#define VMAX 10.0        /* max value of wave amplitude */
 #define TEST_GRADIENT 0 /* print norm squared of gradient */
 #define REFRESH_B (ZPLOT_B != ZPLOT)||(CPLOT_B != CPLOT)    /* to save computing time, to be improved */
 #define COMPUTE_WRAP_ANGLE ((WRAP_ANGLE)&&((cplot == Z_ANGLE_GRADIENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_ARGUMENT)||(cplot == Z_ANGLE_GRADIENTX)||(cplot == Z_EULER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_ANGLE)))
-#define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)||(PRINT_AVERAGE_SPEED))
+#define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(PRINT_A_GS)||(PRINT_B_GS)||(PRINT_PROBABILITIES)||(PRINT_NOISE)||(PRINT_FLOW_SPEED)||(PRINT_AVERAGE_SPEED))
 #define COMPUTE_PRESSURE ((ZPLOT == Z_EULER_PRESSURE)||(CPLOT == Z_EULER_PRESSURE)||(ZPLOT_B == Z_EULER_PRESSURE)||(CPLOT_B == Z_EULER_PRESSURE))
 #define ASYM_SPEED_COLOR (VMEAN_SPEED == 0.0)
@@ -1301,6 +1329,8 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
        first = 0;
    }
 //     printf("[evolve_wave_half], equation = %i\n", RDE_EQUATION); 
    /* smooth vector fields at poles */
 //     if ((SPHERE)&&(!SMOOTHBLOCKS)) for (k=0; k<NFIELDS; k++) smooth_poles(phi_in[k]);
 //     else for (k=0; k<NFIELDS; k++) block_poles(phi_in[k]);
@@ -1476,6 +1506,15 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
            compute_grad_product_grid(phi_in[2], chi_u, prod_gradients, wsphere, ngridpoints);            
            break;
        }
        case (E_GRAY_SCOTT_SPHERE):
        {
            delta_u = (double *)malloc(NX*NY*sizeof(double));
            delta_v = (double *)malloc(NX*NY*sizeof(double));
            compute_laplacian_grid(phi_in[1], delta_u, wsphere, ngridpoints);
            compute_laplacian_grid(phi_in[2], delta_v, wsphere, ngridpoints);
            break;
        }
        default: 
        {
            /* do nothing */
@@ -1493,8 +1532,9 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
    }
-    if (RDE_EQUATION != E_KURAMOTO_SPHERE) 
+    if ((RDE_EQUATION != E_KURAMOTO_SPHERE)&&(RDE_EQUATION != E_KELLER_SEGEL_SPHERE)&&(RDE_EQUATION != E_GRAY_SCOTT_SPHERE))
    {
 //         printf("[evolve_wave_half], iteration\n");
        #pragma omp parallel for private(i,j,k,x,y,z,deltax,deltay,deltaz,rho)
        for (i=imin; i<imax; i++){
        for (j=0; j<NY; j++){
@@ -1545,6 +1585,19 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
                    phi_out[1][i*NY+j] = y + intstep*(deltay + KS_RSCALE*(x - A_KS*y));
                    break;
                }
                case (E_GRAY_SCOTT):
                {
 //                     printf("[evolve_wave_half], GS iteration\n");
                    x = phi_in[0][i*NY+j];
                    y = phi_in[1][i*NY+j];
                    deltax = GS_DSCALE*D_GSU*delta_phi[0][i*NY+j];
                    deltay = GS_DSCALE*D_GSV*delta_phi[1][i*NY+j];
                    phi_out[0][i*NY+j] = x + intstep*(deltax + GS_RSCALE*(x*x*y - (a_gs + b_gs)*x));
                    phi_out[1][i*NY+j] = y + intstep*(deltay + GS_RSCALE*(-x*x*y + a_gs*(1.0-y)));
 //                     printf("(u,v) = (%.3lg, %.3lg)\n", x, y); 
                    break;
                }
 //                 case (E_KURAMOTO_SPHERE):
 //                 {
 //                     x = K_KURAMOTO*kuramoto_int[i*NY+j] + OMEGA_KURAMOTO;
@@ -1806,7 +1859,7 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
            else phi_out[1][i] = x;
        }
        /* convert field 1 on simulation grid to field 0 on longitude_latitude grid */
-        convert_fields_from_grids(phi_out[1], phi_out[0], wsphere);
+        convert_fields_from_grids(phi_out[1], phi_out[0], rde, wsphere);
        /* TEST */
 //         for (i=0; i<NX*NY; i++) rde[i].theta = phi_out[0][i];
@@ -1829,7 +1882,22 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
        }
        if (floor_u) printf("u floored to KS_UMAX\n");
        /* convert field 1 on simulation grid to field 0 on longitude_latitude grid */
-        convert_fields_from_grids(phi_out[1], phi_out[0], wsphere);
+        convert_fields_from_grids(phi_out[1], phi_out[0], rde, wsphere);
    }
    else if (RDE_EQUATION == E_GRAY_SCOTT_SPHERE)
    {
        #pragma omp parallel for private(i,x,y,deltax,deltay)
        for (i=0; i<ngridpoints; i++)
        {
            x = phi_in[1][i];
            y = phi_in[2][i];
            deltax = GS_DSCALE*D_GSU*delta_u[i];
            deltay = GS_DSCALE*D_GSV*delta_v[i];
            phi_out[1][i] = x + intstep*(deltax + GS_RSCALE*(x*x*y - (a_gs + b_gs)*x));
            phi_out[2][i] = y + intstep*(deltay + GS_RSCALE*(-x*x*y + a_gs*(1.0-y)));
        }
        convert_fields_from_grids(phi_out[1], phi_out[0], rde, wsphere);
    }
    /* in-flow/out-flow b.c. for incompressible Euler equation */
@@ -1907,6 +1975,12 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
            free(prod_gradients);
            break;
        }
        case (E_GRAY_SCOTT_SPHERE):
        {
            free(delta_u);
            free(delta_v);
            break;
        }
    }
    if (COMPUTE_PRESSURE) 
@@ -2097,11 +2171,13 @@ void print_parameters(double *phi[NFIELDS], t_rde rde[NX*NY], short int xy_in[NX
        if (WINWIDTH > 1280)
        {
            boxheight = 0.035;
-            boxwidth = 0.21;
+            boxwidth = 0.15;
 //             boxwidth = 0.21;
            if (left)
            {
                xbox = XMIN + 0.4;
-                xtext = XMIN + 0.2;
+//                 xtext = XMIN + 0.2;
                xtext = XMIN + 0.3;
            }
            else
            {
@@ -2159,6 +2235,7 @@ void print_parameters(double *phi[NFIELDS], t_rde rde[NX*NY], short int xy_in[NX
        if (PRINT_TIME) sprintf(message, "Time %.3f", time);
        else if (PRINT_VISCOSITY) sprintf(message, "Viscosity %.3f", viscosity);
        else if (PRINT_RPSLZB) sprintf(message, "b = %.3f", rpslzb);
        else if (PRINT_A_GS) sprintf(message, "a = %.3f", a_gs);
        else if (PRINT_NOISE) sprintf(message, "noise %.3f", noise);
        else if (PRINT_FLOW_SPEED) sprintf(message, "Speed %.3f", flow_speed);
        else if (PRINT_AVERAGE_SPEED) 
@@ -2181,6 +2258,22 @@ void print_parameters(double *phi[NFIELDS], t_rde rde[NX*NY], short int xy_in[NX
                write_text(pos[0], pos[1], message2);
            }
        }
        if (PRINT_B_GS) 
        {
            y -= 0.1;
            erase_area_hsl(xbox, y + 0.02, boxwidth, boxheight, 0.0, 0.9, 0.0);
            glColor3f(1.0, 1.0, 1.0);
            sprintf(message, "b = %.3f", b_gs);
            if (PLOT_3D) 
            {
                write_text(xtext, y, message);
            }
            else
            {
                xy_to_pos(xtext, y, pos);
                write_text(pos[0], pos[1], message);
            }
        }
    }
 }
@@ -2247,6 +2340,32 @@ double rpslzb_schedule(int i)
    }
 }
 double a_gs_schedule(int i)
 {
    double ratio;
    if (i < A_GS_INITIAL_TIME) return (A_GS);
    else if (i > NSTEPS - A_GS_FINAL_TIME) return(A_GS - A_GS_CHANGE);
    else 
    {
        ratio = (double)(i - A_GS_INITIAL_TIME)/(double)(NSTEPS - A_GS_INITIAL_TIME - A_GS_FINAL_TIME);
        return (A_GS - ratio*A_GS_CHANGE);
    }
 }
 double b_gs_schedule(int i)
 {
    double ratio;
    if (i < B_GS_INITIAL_TIME) return (B_GS);
    else if (i > NSTEPS - B_GS_FINAL_TIME) return(B_GS - B_GS_CHANGE);
    else 
    {
        ratio = (double)(i - B_GS_INITIAL_TIME)/(double)(NSTEPS - B_GS_INITIAL_TIME - B_GS_FINAL_TIME);
        return (B_GS - ratio*B_GS_CHANGE);
    }
 }
 double flow_speed_schedule(int i)
 {
    double ratio;
@@ -2395,7 +2514,7 @@ void animation()
 //     if (ADD_TRACERS) tracers = (double *)malloc(2*NSTEPS*N_TRACERS*sizeof(double));
    if (ADD_TRACERS) tracers = (double *)malloc(4*NSTEPS*N_TRACERS*sizeof(double));
-    if ((RDE_EQUATION == E_KURAMOTO_SPHERE)||(RDE_EQUATION == E_KELLER_SEGEL_SPHERE))
+    if ((RDE_EQUATION == E_KURAMOTO_SPHERE)||(RDE_EQUATION == E_KELLER_SEGEL_SPHERE)||(RDE_EQUATION == E_GRAY_SCOTT_SPHERE))
        ngridpoints = initialize_simulation_grid_sphere(wsphere);
    dx = (XMAX-XMIN)/((double)NX);
@@ -2474,8 +2593,27 @@ void animation()
 //     init_linear_blob_sphere(0, 1.3*PI, 0.65*PI, 0.0, 0.0, 0.3, 0.1, 0.1, SWATER_MIN_HEIGHT, phi, xy_in, wsphere);
-    init_random(0.5, 0.4, phi, xy_in, wsphere);
+//     init_random(0.5, 0.1, phi, xy_in, wsphere);
-//     init_onefield_random(1, 0.2, 0.1, phi, xy_in, wsphere);
+
 //     init_random_zerosum(0.4, 0.4, phi, xy_in, wsphere);
 //     init_gs_blob(0.0, 0.0, 0.0, 1.0, 0.25, phi, xy_in, wsphere);
 //     init_gs_blob(0.0, 0.0, 0.0, 1.0, 0.05, phi, xy_in, wsphere);
 //     init_gs_blob(0.0, 0.0, 0.0, 1.0, 0.25, phi, xy_in, wsphere);
 //     init_gs_blobs_hex(6, 4, 1.0, 0.01, 0.0025, phi, xy_in, wsphere);
 //     init_gs_blobs_square(1, 1, 1.0, 0.01, 0.002, phi, xy_in, wsphere);
 //     init_gs_blobs_poisson(0.5, 100, 1.0, 0.007, 0.0015, phi, xy_in, wsphere);
 //     init_gs_blobs_sphere(1.0, 0.1, 0.1, phi, xy_in, wsphere, C_SPH_CUBE_B);
    init_gs_blobs_sphere(1.0, 0.05, 0.05, phi, xy_in, wsphere, C_SPH_POISSON);
    //  init_gs_rect_sphere(1.0, PID, PI, PID-0.02, PID+0.02, phi, xy_in, wsphere, 0);
 //     init_onefield_random(0, 0.5, 0.5, phi, xy_in, wsphere);
 //     init_onefield_random(1, 0.5, 0.5, phi, xy_in, wsphere);
 //     init_random_smoothed(0.0, 0.4, phi, xy_in, wsphere);
 //     init_expanding_blob_sphere(0, (279.0/180.0)*PI, (115.0/180.0)*PI, 0.75, 0.04, 0.06, SWATER_MIN_HEIGHT, phi, xy_in, wsphere);
@@ -2503,7 +2641,7 @@ void animation()
    if (ADD_MOON_FORCING) compute_forcing_schedule(0, wsphere);
    if (RDE_EQUATION == E_KURAMOTO_SPHERE)
-        convert_fields_from_grids(phi[1], phi[0], wsphere);
+        convert_fields_from_grids(phi[1], phi[0], rde, wsphere);
    blank();
    glColor3f(0.0, 0.0, 0.0);
@@ -2513,7 +2651,7 @@ void animation()
    printf("Drawing wave\n");
-    draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+    draw_wave_rde(0, 0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
 //     draw_billiard();
    if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, VISCOSITY, noise);
    if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);
@@ -2542,6 +2680,8 @@ void animation()
            }
        }
        if (CHANGE_RPSLZB) rpslzb = rpslzb_schedule(i);
        if (CHANGE_A_GS) a_gs = a_gs_schedule(i);
        if (CHANGE_B_GS) b_gs = b_gs_schedule(i);
        if (CHANGE_FLOW_SPEED) flow_speed = flow_speed_schedule(i); 
        else flow_speed = IN_OUT_FLOW_AMP;
        if (ADD_MOON_FORCING) compute_forcing_schedule(i, wsphere);
@@ -2553,7 +2693,7 @@ void animation()
        }
        printf("Drawing wave %i\n", i);
-        draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+        draw_wave_rde(i, 0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
 //         nvid = (int)((double)NVID*(1.0 + (ACCELERATION_FACTOR - 1.0)*(double)i/(double)NSTEPS));
        /* increase integration step */
@@ -2652,7 +2792,7 @@ void animation()
            if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
            {
-                draw_wave_rde(1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
+                draw_wave_rde(i, 1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
                if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, i, 0, 1.0);
 //                 draw_billiard();
                if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
@@ -2684,7 +2824,7 @@ void animation()
    {
        if (DOUBLE_MOVIE) 
        {
-            draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
+            draw_wave_rde(NSTEPS, 0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 0, 1.0, 1);
            if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, NSTEPS, 0, 1.0);
 //             draw_billiard();
            if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
@@ -2696,7 +2836,7 @@ void animation()
            else for (i=0; i<MID_FRAMES; i++) 
            {
                fade_value = 1.0 - (double)i/(double)MID_FRAMES;
-                draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
+                draw_wave_rde(NSTEPS, 0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
                if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
 //                 draw_billiard();
                if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
@@ -2704,7 +2844,7 @@ void animation()
                if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
                save_frame_counter(NSTEPS + i + 1);
            }
-            draw_wave_rde(1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
+            draw_wave_rde(NSTEPS, 1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 0, 1.0, REFRESH_B);
            if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, NSTEPS, 0, 1.0);
            if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
            if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0); 
@@ -2714,7 +2854,7 @@ void animation()
            else for (i=0; i<END_FRAMES; i++) 
            {
                fade_value = 1.0 - (double)i/(double)END_FRAMES;
-                draw_wave_rde(1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
+                draw_wave_rde(NSTEPS, 1, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT_B, CPLOT_B, COLOR_PALETTE_B, 1, fade_value, 0);
                if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
                if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
                if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 1, fade_value);   
@@ -2728,8 +2868,9 @@ void animation()
            else for (i=0; i<END_FRAMES; i++) 
            {
                fade_value = 1.0 - (double)i/(double)END_FRAMES;
-                draw_wave_rde(0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
+                draw_wave_rde(NSTEPS, 0, phi, xy_in, rde, wsphere, wsphere_hr, potential_field, ZPLOT, CPLOT, COLOR_PALETTE, 1, fade_value, 0);
                if ((ADD_TRACERS)&&(!PLOT_3D)) draw_tracers(phi, tracers, NSTEPS, 1, fade_value);
                if (PRINT_PARAMETERS) print_parameters(phi, rde, xy_in, time, PRINT_LEFT, viscosity_printed, noise);
                if (DRAW_COLOR_SCHEME) draw_color_bar_palette(CPLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value); 
                glutSwapBuffers();
                save_frame_counter(NSTEPS + 1 + counter + i);
@@ -137,7 +137,7 @@ void init_hashcell_coordinates(t_hashgrid hashgrid[HASHX*HASHY])
                /* TODO: correct area */
                hashgrid[n].area = dx_sphere[j]*dy*sin(0.5*(hashgrid[n].y1 + hashgrid[n].y2));
-                printf("Cell %i corners (%.3lg, %.3lg), (%.3lg, %.3lg)\n", n, hashgrid[n].x1, hashgrid[n].y1, hashgrid[n].x2, hashgrid[n].y2); 
+//                 printf("Cell %i corners (%.3lg, %.3lg), (%.3lg, %.3lg)\n", n, hashgrid[n].x1, hashgrid[n].y1, hashgrid[n].x2, hashgrid[n].y2); 
            }
            for (i=hashx_sphere[j]; i<HASHX; i++)
            {
@@ -192,7 +192,7 @@ void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
            if (hashx_sphere[j] > HASHX-1) hashx_sphere[j] = HASHX-1;
            dx_sphere[j] = DPI/(double)hashx_sphere[j];
-            printf("hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]);
+//             printf("hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]);
            fprintf(lj_log, "hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]);
        }
        for (j=HASHY/2 + 1; j<HASHY; j++)
@@ -200,7 +200,7 @@ void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
            hashx_sphere[j] = hashx_sphere[HASHY-1-j];
            dx_sphere[j] = dx_sphere[HASHY-1-j];
-            printf("hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]);
+//             printf("hashx_sphere[%i] = %i, dx_sphere[%i] = %.3lg\n", j, hashx_sphere[j], j, dx_sphere[j]);
        }
    }
@@ -711,19 +711,19 @@ void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
                }
            }
-            for (i=0; i<HASHX*HASHY; i++)
+//             for (i=0; i<HASHX*HASHY; i++)
-            {
+//             {
-                p = hashgrid[i].nneighb;
+//                 p = hashgrid[i].nneighb;
-                printf("Hash cell %i has %i neighbours: ", i, p);
+//                 printf("Hash cell %i has %i neighbours: ", i, p);
-                fprintf(lj_log, "Hash cell %i has %i neighbours: ", i, p);
+//                 fprintf(lj_log, "Hash cell %i has %i neighbours: ", i, p);
-                for (j=0; j<p; j++)
+//                 for (j=0; j<p; j++)
-                {
+//                 {
-                    printf("%i ", hashgrid[i].neighbour[j]);
+//                     printf("%i ", hashgrid[i].neighbour[j]);
-                    fprintf(lj_log, "%i ", hashgrid[i].neighbour[j]);
+//                     fprintf(lj_log, "%i ", hashgrid[i].neighbour[j]);
-                }
+//                 }
-                printf("\n");
+//                 printf("\n");
-                fprintf(lj_log, "\n");
+//                 fprintf(lj_log, "\n");
-            }
+//             }
@@ -342,6 +342,15 @@ void draw_circle_precomp(double x, double y, double r)
    glEnd();
 }
 void draw_ellipse_precomp(double x, double y, double rx, double ry)
 {
    int i;
    glBegin(GL_LINE_LOOP);
    for (i=0; i<=NSEG; i++) glVertex2d(x + rx*cosangle[i], y + ry*sinangle[i]);
    glEnd();
 }
 void draw_colored_circle_precomp(double x, double y, double r, double rgb[3])
 {
    int i;
@@ -565,7 +574,7 @@ int generate_poisson_discs(t_point *point, int nmax, double xmin, double xmax, d
 /* generate a Poisson disc sample in a given rectangle */
 {
    int i, j, k, n_p_active, ncandidates=PDISC_CANDIDATES, naccepted, npoints; 
-    double r, phi, x, y;
+    double r, phi, x, y, z, x1, y1, z1, x2, y2;
    short int active_poisson[nmax], far;
    printf("Generating Poisson disc sample\n");
@@ -587,8 +596,32 @@ int generate_poisson_discs(t_point *point, int nmax, double xmin, double xmax, d
        {
            r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0);
            phi = DPI*(double)rand()/RAND_MAX;
-            x = point[i].xc + r*cos(phi);
+            
-            y = point[i].yc + r*sin(phi);
+            if (sphere)
            {
                /* random point centered at north pole */
                x = sin(r)*cos(phi);
                y = sin(r)*sin(phi);
                z = cos(r);
                /* rotation by theta = point[i].yc */
                x1 = x*cos(point[i].yc) + z*sin(point[i].yc);
                y1 = y;
                z1 = -x*sin(point[i].yc) + z*cos(point[i].yc);
                /* rotation by phi = point[i].xc */
                x2 = x1*cos(point[i].xc) - y1*sin(point[i].xc);
                y2 = x1*sin(point[i].xc) + y1*cos(point[i].xc);
                y = acos(z1);
                x = argument(x2, y2);
                if (x < 0.0) x += DPI;
            }
            else
            {
                x = point[i].xc + r*cos(phi);
                y = point[i].yc + r*sin(phi);
            }
            far = 1;
            for (k=0; k<npoints; k++) if ((k!=i))
            {
@@ -3233,6 +3266,23 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS], t_belt conveyor_belt[N
            }
            break;
        }
        case (S_VERTICAL):
        {
            width = MU;
            add_rectangle_to_segments(-width, YMAX, width, YMIN, segment, 0); 
            add_rectangle_to_segments(XMAX-width, YMAX, XMAX+width, YMIN, segment, 0); 
            cycle = 0;
            for (i=0; i<nsegments; i++) 
            {
                segment[i].concave = 0;
                segment[i].group = 0;
                segment[i].inactivate = 0;
            }
            break;
        }
        case (S_CUP):
        {
            angle = APOLY*PID;
@@ -5786,7 +5836,7 @@ int compute_particle_interaction(int i, int k, double force[2], double *torque,
        {
            charge = particle[i].charge*particle[k].charge;
            f = -KCOULOMB_FACTOR*krepel*charge/(1.0e-12 + distance*distance);
-            if (charge <= 0.0) 
+            if ((charge <= 0.0)||(COULOMB_ALWAYS_REPEL))
                f1 = COULOMB_LJ_FACTOR*krepel*lennard_jones_force(distance, particle[k]);
            else f1 = 0.0;
            force[0] = f1*ca - f*sa;
@@ -8130,7 +8180,12 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
    }
    else if (cont) 
    {
-        if (nsides == NSEG) draw_colored_circle_precomp(xc1, yc1, radius, rgb);
+        if (nsides == NSEG) 
        {
            if ((SPHERE)&&(DRAW_ELLIPSES_ON_SPHERE))
                draw_colored_ellipse_precomp(xc1, yc1, radius/(sin(yc1) + SIN_THETA_REG), radius, rgb);
            else draw_colored_circle_precomp(xc1, yc1, radius, rgb);
        }
        else draw_colored_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID, rgb);
    }
@@ -8193,7 +8248,12 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
    }
    else if (cont) 
    {
-        if (nsides == NSEG) draw_circle_precomp(xc1, yc1, radius);
+        if (nsides == NSEG) 
        {
            if ((SPHERE)&&(DRAW_ELLIPSES_ON_SPHERE))
                draw_ellipse_precomp(xc1, yc1, radius/(sin(yc1) + SIN_THETA_REG), radius);
            else draw_circle_precomp(xc1, yc1, radius);
        }
        else draw_polygon(xc1, yc1, radius, nsides, angle + APOLY*PID); 
    }
@@ -8786,44 +8846,12 @@ void draw_background(t_particle particle[NMAXCIRCLES], t_segment segment[NMAXSEG
    double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue, valx, valy, lum;
    static int first = 1, counter = 0;
    static double area_factor;
-    
+      
 //     if (first)
 //     {
 //         area_factor = hashgrid[mhash(0,HASHY/2)].area;
 //         first = 0;
 //     }
 //     p1 = 0.75;
 //     p2 = 1.0 - p1;
 // //     pp1 = 0.95;
 //     pp1 = 0.99;
 //     pp2 = 1.0 - pp1;
    glBegin(GL_QUADS);
    for (i=0; i<HASHX; i++)
        for (j=0; j<HASHY; j++)
        {
            n = mhash(i, j);
 //             if (first) 
 //             {
 //                 hashgrid[n].hue1 = 180.0;
 //                 hashgrid[n].hue2 = 180.0;
 //             }
 //             /* set two old values for option DOUBLE_MOVIE */
 //             if (DOUBLE_MOVIE)
 //             {
 //                 if (counter) oldhue = hashgrid[n].hue1;
 //                 else oldhue = hashgrid[n].hue2;
 //             }
 //             else oldhue = hashgrid[n].hue1;
 //             if (DOUBLE_MOVIE)
 //             {
 //                 if (counter) hashgrid[n].hue1 = hue;
 //                 else hashgrid[n].hue2 = hue;
 //             }
 //             else hashgrid[n].hue1 = hue;
            glColor3f(hashgrid[n].r, hashgrid[n].g, hashgrid[n].b);
            glVertex2d(hashgrid[n].x1, hashgrid[n].y1);
@@ -8832,10 +8860,51 @@ void draw_background(t_particle particle[NMAXCIRCLES], t_segment segment[NMAXSEG
            glVertex2d(hashgrid[n].x1, hashgrid[n].y2);
        }
    glEnd();
 //     first = 0;
 //     counter = 1 - counter;
 }
 void draw_potential(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
 /* color background according to potential on sphere */
 {
    int i, j, k, cell, res = 3;
    double pot, value, phi, theta, hue, rgb[3];
    static double dphi, dtheta, dphi_res, dtheta_res;
    static int first = 1;
    if (first)
    {
        dphi = (XMAX-XMIN)/(double)NX_SPHERE;
        dtheta = (YMAX-YMIN)/(double)NY_SPHERE;
        dphi_res = (double)res*dphi;
        dtheta_res = (double)res*dtheta;
        first = 0;
    }
    glBegin(GL_QUADS);
    for (i=0; i<NX_SPHERE; i+=res)
        for (j=0; j<NY_SPHERE; j+=res)
        {
            cell = i*NY_SPHERE + j;
            pot = wsphere[cell].potential;
            value = atan(BG_POTENTIAL_SLOPE*pot)/PI + 0.5;
 //             value = 0.5*(tanh(BG_POTENTIAL_SLOPE*pot)/PI + 1.0);
            hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*value;
            hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb);
            if (SHADE_BG_COLOR_2D) 
                for (k=0; k<3; k++) rgb[k] *= wsphere[cell].cos_angle_pot;
            phi = XMIN + (double)i*dphi;
            theta = YMIN + (double)j*dtheta;
            glColor3f(rgb[0], rgb[1], rgb[2]);
            glVertex2d(phi, theta);
            glVertex2d(phi + dphi_res, theta);
            glVertex2d(phi + dphi_res, theta + dtheta_res);
            glVertex2d(phi, theta + dtheta_res);
        }
    glEnd();
 }
 void draw_particles(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], int plot, double beta, t_collision *collisions, int ncollisions, int bg_color, t_hashgrid hashgrid[HASHX*HASHY], t_lj_parameters params)
 {
@@ -8869,7 +8938,7 @@ void draw_particles(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRC
                    else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio;
                }
                else hue = 0.25*PARTICLE_HUE_MIN + 0.75*PARTICLE_HUE_MAX;
-                erase_area_hsl_turbo(0.0, YMIN, 2.0*PARTIAL_THERMO_WIDTH, PARTIAL_THERMO_HEIGHT*(YMAX - YMIN),  hue, 0.9, 0.15);
+                erase_area_hsl_turbo(0.5*(XMIN + XMAX), YMIN, 2.0*PARTIAL_THERMO_WIDTH, PARTIAL_THERMO_HEIGHT*(YMAX - YMIN),  hue, 0.9, 0.15);
                break;
            }
            case(TH_LAYER): 
@@ -10055,6 +10124,13 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
        sprintf(message, "B field = %.2f", 25.0*NVID*DT_PARTICLE*params.bfield);
        write_text(xtext + 0.08, y, message);
    }   
    if (CHANGE_OMEGA_SPHERE)    /* print angular frequency of sphere */
    {
        erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0);
        glColor3f(1.0, 1.0, 1.0);
        sprintf(message, "Omega = %.2f", params.omega_sphere);
        write_text(xtext + 0.08, y, message);
    }   
    if (PRINT_NPARTICLES)   /* print number of particles */
    {
@@ -10135,10 +10211,20 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
        erase_area_hsl(xmid, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0);
        glColor3f(1.0, 1.0, 1.0);
        sprintf(message, "Gravity %.2f", params.gravity/GRAVITY);
 //         write_text(xmidtext + 0.1, y, message);
        write_text(xmidtext, y, message);
        y -= 0.1;
    } 
    if (SPHERE_GRAVITY > 0)
    {
        erase_area_hsl(xmid, y + 0.025*scale, 0.22*scale, 0.05*scale, 0.0, 0.9, 0.0);
        glColor3f(1.0, 1.0, 1.0);
        if (SPHERE_GRAVITY == 1)
            sprintf(message, "g angle %.2f", params.g_angle*180.0/PI - 90.0);
        else sprintf(message, "g angle %.2f", params.g_angle*180.0/PI);
        write_text(xmidtext, y, message);
    } 
    if (CHANGE_RADIUS)
    {
        erase_area_hsl(xmid, y + 0.025*scale, 0.3*scale, 0.05*scale, 0.0, 0.9, 0.0);
@@ -11685,7 +11771,7 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES])
 /* initialize all particles, obstacles, and the hashgrid */
 {
    int i, j, k, n, type, nactive = 0, hashcell;
-    double x, y, h, xx, yy, rnd, angle, dx, dy;
+    double x, y, h, xx, yy, rnd, angle, dx, dy, mratio;
    printf("Initializing configuration\n");
@@ -11939,13 +12025,57 @@ t_segment segment[NMAXSEGMENTS], t_molecule molecule[NMAXCIRCLES])
 //     }
    /* position-dependent particle type */
-    if (POSITION_DEPENDENT_TYPE) for (i=0; i<ncircles; i++)
+    switch (POSITION_DEPENDENT_TYPE) {
-        if (((!POSITION_Y_DEPENDENCE)&&((particle[i].xc - POSITION_DEP_X)*POSITION_DEP_SIGN < 0.0))||((POSITION_Y_DEPENDENCE)&&(particle[i].yc*POSITION_DEP_SIGN < 0.0))) 
+        case (PDIC_TWO):
        {
-            particle[i].type = 2;
+            for (i=0; i<ncircles; i++)
-            particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
+            {
-            particle[i].radius = MU_B;
+                if (((!POSITION_Y_DEPENDENCE)&&((particle[i].xc - POSITION_DEP_X)*POSITION_DEP_SIGN < 0.0))||((POSITION_Y_DEPENDENCE)&&((particle[i].yc - POSITION_DEP_Y)*POSITION_DEP_SIGN < 0.0))) 
                {
                    particle[i].type = 2;
                    particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
                    particle[i].radius = MU_B;
                }
            }
            break;
        }
        case (PDIC_THREE):
        {
            for (i=0; i<ncircles; i++)
            {
                if (POSITION_Y_DEPENDENCE)
                {
                    if (particle[i].yc > POSITION_DEP_Y)
                    {
                        particle[i].type = 3;
                        particle[i].mass_inv = 1.0/PARTICLE_MASS_C;
                        particle[i].radius = MU_C;
                    }
                    else if (particle[i].yc > YMAX - POSITION_DEP_Y)
                    {
                        particle[i].type = 2;
                        particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
                        particle[i].radius = MU_B;
                    }
                }
                else
                {
                    /* TODO */
                }
            }
            break;
        }
        case (PDIC_MASSY):
        {
            for (i=0; i<ncircles; i++)
            {
                mratio = (particle[i].yc - YMIN)/(YMAX - YMIN);
                particle[i].mass_inv = 1.0/(PARTICLE_MASS*(1.0 + mratio*(POSITION_DEP_MASS_RATIO-1.0)));
            }
            break;
        }
    }
    /* add copies in case of particle pairing */
    if (PAIR_PARTICLES) init_particle_pairs(particle, molecule);
@@ -12386,7 +12516,7 @@ int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin, t
 /* only couple particles satisfying condition PARTIAL_THERMO_REGION to thermostat */
 {
    int condition, i, nthermo = 0;
-    double x, y, height, a, b;
+    double x, y, height, a, b, xmid;
    static double maxheight;
    static int first = 1;
@@ -12430,7 +12560,8 @@ int partial_thermostat_coupling(t_particle particle[NMAXCIRCLES], double xmin, t
            {
                x = particle[i].xc;
                y = particle[i].yc;
-                condition = ((y < YMIN + PARTIAL_THERMO_HEIGHT*(YMAX - YMIN))&&(vabs(x) < PARTIAL_THERMO_WIDTH*XMAX));
+                xmid = 0.5*(XMIN + XMAX);
                condition = ((y < YMIN + PARTIAL_THERMO_HEIGHT*(YMAX - YMIN))&&(vabs(x - xmid) < PARTIAL_THERMO_WIDTH*XMAX));
                break;
            }
            case (TH_LAYER):
@@ -16964,6 +17095,23 @@ double efield_schedule(int i)
 }
 double omega_sphere_schedule(int i)
 {
    static double omagafactor;
    static int first = 1;
    double omega;
    if (first) 
    {
        omagafactor = (OMEGA_SPHERE_FACTOR-1.0)/(double)(NSTEPS);
        first = 0;
    }
    if (i < INITIAL_TIME) omega = OMEGA_SPHERE;
    else omega = OMEGA_SPHERE*(1.0 + (double)(i-INITIAL_TIME)*omagafactor);
    return(omega);
 }
 void draw_current_plot(double *currents, int i)
 /* draw plot of current vs E field */
 {
@@ -17705,18 +17853,19 @@ int wrap_particles_sphere_bc(t_particle particle[NMAXCIRCLES], t_molecule molecu
    return(move);
 }
-int init_absorbers(t_absorber absorber[NMAX_ABSORBERS])
+int init_absorbers(t_absorber absorber[NMAX_ABSORBERS], int pattern, double radius)
 /* initialise absorbing discs */
 {
-    int k, nabsorb = 0;
+    int k, nabsorb = 0, ndiscs;
    double alpha, beta, gamma;
    t_point *point;
-    switch (ABSORBER_PATTERN) {
+    switch (pattern) {
        case (AP_ONE):
        {
            absorber[nabsorb].xc = ABSORBER_X;
            absorber[nabsorb].yc = ABSORBER_Y;
-            absorber[nabsorb].radius = ABSORBER_R;
+            absorber[nabsorb].radius = radius;
            nabsorb++;
            break;
        }
@@ -17724,11 +17873,23 @@ int init_absorbers(t_absorber absorber[NMAX_ABSORBERS])
        {
            absorber[nabsorb].xc = ABSORBER_X;
            absorber[nabsorb].yc = ABSORBER_Y;
-            absorber[nabsorb].radius = ABSORBER_R;
+            absorber[nabsorb].radius = radius;
            nabsorb++;
            absorber[nabsorb].xc = ABSORBER_X;
            absorber[nabsorb].yc = YMAX - ABSORBER_Y;
-            absorber[nabsorb].radius = ABSORBER_R;
+            absorber[nabsorb].radius = radius;
            nabsorb++;
            break;
        }
        case (AP_POLES):
        {
            absorber[nabsorb].xc = 0.0;
            absorber[nabsorb].yc = 0.0;
            absorber[nabsorb].radius = radius;
            nabsorb++;
            absorber[nabsorb].xc = 0.0;
            absorber[nabsorb].yc = PI;
            absorber[nabsorb].radius = radius;
            nabsorb++;
            break;
        }
@@ -17739,11 +17900,11 @@ int init_absorbers(t_absorber absorber[NMAX_ABSORBERS])
            {
                absorber[nabsorb].xc = (double)k*PID;
                absorber[nabsorb].yc = alpha;
-                absorber[nabsorb].radius = ABSORBER_R;
+                absorber[nabsorb].radius = radius;
                nabsorb++;
                absorber[nabsorb].xc = (double)k*PID;
                absorber[nabsorb].yc = PI - alpha;
-                absorber[nabsorb].radius = ABSORBER_R;
+                absorber[nabsorb].radius = radius;
                nabsorb++;
            }
            break;
@@ -17798,7 +17959,7 @@ int init_absorbers(t_absorber absorber[NMAX_ABSORBERS])
            for (k=0; k<nabsorb; k++) 
            {
-                absorber[k].radius = ABSORBER_R;
+                absorber[k].radius = radius;
                absorber[k].xc += ABSORBER_X;
                if (absorber[k].xc > DPI) absorber[k].xc -= DPI;
            }
@@ -17832,12 +17993,31 @@ int init_absorbers(t_absorber absorber[NMAX_ABSORBERS])
            for (k=0; k<nabsorb; k++) 
            {
-                absorber[k].radius = ABSORBER_R;
+                absorber[k].radius = radius;
                absorber[k].xc += ABSORBER_X;
                if (absorber[k].xc > DPI) absorber[k].xc -= DPI;
            }
            break;
        }
        case (AP_PDISC):
        {
            point = (t_point *)malloc(NMAX_ABSORBERS*sizeof(t_point));
            ndiscs = generate_poisson_discs(point, NMAX_ABSORBERS, PI+0.1, DPI-0.1, 0.1, PI-0.1, ABSORBER_PDIST, 1);
            for (k=0; k<ndiscs; k++) 
            {
                if ((point[k].xc > PI)&&(point[k].xc < DPI))
                {
                    absorber[nabsorb].xc = point[k].xc;
                    absorber[nabsorb].yc = point[k].yc;
                    absorber[nabsorb].radius = radius;
                    nabsorb++;
                }
            }
            free(point);
            break;
        }
    }
    return(nabsorb);
 }
@@ -17871,7 +18051,7 @@ void draw_frame_2d(int i, int plot, int bg_color, int ncollisions, int traj_posi
                t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], 
                t_group_data *group_speeds, t_group_segments *segment_group, 
                t_belt *conveyor_belt, int *tracer_n,
-                t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES], t_absorber absorber[NMAX_ABSORBERS])
+                t_otriangle otriangle[NMAX_TRIANGLES_PER_OBSTACLE*NMAXOBSTACLES], t_ofacet ofacet[NMAXOBSTACLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], t_absorber absorber[NMAX_ABSORBERS])
 /* draw a movie frame */
 {
    printf("Drawing frame\n");
@@ -17879,8 +18059,9 @@ void draw_frame_2d(int i, int plot, int bg_color, int ncollisions, int traj_posi
    compute_all_particle_colors(particle, cluster, plot);
    if ((COLOR_BACKGROUND)&&(bg_color > 0)) 
    {
-        compute_background_color(particle, segment, obstacle, bg_color, hashgrid);
+        compute_background_color(particle, segment, obstacle, bg_color, hashgrid, wsphere);
-        draw_background(particle, segment, obstacle, bg_color, hashgrid);
+        if (bg_color == BG_POTENTIAL) draw_potential(wsphere);
        else draw_background(particle, segment, obstacle, bg_color, hashgrid);
    }
 //     else if (!TRACER_PARTICLE) blank();
@@ -17934,7 +18115,7 @@ void draw_frame_3d(int i, int plot, int bg_color, int ncollisions, int traj_posi
    /* initialize background colors */
    if ((COLOR_BACKGROUND)&&(bg_color > 0))
-        compute_background_color(particle, segment, obstacle, bg_color, hashgrid);
+        compute_background_color(particle, segment, obstacle, bg_color, hashgrid, wsphere);
    /* compute all particle colors */
    compute_all_particle_colors(particle, cluster, plot);
@@ -17970,5 +18151,5 @@ void draw_frame(int i, int plot, int bg_color, int ncollisions, int traj_positio
    if (DRAW_SPHERE) draw_frame_3d(i, plot, bg_color, ncollisions, traj_position,
            traj_length, wall, pressure, pleft, pright, currents, particle_numbers, refresh, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere, absorber);
    else draw_frame_2d(i, plot, bg_color, ncollisions, traj_position, traj_length, 
-            wall, pressure, pleft, pright, currents, particle_numbers, refresh, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, absorber);
+            wall, pressure, pleft, pright, currents, particle_numbers, refresh, params, particle, cluster, collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt, tracer_n, otriangle, ofacet, wsphere, absorber);
 }
@@ -690,11 +690,11 @@ void compute_all_particle_colors(t_particle particle[NMAXCIRCLES], t_cluster clu
 }
-void compute_background_color(t_particle particle[NMAXCIRCLES], t_segment segment[NMAXSEGMENTS], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY])
+void compute_background_color(t_particle particle[NMAXCIRCLES], t_segment segment[NMAXSEGMENTS], t_obstacle obstacle[NMAXOBSTACLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
 /* color background according to particle properties */
 {
-    int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs;
+    int i, j, k, n, p, q, m, nnb, number, avrg_fact, obs, i0, j0, cell, i1, j1, cell1;
-    double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue, valx, valy, lum;
+    double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue, valx, valy, lum, nablax, nablay, norm, pscal, ca;
    static int first = 1, counter = 0;
    static double area_factor;
@@ -1060,6 +1060,22 @@ void compute_background_color(t_particle particle[NMAXCIRCLES], t_segment segmen
                    hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIRECTION);
                    break;
                }
                case (BG_POTENTIAL):
                {
                    i0 = (int)(0.5*(hashgrid[n].x1+hashgrid[n].x2)*(double)NX_SPHERE/DPI);
                    j0 = (int)(0.5*(hashgrid[n].y1+hashgrid[n].y2)*(double)NY_SPHERE/PI);
                    cell = i0*NY_SPHERE+j0;
                    if (cell < 0) cell = 0;
                    value = atan(BG_POTENTIAL_SLOPE*wsphere[cell].potential)/PI + 0.5;
                    hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*value;
                    hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb);
                    if (SHADE_BG_COLOR_2D)
                        for (k=0; k<3; k++) rgb[k] *= wsphere[cell].cos_angle_pot;
                    break;
                }
            }
            if (DOUBLE_MOVIE)
@@ -1286,6 +1302,7 @@ void init_lj_sphere(t_lj_sphere *wsphere)
            wsphere[i*NY_SPHERE+j].z = -wsphere[i*NY_SPHERE+j].ctheta;
            wsphere[i*NY_SPHERE+j].radius = 1.0;
            wsphere[i*NY_SPHERE+j].initial_radius = 1.0;
            wsphere[i*NY_SPHERE+j].cos_angle_sphere = wsphere[i*NY_SPHERE+j].x*light[0] + wsphere[i*NY_SPHERE+j].y*light[1] + wsphere[i*NY_SPHERE+j].z*light[2];
@@ -1361,7 +1378,7 @@ double test_distance(double phi, double psi, int i, t_particle* particle)
 }
-void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRCLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
+void add_particle_to_sphere(int i, int j, int i0, int j0, int part, t_particle particle[NMAXCIRCLES], t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
 /* compute the effect at point (i,j) of adding a particle */
 {
    int draw_normal_particle = 1, k;
@@ -1554,7 +1571,8 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC
        if (dist < r)
        {
-            h = 1.0 + 1.5*sqrt(r*r - dist*dist);
+//             h = wsphere[i*NY_SPHERE+j].initial_radius + 1.5*sqrt(r*r - dist*dist);
            h = wsphere[i0*NY_SPHERE+j0].initial_radius + 1.5*sqrt(r*r - dist*dist);
            if ((h > wsphere[i*NY_SPHERE+j].radius)&&(!wsphere[i*NY_SPHERE+j].locked))
            {
                wsphere[i*NY_SPHERE+j].r = particle[part].rgb[0];
@@ -1564,6 +1582,7 @@ void add_particle_to_sphere(int i, int j, int part, t_particle particle[NMAXCIRC
            }
        }
    }
 }
 void draw_trajectory_sphere(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_hashgrid hashgrid[HASHX*HASHY], int traj_position, int traj_length, t_particle *particle, t_cluster *cluster, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot)
@@ -1758,7 +1777,7 @@ void draw_segments_sphere(t_segment segment[NMAXSEGMENTS], t_lj_sphere wsphere[N
 /* draw the repelling segments on the sphere */
 {
    int s, i, cell, npoints, i0, j0, i1, j1, p, q;
-    double x1, y1, x2, y2, x, y, dt, length;
+    double x1, y1, x2, y2, x, y, dt, length, ca;
    static double dphi, dtheta;
    static int first = 1;
@@ -1796,11 +1815,19 @@ void draw_segments_sphere(t_segment segment[NMAXSEGMENTS], t_lj_sphere wsphere[N
                    if (j1 < 0) j1 = 0;
                    if (j1 >= NY_SPHERE) j1 = NY_SPHERE-1;
                    cell = i1*NY_SPHERE+j1;
                    wsphere[cell].r = 0.0;
                    wsphere[cell].g = 0.0;
                    wsphere[cell].b = 0.0;
-                    wsphere[cell].radius += 0.005;
+                    if (wsphere[cell].locked == 0)
                    {
                        ca = wsphere[cell].cos_angle_sphere;
                        ca = (ca + 1.0)*0.4 + 0.2;
                        wsphere[cell].r = ca;
                        wsphere[cell].g = ca;
                        wsphere[cell].b = ca;
                        wsphere[cell].radius += 0.05;
                        wsphere[cell].locked = 1;
                    }
                }
        }
    }
@@ -1840,11 +1867,131 @@ void draw_absorbers_sphere(t_absorber absorber[NMAX_ABSORBERS], t_lj_sphere wsph
                {
                    cell = i*NY_SPHERE+j;
                    wsphere[cell].locked = 1;
-                    wsphere[cell].radius = 0.95;
+//                     wsphere[cell].radius = 0.95;
                    wsphere[cell].radius = wsphere[cell].initial_radius - 0.05;
                }
            }
 }
 void compute_gradient_potential_sphere(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE])
 /* compute the gradient of the potential on the sphere */
 {
    int i, j, cell;
    double dphi_inv, dtheta_inv, stheta_inv;
    dphi_inv = 0.5*(double)NX_SPHERE/DPI;
    dtheta_inv = 0.5*(double)NY_SPHERE/PI;
    /* bulk */
    for (j=1; j<NX_SPHERE-1; j++)
    {
        stheta_inv = wsphere[j].stheta;
        for (i=1; i<NX_SPHERE-1; i++)
        {
            cell = i*NY_SPHERE+j;
            wsphere[cell].nablax = (wsphere[cell+NY_SPHERE].potential - wsphere[cell-NY_SPHERE].potential)*dphi_inv*stheta_inv;
            wsphere[cell].nablay = (wsphere[cell+1].potential - wsphere[cell-1].potential)*dtheta_inv;
        }
        /* i=0 */
        cell = j;
        wsphere[cell].nablax = (wsphere[cell+NY_SPHERE].potential - wsphere[(NX_SPHERE-1)+j].potential)*dphi_inv*stheta_inv;
        wsphere[cell].nablay = (wsphere[cell+1].potential - wsphere[cell-1].potential)*dtheta_inv;
        /* i=NX_SPHERE-1 */
        cell = (NX_SPHERE-1)*NY_SPHERE+j;
        wsphere[cell].nablax = (wsphere[j].potential - wsphere[cell-NY_SPHERE].potential)*dphi_inv*stheta_inv;
        wsphere[cell].nablay = (wsphere[cell+1].potential - wsphere[cell-1].potential)*dtheta_inv;
    }
    /* poles */
    stheta_inv = wsphere[1].stheta;
    for (i=0; i<NX_SPHERE; i++)
    {
        /* j=0 */
        cell = i*NY_SPHERE;
        wsphere[cell].nablax = 0.0;
        wsphere[cell].nablay = 2.0*(wsphere[cell+1].potential - wsphere[cell].potential)*dtheta_inv;
        /* j=NY_SPHERE-1 */
        cell = i*NY_SPHERE + NY_SPHERE - 1;
        wsphere[cell].nablax = 0.0;
        wsphere[cell].nablay = 2.0*(wsphere[cell].potential - wsphere[cell-1].potential)*dtheta_inv;
    }
 }
 void init_potential_sphere(t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], t_absorber wells[NMAX_ABSORBERS])
 /* initialise the potential and its gradient, for option ADD_POTENTIAL_SPHERE */
 {
    int i, j;
    double phi, theta, phi1, theta1, dist, pot, x, y, z, max, nablax, nablay, norm, pscal, ca, shade_scale2;
    printf("Initialising potential on sphere\n");
    for (i=0; i<NX_SPHERE*NY_SPHERE; i++) wsphere[i].potential = 0.0; 
    switch(SPHERE_POTENTIAL) {
        case (SPP_WELLS):
        {
            for (i=0; i<NX_SPHERE*NY_SPHERE; i++) 
            {
                phi = wsphere[i].phi;
                theta = wsphere[i].theta;
                for (j=0; j<nwells; j++)
                {
                    phi1 = wells[j].xc;
                    theta1 = wells[j].yc;
                    dist = dist_sphere(phi, theta, phi1, theta1)/wells[j].radius;
                    pot = -POT_SPHERE_AMP/(POT_SPHERE_SMOOTH + dist*dist);
                    wsphere[i].potential += pot;
                }
            }
            break;
        }
        case (SPP_CUBE):
        {
            for (i=0; i<NX_SPHERE*NY_SPHERE; i++) 
            {
                max = 0.0;
                x = vabs(wsphere[i].x);
                y = vabs(wsphere[i].y);
                z = vabs(wsphere[i].z);
                if (x > max) max = x;
                if (y > max) max = y;
                if (z > max) max = z;
                wsphere[i].potential = 1.0/max - 1.0; 
            }
            break;
        }
    }
    for (i=0; i<NX_SPHERE*NY_SPHERE; i++) 
    {
        wsphere[i].initial_radius = 1.0 + RSCALE_POTENTIAL*wsphere[i].potential;
        if (wsphere[i].initial_radius < 0.5) wsphere[i].initial_radius = 0.5;
    }
    compute_gradient_potential_sphere(wsphere);
    if (SHADE_BG_COLOR_2D)
    {
        shade_scale2 = SHADE_SCALE_BG_2D*SHADE_SCALE_BG_2D;
        for (i=0; i<NX_SPHERE*NY_SPHERE; i++) 
        {
            nablax = wsphere[i].nablax;
            nablay = wsphere[i].nablay;
            norm = sqrt(shade_scale2 + nablax*nablax + nablay*nablay);
            pscal = -nablax*light[0] - nablay*light[1] + SHADE_SCALE_BG_2D;
            ca = pscal/norm;
            wsphere[i].cos_angle_pot = (ca + 1.0)*0.4 + 0.2;
 //             printf("cos_angle_pot[%i] = %.3lg\n", i, wsphere[i].cos_angle_pot);
        }
    }
 }
 void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_cluster cluster[NMAXCIRCLES], t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES], t_segment segment[NMAXSEGMENTS], int traj_position, int traj_length, t_lj_sphere wsphere[NX_SPHERE*NY_SPHERE], int *tracer_n, int plot, int bg_color, t_absorber absorber[NMAX_ABSORBERS])
 /* initialize sphere radius and colors from particles */
 {
@@ -1864,7 +2011,7 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
    #pragma omp parallel for private(i)
    for (i=0; i<NX_SPHERE*NY_SPHERE; i++) if (!wsphere[i].locked)
    {
-        wsphere[i].radius = 1.0;
+        wsphere[i].radius = wsphere[i].initial_radius;
        if ((COLOR_BACKGROUND)&&(bg_color > 0))
        {
@@ -1885,8 +2032,8 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
    if (TRACER_PARTICLE) 
        draw_trajectory_sphere(trajectory, hashgrid, traj_position, traj_length, particle, cluster, wsphere, tracer_n, plot);
-    if (ADD_FIXED_SEGMENTS)
+//     if (ADD_FIXED_SEGMENTS)
-        draw_segments_sphere(segment, wsphere);
+//         draw_segments_sphere(segment, wsphere);
 //     if (ADD_ABSORBERS)
 //         draw_absorbers_sphere(absorber, wsphere); 
@@ -1910,13 +2057,14 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
        if (particle[part].active)
        {
            deltaj = (int)(2.5*particle[part].radius/dtheta);
            i0 = 0;
            j0 = (int)(particle[part].yc/dtheta);
            jmax = j0 + deltaj;
            #pragma omp parallel for private(i,j,x,y,dist,r)
            for (j=0; j<jmax; j++)
                for (i=0; i<NX_SPHERE; i++)
-                    add_particle_to_sphere(i, j, part, particle, wsphere);
+                    add_particle_to_sphere(i, j, i0, j0, part, particle, wsphere);
        }
    }
@@ -1929,13 +2077,14 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
        if (particle[part].active)
        {
            deltaj = (int)(2.5*particle[part].radius/dtheta);
            i0 = 0;
            j0 = (int)(particle[part].yc/dtheta);
            jmin = j0 - deltaj;
            #pragma omp parallel for private(i,j,x,y,dist,r)
            for (j=jmin; j<NY_SPHERE; j++)
                for (i=0; i<NX_SPHERE; i++)
-                    add_particle_to_sphere(i, j, part, particle, wsphere);
+                    add_particle_to_sphere(i, j, i0, j0, part, particle, wsphere);
        }
    }
@@ -1971,7 +2120,7 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
                    for (j=jmin; j<jmax; j++)
                        for (i=imin; i<imax; i++)
                        {
-                            add_particle_to_sphere(i, j, part, particle, wsphere);
+                            add_particle_to_sphere(i, j, i0, j0, part, particle, wsphere);
                        }
                }
            }
@@ -2002,9 +2151,9 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
                for (j=jmin; j<jmax; j++)
                {
                    for (i=0; i<imax; i++)
-                        add_particle_to_sphere(i, j, part, particle, wsphere);
+                        add_particle_to_sphere(i, j, i0, j0, part, particle, wsphere);
                    for (i=imin; i<NX_SPHERE; i++)
-                        add_particle_to_sphere(i, j, part, particle, wsphere);
+                        add_particle_to_sphere(i, j, i0, j0, part, particle, wsphere);
                }
            }
        }
@@ -2034,9 +2183,9 @@ void init_sphere_radius(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
                for (j=jmin; j<jmax; j++)
                {
                    for (i=imin; i<NX_SPHERE; i++)
-                        add_particle_to_sphere(i, j, part, particle, wsphere);
+                        add_particle_to_sphere(i, j, i0, j0, part, particle, wsphere);
                    for (i=0; i<imax; i++)
-                        add_particle_to_sphere(i, j, part, particle, wsphere);
+                        add_particle_to_sphere(i, j, i0, j0, part, particle, wsphere);
                }
            }
        }
@@ -4,6 +4,8 @@ double intstep;   /* integration step */
 double intstep1;  /* integration step used in absorbing boundary conditions */
 double viscosity; /* viscosity (parameter in front of Laplacian) */
 double rpslzb = RPSLZB;    /* second parameter in Rock-Paper-Scissors-Lizard-Spock equation */
 double a_gs = A_GS;  /* first parameter in Gray-Scott model */
 double b_gs = B_GS;  /* second parameter in Gray-Scott model */
 double flow_speed;  /* flow speed for laminar boundary conditions in Euler equation */
 // double gaussian()
@@ -33,6 +35,7 @@ double flow_speed;  /* flow speed for laminar boundary conditions in Euler equat
 //     return X;
 // }
 void init_random(double mean, double amplitude, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
 /* initialise field with gaussian at position (x,y) */
 {
@@ -216,6 +219,48 @@ void add_random_onefield_smoothed(int field, double mean, double amplitude, doub
    printf("Fields initialised\n");
 }
 void init_random_zerosum(double mean, double amplitude, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
 /* initialise two fields randomly with u + v = 1 */
 {
    int i, j, k, in;
    double xy[2], dist2, module, phase, scale2;    
    printf("Initialising xy_in\n");
    #pragma omp parallel for private(i,j)
    for (i=0; i<NX; i++)
    {
        if (i%100 == 0) printf("Column %i of %i\n", i, NX);
        for (j=0; j<NY; j++)
        {
            ij_to_xy(i, j, xy);
            if (SPHERE) xy_in[i*NY+j] = xy_in_billiard_sphere(i, j, wsphere);
            else xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
 //             xy_in[i*NY+j] = 1;
        }
    }
    printf("Initialising fields\n");
    for (i=0; i<NX; i++)
    {
        if (i%100 == 0) printf("Field %i of %i, column %i of %i\n", k, NFIELDS, i, NX);
        for (j=0; j<NY; j++)
        {
            if (xy_in[i*NY+j]) 
            {
                phi[0][i*NY+j] = mean + amplitude*(2.0*(double)rand()/RAND_MAX - 1.0);
                phi[1][i*NY+j] = 1.0 - phi[0][i*NY+j];
            }
            else 
            {
                phi[0][i*NY+j] = 0.0;
                phi[1][i*NY+j] = 1.0;
            }
        }
    }
    printf("Fields initialised\n");
 }
 void init_gaussian(double x, double y, double mean, double amplitude, double scalex, 
                   double *phi[NX], short int xy_in[NX*NY])
 /* initialise field with gaussian at position (x,y) */
@@ -1843,6 +1888,329 @@ void init_tidal_state(int add, double amp, double min, double *phi[NFIELDS], sho
 }
 void init_gs_blob(double x, double y, double mean, double amplitude, double scalex, 
                   double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
 /* initialise field with gaussian at position (x,y) */
 {
    int i, j, k;
    double xy[2], dist2, module, phase, scale2;    
    scale2 = scalex*scalex;
    printf("Initialising field\n");
    for (i=0; i<NX; i++)
    {
        if (i%100 == 0) printf("Field %i of %i, column %i of %i\n", k, NFIELDS, i, NX);
        for (j=0; j<NY; j++)
        {
            ij_to_xy(i, j, xy);
 	    xy_in[i*NY+j] = 1;
            dist2 = 0.5*(xy[0]-x)*(xy[0]-x) + (xy[1]-y)*(xy[1]-y);
            if (dist2 < scale2) module = 1.0;
            else module = 0.0;
            phi[0][i*NY+j] = module;
            phi[1][i*NY+j] = 1.0 - 0.5*module;
        }
    }
 }
 void init_gs_blobs_hex(int nx, int ny, double amplitude, double scalex, double scaley, 
                   double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
 /* initialise Gray-Scott model with grid of ellipses */
 {
    int i, j, p, q;
    double x, y, xy[2], distmin, dist2, module, scalex2, scaley2, dx, dy;    
    scalex2 = scalex*scalex;
    scaley2 = scaley*scaley;
    dx = (XMAX - XMIN)/(double)nx;
    dy = (YMAX - YMIN)/(double)ny;
    printf("Initialising field\n");
    for (i=0; i<NX; i++)
    {
        if (i%100 == 0) printf("Fields, column %i of %i\n", i, NX);
        for (j=0; j<NY; j++)
        {
            ij_to_xy(i, j, xy);
 	    xy_in[i*NY+j] = 1;
            distmin = 10.0;
            for (p=0; p<nx+2; p++)
                for (q=0; q<ny+2; q++)
                {
                    x = XMIN + ((double)p - 0.5)*dx;
                    y = YMIN + ((double)q - 0.5)*dy;
                    if (q%2 == 1) x += 0.5*dx;
                    dist2 = (xy[0]-x)*(xy[0]-x)/scalex + (xy[1]-y)*(xy[1]-y)/scaley;
                    if (dist2 < distmin) distmin = dist2;
                }
            if (distmin < 1.0) module = 1.0;
            else module = 0.0;
            phi[0][i*NY+j] = module;
            phi[1][i*NY+j] = 1.0 - 0.5*module;
        }
    }
 }
 void init_gs_blobs_square(int nx, int ny, double amplitude, double scalex, double scaley, 
                   double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
 /* initialise Gray-Scott model with grid of ellipses */
 {
    int i, j, p, q;
    double x, y, xy[2], distmin, dist2, module, scalex2, scaley2, dx, dy;    
    scalex2 = scalex*scalex;
    scaley2 = scaley*scaley;
    dx = (XMAX - XMIN)/(double)nx;
    dy = (YMAX - YMIN)/(double)ny;
    printf("Initialising field\n");
    for (i=0; i<NX; i++)
    {
        if (i%100 == 0) printf("Fields, column %i of %i\n", i, NX);
        for (j=0; j<NY; j++)
        {
            ij_to_xy(i, j, xy);
 	    xy_in[i*NY+j] = 1;
            distmin = 10.0;
            for (p=0; p<nx+2; p++)
                for (q=0; q<ny+2; q++)
                {
                    x = XMIN + ((double)p - 0.5)*dx;
                    y = YMIN + ((double)q - 0.5)*dy;
                    dist2 = (xy[0]-x)*(xy[0]-x)/scalex + (xy[1]-y)*(xy[1]-y)/scaley;
                    if (dist2 < distmin) distmin = dist2;
                }
            if (distmin < 1.0) module = 1.0;
            else module = 0.0;
            phi[0][i*NY+j] = module;
            phi[1][i*NY+j] = 1.0 - 0.5*module;
        }
    }
 }
 int generate_poisson_sample(double xmin, double xmax, double ymin, double ymax, double dpoisson, int nmaxpoints, double *px, double *py)
 /* generate a Poisson disc sample in a given rectangle */
 {
    int i, j, k, n_p_active, ncandidates=500, naccepted, p, q, npoints; 
    double r, phi, x, y, x1, y1, distance, dist1;
    short int active_poisson[NX*NY], far;
    printf("Generating Poisson disc sample\n");
    /* generate first point */
    px[0] = xmin + (xmax - xmin)*(double)rand()/RAND_MAX;
    py[0] = ymin + (ymax - ymin)*(double)rand()/RAND_MAX;
    active_poisson[0] = 1;
    n_p_active = 1;
    npoints = 1;
    while ((n_p_active > 0)&&(npoints < nmaxpoints))
    {
        /* randomly select an active circle */
        i = rand()%(npoints);
        while (!active_poisson[i]) i = rand()%(npoints);                 
        /* generate new candidates */
        naccepted = 0;
        for (j=0; j<ncandidates; j++)
        {
            r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0);
            phi = DPI*(double)rand()/RAND_MAX;
            x = px[i] + r*cos(phi);
            y = py[i] + r*sin(phi);
            far = 1;
            dist1 = 10.0;
            for (k=0; k<npoints; k++) /*if ((k!=i))*/
            {
                for (p=-1; p<2; p++)
                    for (q=-1; q<2; q++)
                    {
                        /* new point is far away from point k */
                        distance = module2(x - px[k] - (double)p*(XMAX-XMIN), y - py[k] - (double)q*(YMAX-YMIN));
                        far = far*(distance >= dpoisson);
                        if (distance < dist1) dist1 = distance;
                    }
            }
            if ((far)&&(npoints < nmaxpoints))    /* accept new point */
            {
                printf("New grid point at (%.3f,%.3f) accepted\n", x, y);
                printf("Minimal distance = %.3lg\n", dist1); 
                printf("%i grid points generated\n", npoints);
                px[npoints] = x;
                py[npoints] = y;
                active_poisson[npoints] = 1;
                npoints++;
                n_p_active++;
                naccepted++;
            }
        }
        if (naccepted == 0)    /* inactivate circle i */ 
        {
            active_poisson[i] = 0;
            n_p_active--;
        }
        printf("%i active points\n", n_p_active);
    }
    printf("Generated %i points\n", npoints);
    return(npoints);
 }
 void init_gs_blobs_poisson(double dpoisson, double nmax, double amplitude, double scalex, double scaley, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
 /* initialise Gray-Scott model with Poisson disc sample of ellipses */
 {
    int i, j, k, npoints, p, q;
    double x, y, xy[2], distmin, dist2, module, scalex2, scaley2, dx, dy;  
    double *px, *py;
    px = (double *)malloc(nmax*sizeof(double));
    py = (double *)malloc(nmax*sizeof(double));
    scalex2 = scalex*scalex;
    scaley2 = scaley*scaley;
    printf("Initialising field\n");
    npoints = generate_poisson_sample(XMIN, XMAX, YMIN, YMAX, dpoisson, nmax, px, py);
    for (i=0; i<NX; i++)
    {
        if (i%100 == 0) printf("Fields, column %i of %i\n", i, NX);
        for (j=0; j<NY; j++)
        {
            ij_to_xy(i, j, xy);
 	    xy_in[i*NY+j] = 1;
            distmin = 10.0;
            for (k = 0; k < npoints; k++)
            {
                for (p=-1; p<2; p++)
                    for (q=-1; q<2; q++)
                    {
                        x = px[k] + (double)p*(XMAX - XMIN);
                        y = py[k] + (double)q*(YMAX - YMIN);
                        dist2 = (xy[0]-x)*(xy[0]-x)/scalex + (xy[1]-y)*(xy[1]-y)/scaley;
                        if (dist2 < distmin) distmin = dist2;
                    }
            }
            if (distmin < 1.0) module = 1.0;
            else module = 0.0;
            phi[0][i*NY+j] = module;
            phi[1][i*NY+j] = 1.0 - 0.5*module;
        }
    }
    free(px);
    free(py);
 }
 void init_gs_blobs_sphere(double amplitude, double scalex, double scaley, 
                   double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY], int pattern)
 /* initialise Gray-Scott model with grid of ellipses */
 {
    int i, j, p, q, nircles;
    double x, y, xy[2], distmin, dist2, module, scalex2, scaley2, dx, dy, phix, thetay;    
    t_circles_sphere circles[NMAXCIRC_SPHERE];
 //     scalex2 = scalex*scalex;
 //     scaley2 = scaley*scaley;
    printf("Initialising field\n");
    for (i=0; i<NX*NY; i++)
        xy_in[i] = 1;
    nircles = init_circle_sphere(circles, pattern);
    for (i=0; i<ngridpoints; i++)
    {
        if (i%10000 == 0) printf("Fields, point %i of %i\n", i, ngridpoints);
 //         xy_in[i] = 1;
        phix = wsphere[i].phigrid;
        thetay = wsphere[i].thetagrid;
        distmin = 10.0;
        for (p=0; p<nircles; p++)
        {
            x = circles[p].phi;
            y = circles[p].theta;
            dist2 = dist_sphere(phix, thetay, x, y)/scalex;
            if (dist2 < distmin) distmin = dist2;
        }
        if (distmin < 1.0) module = 1.0;
        else module = 0.0;
        phi[0][i] = 0.0;
        phi[1][i] = module*amplitude;
        phi[2][i] = 1.0 - 0.5*module*amplitude;
    }
 }
 void init_gs_rect_sphere(double amplitude, double phimin, double phimax, double thetamin, double thetamax, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY], short int add)
 /* initialise Gray-Scott model with grid of ellipses */
 {
    int i, j, p, q, nircles;
    double x, y, xy[2], distmin, dist2, module, scalex2, scaley2, dx, dy, phix, thetay;    
    printf("Initialising field\n");
    for (i=0; i<NX*NY; i++)
        xy_in[i] = 1;
    if (add) for (i=0; i<ngridpoints; i++)
    {
        if (i%10000 == 0) printf("Fields, point %i of %i\n", i, ngridpoints);
        phix = wsphere[i].phigrid;
        thetay = wsphere[i].thetagrid;
        if ((phix > phimin)&&(phix < phimax)&&(thetay > thetamin)&&(thetay < thetamax)) 
        {
            phi[0][i] = 0.0;
            phi[1][i] = amplitude;
            phi[2][i] = 1.0 - 0.5*amplitude;
        }
    }
    else for (i=0; i<ngridpoints; i++)
    {
        if (i%10000 == 0) printf("Fields, point %i of %i\n", i, ngridpoints);
        phix = wsphere[i].phigrid;
        thetay = wsphere[i].thetagrid;
        if ((phix > phimin)&&(phix < phimax)&&(thetay > thetamin)&&(thetay < thetamax)) module = 1.0;
        else module = 0.0;
        {
            phi[0][i] = 0.0;
            phi[1][i] = module*amplitude;
            phi[2][i] = 1.0 - 0.5*module*amplitude;
        }
    }
 }
 double distance_to_segment(double x, double y, double x1, double y1, double x2, double y2)
 /* distance of (x,y) to segment from (x1,y1) to (x2,y2) */
 {
@@ -4235,6 +4603,27 @@ void adjust_height(double *phi[NFIELDS], t_rde rde[NX*NY])
        }
 }
 void update_time_slice(int time, double *phi[NFIELDS], t_rde rde[NX*NY])
 /* adjust height of field */
 {
    int i, j; 
    static int first = 1;
    if (first)
    {
        #pragma omp parallel for private(i)
        for (i=0; i<NX*NY; i++) rde[i].time = 1.0;
        first = 0;
    }
    #pragma omp parallel for private(i,j)
    for (i=0; i<NX*NY; i++)
    {
        if (phi[0][i] < TIME_SLICE_THRESHOLD)
            rde[i].time = 1.0 - (double)time/(double)NSTEPS;
    }
 }
 void compute_direction(double *phi[NFIELDS], t_rde rde[NX*NY])
 /* compute the direction of a field */
 {
@@ -5958,6 +6347,12 @@ void compute_laplacian_grid(double phi_in[NX*NY], double phi_out[NX*NY], t_wave_
            phi_out[i] += phi_in[n];
        }
    }
    if (WEIGHT_GRID)
    {
        #pragma omp parallel for private(i)
        for (i=0; i<ngridpoints; i++) phi_out[i] *= wsphere[i].laplace_weight;
    }
 }
@@ -6781,7 +7176,11 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3]
            hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette);
            break;
        }
-        
+        case (Z_TIME_SLICE):
        {
            color_scheme_asym_palette(COLOR_SCHEME, palette, VSCALE_TIMESLICE*value, 1.0, 0, rgb);
            break;
        }
    }
 }
@@ -6923,7 +7322,7 @@ double compute_depth_colors_rde(int i, int j, t_rde rde[NX*NY], double potential
 }
-void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, int cplot, t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
+void compute_rde_fields(int time, double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, int cplot, t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
 /* compute the necessary auxiliary fields */
 {
    int i, j;
@@ -7043,6 +7442,18 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot,
                compute_direction(phi, rde);
            break;
        }
        case (E_GRAY_SCOTT):
        {
            if (zplot == Z_TIME_SLICE)
                update_time_slice(time, phi, rde); 
            break;
        }
        case (E_GRAY_SCOTT_SPHERE):
        {
            if (zplot == Z_TIME_SLICE)
                update_time_slice(time, phi, rde); 
            break;
        }
 //         case (E_KELLER_SEGEl):
 //         {
 //             
@@ -7257,7 +7668,14 @@ void init_zfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot, t_
            for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_zfield[movie] = &rde[i*NY+j].field_norm;
            break;
        }
-        
+        case (Z_TIME_SLICE):
        {
 //             #pragma omp parallel for private(i,j)
 //             for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_zfield[movie] = &rde[i*NY+j].time;
            #pragma omp parallel for private(i)
            for (i=0; i<NX*NY; i++) rde[i].p_zfield[movie] = &rde[i].time;
            break;
        }
    }
 }
@@ -7466,6 +7884,14 @@ void init_cfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int cplot, t_
            for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_cfield[movie] = &rde[i*NY+j].field_arg;
            break;
        }
        case (Z_TIME_SLICE):
        {
 //             #pragma omp parallel for private(i,j)
 //             for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_cfield[movie] = &rde[i*NY+j].time;
            #pragma omp parallel for private(i)
            for (i=0; i<NX*NY; i++) rde[i].p_cfield[movie] = &rde[i].time;
            break;
        }
    }
 }
@@ -8390,7 +8816,7 @@ void draw_periodicised_wave_3d(int movie, double *phi[NFIELDS], short int xy_in[
 }
-void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rde rde[NX*NY], 
+void draw_wave_rde(int time, 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], 
                   double potential[NX*NY], int zplot, int cplot, int palette, int fade, 
                   double fade_value, int refresh)
@@ -8404,7 +8830,7 @@ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rd
    if (refresh)
    {
        printf("Computing fields\n");
-        compute_rde_fields(phi, xy_in, zplot, cplot, rde, wsphere);      
+        compute_rde_fields(time, phi, xy_in, zplot, cplot, rde, wsphere);      
        printf("Computed fields\n");
        if ((PLOT_3D)&&(SHADE_3D)) 
        {
@@ -8831,17 +9257,6 @@ void block_poles(double phi[NX*NY])
 }
 double dist_sphere(double phi1, double theta1, double phi2, double theta2)
 /* returns angle between points given in spherical coordinates */
 {
    double dist;
    dist = sin(theta1)*sin(theta2)*cos(phi1 - phi2);
    dist += cos(theta1)*cos(theta2);
    return(acos(dist));
 }
 int generate_poisson_grid_sphere(t_wave_sphere wsphere[NX*NY], double dpoisson)
 /* generate a Poisson disc sample in a given rectangle */
 {
@@ -8914,6 +9329,7 @@ int generate_poisson_grid_sphere(t_wave_sphere wsphere[NX*NY], double dpoisson)
    return(npoints);
 }
 int generate_equilatitude_grid_sphere(t_wave_sphere wsphere[NX*NY])
 /* generate simulation gird with equally spaced longitudes and latitudes */
 {
@@ -9506,6 +9922,42 @@ void optimize_cubic_grid_sphere(t_wave_sphere wsphere[NX*NY], int npoints, int n
    free(ftheta);
 }
 void compute_cubic_grid_sphere_weights(t_wave_sphere wsphere[NX*NY], int npoints)
 /* compute weights of corrected Laplacian on cubic grid */
 {
    int i, k, p;
    double dist, wmax, ratio;
    printf("Computing grid Laplacian weights\n"); 
    wmax = 0.0;
    for (i=0; i<npoints; i++) 
    {
        dist = 0.0;
        for (k=0; k<4; k++)
        {
            p = wsphere[i].neighbor[k];
            dist += dist_sphere(wsphere[i].phigrid, wsphere[i].thetagrid, wsphere[p].phigrid, wsphere[p].thetagrid);
        }
        dist *= 0.25;
        wsphere[i].laplace_weight = 1.0/(dist*dist);
        if (wsphere[i].laplace_weight > wmax) wmax = wsphere[i].laplace_weight;
    }
    ratio = 1.0/wmax;
    printf("Ratio = %.3lg\n", ratio); 
    for (i=0; i<npoints; i++) 
    {
        wsphere[i].laplace_weight*= ratio;
 //         if (wsphere[i].laplace_weight > 1.5) wsphere[i].laplace_weight = 1.5;
 //         if (wsphere[i].laplace_weight < 0.5) wsphere[i].laplace_weight = 0.5;
        printf("Grid point %i weight = %.3lg\n", i, wsphere[i].laplace_weight); 
    }
 }
 int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY])
 /* initialize Poisson random simulation grid on a sphere */
 {
@@ -9526,7 +9978,8 @@ int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY])
        case (GRID_CUBIC):
        {
            npoints = generate_cubic_grid_sphere(wsphere);
-            if (OPTIMIZE_GRID) optimize_cubic_grid_sphere(wsphere, npoints, 600);
+            if (OPTIMIZE_GRID) optimize_cubic_grid_sphere(wsphere, npoints, OPTIMIZE_STEPS);
            if (WEIGHT_GRID) compute_cubic_grid_sphere_weights(wsphere, npoints);
            break;
        }
    }
@@ -9541,7 +9994,7 @@ int initialize_simulation_grid_sphere(t_wave_sphere wsphere[NX*NY])
    return(npoints);
 }
-void convert_fields_from_grids(double phi_in[NX*NY], double phi_out[NX*NY], t_wave_sphere wsphere[NX*NY])
+void convert_fields_from_grids(double phi_in[NX*NY], double phi_out[NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
 /* convert field phi_in in grid coordinates to field phi_out in latitude-longitude coordinates */
 {
 //     int i;
@@ -353,21 +353,21 @@ double gaussian()
-// int in_polygon(double x, double y, double r, int npoly, double apoly)
+int in_polygon(double x, double y, double r, int npoly, double apoly)
-// /* test whether (x,y) is in regular polygon of npoly sides inscribed in circle of radious r, turned by apoly Pi/2 */
+/* test whether (x,y) is in regular polygon of npoly sides inscribed in circle of radious r, turned by apoly Pi/2 */
-// {
+{
-//     int condition = 1, k;
+    int condition = 1, k;
-//     double omega, cw, angle; 
+    double omega, cw, angle; 
-//     
+    
-//     omega = DPI/((double)npoly);
+    omega = DPI/((double)npoly);
-//     cw = cos(omega*0.5);
+    cw = cos(omega*0.5);
-//     for (k=0; k<npoly; k++)  
+    for (k=0; k<npoly; k++)  
-//     {
+    {
-//         angle = apoly*PID + ((double)k+0.5)*omega;
+        angle = apoly*PID + ((double)k+0.5)*omega;
-//         condition = condition*(x*cos(angle) + y*sin(angle) < r*cw);
+        condition = condition*(x*cos(angle) + y*sin(angle) < r*cw);
-//     }
+    }
-//     return(condition);
+    return(condition);
-// }
+}
@@ -3621,6 +3621,85 @@ int compute_circle_pairs_coordinates(t_vertex polyline[NMAXPOLY], t_arc polyarc[
    return(nlines);
 }
 int compute_circle_poly_coordinates(t_vertex polyline[NMAXPOLY], t_arc polyarc[NMAXPOLY], t_circle circle[NMAXCIRCLES], int *npolyarc, int *ncircles, short int top)
 /* initialise lines for D_CIRCLE_POLYGON domain */
 {
    double dx, alpha, alphab, h1, h2, h1b, h2b, x1, y1, pos[2], mu, beta;
    int i, j, k, nlines, narcs, ncirc;
    nlines = 0;
    narcs = 0;
    ncirc = 0;
    dx = (XMAX - XMIN)/(double)NGRIDX;
    h1 = sqrt(MU*MU - 0.25*WALL_WIDTH*WALL_WIDTH);
    h2 = LAMBDA - h1;
    alpha = atan(0.5*WALL_WIDTH/h1);
    h1b = sqrt(MU_B*MU_B - 0.25*WALL_WIDTH*WALL_WIDTH);
 //     h2b = LAMBDA - h1b;
    h2b = MU_B*cos(PI/(double)NPOLY) - LAMBDA;
    alphab = atan(0.5*WALL_WIDTH/h1b);
    for (i=0; i<NGRIDX; i++)
    {
        x1 = XMIN + ((double)i+0.5)*dx;
        y1 = LAMBDA;
        polyarc[narcs].xc = x1;
        polyarc[narcs].yc = y1;
        polyarc[narcs].r  = MU;
        polyarc[narcs].angle1 = -PID + alpha;
        polyarc[narcs].dangle = DPI - 2.0*alpha;
        narcs++;
        polyline[nlines].x = x1 + 0.5*WALL_WIDTH;
        polyline[nlines].y = h2;
        nlines++;
        polyline[nlines].x = x1 + 0.5*WALL_WIDTH;
        polyline[nlines].y = h2b;
        nlines++;
        polyline[nlines].x = x1 + 0.5*WALL_WIDTH;
        polyline[nlines].y = h2b;
        nlines++;
        for (j=0; j<NPOLY; j++)
        {
            beta = PID + PI*(double)(-1-2*j)/(double)NPOLY;
            polyline[nlines].x = x1 + MU_B*cos(beta);
            polyline[nlines].y = -y1 + MU_B*sin(beta);
            nlines++;
            polyline[nlines].x = x1 + MU_B*cos(beta);
            polyline[nlines].y = -y1 + MU_B*sin(beta);
            nlines++;
        }
        polyline[nlines].x = x1 - 0.5*WALL_WIDTH;
        polyline[nlines].y = h2b;
        nlines++;
        polyline[nlines].x = x1 - 0.5*WALL_WIDTH;
        polyline[nlines].y = h2b;
        nlines++;
        polyline[nlines].x = x1 - 0.5*WALL_WIDTH;
        polyline[nlines].y = h2;
        nlines++;
    }
    for (i=0; i<nlines; i++)
    {
        xy_to_pos(polyline[i].x, polyline[i].y, pos);        
        polyline[i].posi = pos[0];
        polyline[i].posj = pos[1];
    }
    *npolyarc = narcs;
    *ncircles = ncirc;
    return(nlines);
 }
 int compute_disc_honeycomb_lattice_coordinates(t_vertex polyline[NMAXPOLY], t_rect_rotated polyrect[NMAXPOLY], t_arc polyarc[NMAXPOLY], int *npolyrect_rot, int *npolyarc)
 /* initialise lines for D_CIRCLE_LATTICE_HEX domain */
 {
@@ -5060,6 +5139,10 @@ int init_poly(int depth, t_vertex polyline[NMAXPOLY], t_rectangle polyrect[NMAXP
        {
            return(compute_circle_pairs_coordinates(polyline, polyarc, circles, npolyarc, ncircles, top));
        }
        case (D_CIRCLE_POLYGON):
        {
            return(compute_circle_poly_coordinates(polyline, polyarc, circles, npolyarc, ncircles, top));
        }
        case (D_TWO_PARABOLAS_ASYM_GUIDE):
        {
            polyrect[0].x1 = XMIN - WALL_WIDTH;
@@ -7208,6 +7291,19 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
            }
            return(0);
        }
        case (D_CIRCLE_POLYGON):
        {
            dx = (XMAX - XMIN)/(double)NGRIDX;
            for (k=0; k<NGRIDX; k++)
            {
                x1 = XMIN + dx*((double)k + 0.5);
                y1 = LAMBDA;
                if (module2(x-x1, y-y1) < MU) return(1);
                if (in_polygon(x-x1, y+y1, MU_B, NPOLY, 1.0 - 2.0/(double)NPOLY)) return(1);
                if ((vabs(y) < LAMBDA)&&(vabs(x-x1) < 0.5*WALL_WIDTH)) return(1);
            }
            return(0);
        }
        case (D_MENGER):       
        {
            x1 = 0.5*(x+1.0);
@@ -10269,6 +10365,21 @@ void draw_billiard(int fade, double fade_value)      /* draws the billiard bound
            }
            break;
        }
        case (D_CIRCLE_POLYGON):
        {
            glBegin(GL_LINES);
            for (i=0; i<npolyline; i++) 
            {
                glVertex2d(polyline[i].posi, polyline[i].posj);
            }
            glEnd();
            for (i=0; i<npolyarc; i++)
            {
                draw_circle_arc(polyarc[i].xc, polyarc[i].yc, polyarc[i].r, polyarc[i].angle1, polyarc[i].dangle, NSEG);
            }
            break;
        }
        case (D_MENGER):
        {
            glLineWidth(3);
@@ -98,6 +98,303 @@ void init_wave_sphere_rde(t_wave_sphere *wsphere, int res)
    }
 }
 double dist_sphere(double phi1, double theta1, double phi2, double theta2)
 /* returns angle between points given in spherical coordinates */
 {
    double dist;
    dist = sin(theta1)*sin(theta2)*cos(phi1 - phi2);
    dist += cos(theta1)*cos(theta2);
    return(acos(dist));
 }
 int generate_poisson_sample_sphere(double dpoisson, int nmaxpoints, double *px, double *py)
 /* generate a Poisson disc sample on the sphere */
 {
    int i, j, k, n_p_active, ncandidates=500, naccepted, p, q, npoints; 
    double r, phi, x, y, z, x1, y1, z1, x2, y2, distance, dist1, polar_padding;
    short int active_poisson[NX*NY], far;
    printf("Generating Poisson disc sample\n");
    /* generate first point */
    px[0] = DPI*(double)rand()/RAND_MAX;
    py[0] = polar_padding + (PI - 2.0*polar_padding)*(double)rand()/RAND_MAX;
    active_poisson[0] = 1;
    n_p_active = 1;
    npoints = 1;
    polar_padding = dpoisson;
    while ((n_p_active > 0)&&(npoints < nmaxpoints))
    {
        /* randomly select an active circle */
        i = rand()%(npoints);
        while (!active_poisson[i]) i = rand()%(npoints);                 
        /* generate new candidates */
        naccepted = 0;
        for (j=0; j<ncandidates; j++)
        {
            r = dpoisson*(2.0*(double)rand()/RAND_MAX + 1.0);
            phi = DPI*(double)rand()/RAND_MAX;
            /* random point centered at north pole */
            x = sin(r)*cos(phi);
            y = sin(r)*sin(phi);
            z = cos(r);
            /* rotation by theta = py[i] */
            x1 = x*cos(py[i]) + z*sin(py[i]);
            y1 = y;
            z1 = -x*sin(py[i]) + z*cos(py[i]);
            /* rotation by phi = px[i] */
            x2 = x1*cos(px[i]) - y1*sin(px[i]);
            y2 = x1*sin(px[i]) + y1*cos(px[i]);
            y = acos(z1);
            x = argument(x2, y2);
            if (x < 0.0) x += DPI;
 //             printf("current point (%.3lg, %.3lg), new point (%.3lg, %.3lg), distance = %.3lg\n", px[i], py[i], x, y, dist_sphere(x, y, px[i], py[i])); 
 //             x = px[i] + r*cos(phi);
 //             y = py[i] + r*sin(phi);
 //             if (x > DPI) x -= DPI;
 //             if (x < 0.0) x += DPI;
 //             
 //             far = 1;
 //             if (y > PI - polar_padding) far = 0;
 //             if (y < polar_padding) far = 0;
            far = 1;
            dist1 = 10.0;
            for (k=0; k<npoints; k++) /*if ((k!=i))*/
            {
                /* new point is far away from point k */
                distance = dist_sphere(x, y, px[k], py[k]);
                far = far*(distance >= dpoisson);
                if (distance < dist1) dist1 = distance;
            }
            if ((far)&&(npoints < nmaxpoints))    /* accept new point */
            {
                printf("New grid point at (%.3f,%.3f) accepted\n", x, y);
                printf("Minimal distance = %.3lg\n", dist1); 
                printf("%i grid points generated\n", npoints);
                px[npoints] = x;
                py[npoints] = y;
                active_poisson[npoints] = 1;
                npoints++;
                n_p_active++;
                naccepted++;
            }
        }
        if (naccepted == 0)    /* inactivate circle i */ 
        {
            active_poisson[i] = 0;
            n_p_active--;
        }
        printf("%i active points\n", n_p_active);
    }
    printf("Generated %i points\n", npoints);
    return(npoints);
 }
 int init_circle_sphere(t_circles_sphere *circles, int circle_pattern)
 /* initialise circles on sphere */
 /* for billiard shape D_SPHERE_CIRCLES */
 {
    int ncircles, k;
    double alpha, beta, gamma;
    double px[NMAXCIRC_SPHERE], py[NMAXCIRC_SPHERE];
    switch (circle_pattern) {
        case (C_SPH_DODECA):
        {
            alpha = acos(sqrt(5.0)/3.0);
            beta = acos(1.0/3.0);
            gamma = asin(sqrt(3.0/8.0));
            circles[0].theta = 0.0;
            circles[0].phi = 0.0;
            for (k=0; k<3; k++)
            {
                circles[1+k].theta = alpha;
                circles[1+k].phi = (double)k*DPI/3.0;
            }
            for (k=0; k<3; k++)
            {
                circles[4+k].theta = beta;
                circles[4+k].phi = (double)k*DPI/3.0 + gamma;
                circles[7+k].theta = beta;
                circles[7+k].phi = (double)k*DPI/3.0 - gamma;
            }
            for (k=0; k<3; k++)
            {
                circles[10+k].theta = PI - beta;
                circles[10+k].phi = (double)k*DPI/3.0 + PI/3.0 + gamma;
                circles[13+k].theta = PI - beta;
                circles[13+k].phi = (double)k*DPI/3.0 + PI/3.0 - gamma;
            }
            for (k=0; k<3; k++)
            {
                circles[16+k].theta = PI - alpha;
                circles[16+k].phi = (double)k*DPI/3.0 + PI/3.0;
            }
            circles[19].theta = PI;
            circles[19].phi = 0.0;
            for (k=0; k<20; k++)
            {
                circles[k].radius = MU;
                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
                circles[k].z = cos(circles[k].theta);
            }
            ncircles = 20;
            break;
        }
        case (C_SPH_ICOSA):
        {
            alpha = acos(1.0/sqrt(5.0));
            circles[0].theta = 0.0;
            circles[0].phi = 0.0;
            for (k=0; k<5; k++)
            {
                circles[1+k].theta = alpha;
                circles[1+k].phi = (double)k*DPI/5.0;
            }
            for (k=0; k<5; k++)
            {
                circles[6+k].theta = PI - alpha;
                circles[6+k].phi = (double)k*DPI/5.0 + PI/5.0;
            }
            circles[11].theta = PI;
            circles[11].phi = 0.0;
            for (k=0; k<12; k++) 
            {
                circles[k].radius = MU;
                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
                circles[k].z = cos(circles[k].theta);
            }
            ncircles = 12;
            break;
        }
        case (C_SPH_OCTA):
        {
            circles[0].theta = 0.0;
            circles[0].phi = 0.0;
            for (k=0; k<4; k++)
            {
                circles[1+k].theta = PID;
                circles[1+k].phi = (double)k*PID;
            }
            circles[5].theta = PI;
            circles[5].phi = 0.0;
            for (k=0; k<6; k++) 
            {
                circles[k].radius = MU;
                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
                circles[k].z = cos(circles[k].theta);
            }
            ncircles = 6;
            break;
        }
        case (C_SPH_CUBE):
        {
            alpha = acos(1.0/3.0);
            circles[0].theta = 0.0;
            circles[0].phi = 0.0;
            for (k=0; k<3; k++)
            {
                circles[1+k].theta = alpha;
                circles[1+k].phi = (double)k*DPI/3.0;
                circles[4+k].theta = PI - alpha;
                circles[4+k].phi = (double)k*DPI/3.0 + PI/3.0;
            }
            circles[7].theta = PI;
            circles[7].phi = 0.0;
            for (k=0; k<8; k++) 
            {
                circles[k].radius = MU;
                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
                circles[k].z = cos(circles[k].theta);
            }
            ncircles = 8;
            break;
        }
        case (C_SPH_CUBE_B):
        {
            alpha = acos(1.0/sqrt(3.0));
            for (k=0; k<4; k++)
            {
                circles[k].theta = alpha;
                circles[k].phi = ((double)k + 0.5)*PID;
                circles[4+k].theta = PI - alpha;
                circles[4+k].phi = ((double)k + 0.5)*PID;
            }
            for (k=0; k<8; k++) 
            {
                circles[k].radius = MU;
                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
                circles[k].z = cos(circles[k].theta);
            }
            ncircles = 8;
            break;
        }
        case (C_SPH_POISSON):
        {
            ncircles = generate_poisson_sample_sphere(PDISC_DISTANCE, NMAXCIRC_SPHERE, px, py);
            for (k=0; k<ncircles; k++) 
            {
                circles[k].phi = px[k];
                circles[k].theta = py[k];
                circles[k].radius = MU;
                circles[k].x = sin(circles[k].theta)*cos(circles[k].phi);
                circles[k].y = sin(circles[k].theta)*sin(circles[k].phi);
                circles[k].z = cos(circles[k].theta);
            }
            break;
        }
        default: 
        {
            printf("Function init_circle_sphere not defined for this pattern \n");
        }
    }
    return(ncircles);
 }
 void read_negative_dem_values_rde(double *height_values, t_wave_sphere *wsphere)
 /* init bathymetric data */
@@ -44,7 +44,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 */
@@ -72,7 +72,7 @@
 /* Choice of the billiard table */
-#define B_DOMAIN 421         /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 991         /* choice of domain shape, see list in global_pdes.c */
 #define CIRCLE_PATTERN 202  /* pattern of circles or polygons, see list in global_pdes.c */
 #define IMAGE_FILE 5        /* for option D_IMAGE */
@@ -87,17 +87,17 @@
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 #define PDISC_CONNECT_FACTOR 1.5    /* controls which discs are connected for D_CIRCLE_LATTICE_POISSON domain */
-#define LAMBDA 0.9	    /* parameter controlling the dimensions of domain */
+#define LAMBDA 0.6	    /* parameter controlling the dimensions of domain */
-#define MU -0.2124612       /* parameter controlling the dimensions of domain */
+#define MU 0.15             /* parameter controlling the dimensions of domain */
-#define MU_B 0.3            /* parameter controlling the dimensions of domain */
+#define MU_B 0.45           /* parameter controlling the dimensions of domain */
-#define NPOLY 10             /* number of sides of polygon */
+#define NPOLY 5             /* number of sides of polygon */
-#define APOLY 0.0          /* angle by which to turn polygon, in units of Pi/2 */ 
+#define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define MDEPTH 6            /* depth of computation of Menger gasket */
 #define MRATIO 3            /* ratio defining Menger gasket */
 #define MANDELLEVEL 1000    /* iteration level for Mandelbrot set */
 #define MANDELLIMIT 10.0    /* limit value for approximation of Mandelbrot set */
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
-#define NGRIDX 6            /* number of grid point for grid of disks */
+#define NGRIDX 4            /* number of grid point for grid of disks */
 #define NGRIDY 14           /* number of grid point for grid of disks */
 #define WALL_WIDTH 0.075    /* width of wall separating lenses */
 #define WALL_WIDTH_B 0.1    /* width of wall separating lenses */
@@ -122,7 +122,7 @@
 /* Physical parameters of wave equation */
 #define TWOSPEEDS 0         /* set to 1 to replace hardcore boundary by medium with different speed */
-#define OSCILLATE_LEFT 1    /* set to 1 to add oscilating boundary condition on the left */
+#define OSCILLATE_LEFT 0    /* set to 1 to add oscilating boundary condition on the left */
 #define OSCILLATE_TOPBOT 0   /* set to 1 to enforce a planar wave on top and bottom boundary */
 #define OSCILLATION_SCHEDULE 42  /* oscillation schedule, see list in global_pdes.c */
 #define OSCIL_YMAX 0.0375          /* defines oscilling beam range */
@@ -134,7 +134,7 @@
 #define AMPLITUDE 2.0      /* amplitude of periodic excitation */ 
 #define ACHIRP 0.25        /* acceleration coefficient in chirp */
 #define DAMPING 0.0        /* damping of periodic excitation */
-#define COURANT 0.08       /* Courant number in medium B */
+#define COURANT 0.12        /* Courant number in medium B */
 #define COURANTB 0.25      /* Courant number */
 #define GAMMA 0.0          /* damping factor in wave equation */
 #define GAMMAB 0.0         /* damping factor in wave equation */
@@ -149,10 +149,10 @@
 /* Increasing COURANT speeds up the simulation, but decreases accuracy */
 /* For similar wave forms, COURANT^2*GAMMA should be kept constant */
-#define ADD_OSCILLATING_SOURCE 0        /* set to 1 to add an oscillating wave source */
+#define ADD_OSCILLATING_SOURCE 1        /* set to 1 to add an oscillating wave source */
-#define OSCILLATING_SOURCE_PERIOD 6.0   /* period of oscillating source */
+#define OSCILLATING_SOURCE_PERIOD 12.0   /* period of oscillating source */
 #define ALTERNATE_OSCILLATING_SOURCE 1  /* set to 1 to alternate sign of oscillating source */
-#define N_SOURCES 1                     /* number of sources, for option draw_sources */
+#define N_SOURCES 4                     /* number of sources, for option draw_sources */
 #define ALTERNATE_SOURCE_PHASES 0       /* set to 1 to alternate initial phases of sources */
 #define MAX_PULSING_TIME 5000           /* max time for adding pulses */
@@ -169,10 +169,10 @@
 /* Parameters for length and speed of simulation */
-#define NSTEPS 2700          /* number of frames of movie */
+#define NSTEPS 3200          /* number of frames of movie */
 #define NVID 12              /* number of iterations between images displayed on screen */
 #define NSEG 1000            /* number of segments of boundary */
-#define INITIAL_TIME 100      /* time after which to start saving frames */
+#define INITIAL_TIME 0      /* time after which to start saving frames */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
 #define PRINT_SPEED 0       /* print speed of moving source */
 #define PRINT_FREQUENCY 0       /* print frequency (for phased array) */
@@ -199,7 +199,7 @@
 /* Color schemes */
-#define COLOR_PALETTE 18     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 15     /* Color palette, see list in global_pdes.c  */
 #define COLOR_PALETTE_B 18   /* Color palette, see list in global_pdes.c  */
 #define BLACK 1          /* background */
@@ -212,8 +212,8 @@
 #define PHASE_FACTOR 1.0       /* factor in computation of phase in color scheme P_3D_PHASE */
 #define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
 #define ATTENUATION 0.0   /* exponential attenuation coefficient of contrast with time */
-#define VSHIFT_AMPLITUDE -0.5    /* additional shift for wave amplitude */
+#define VSHIFT_AMPLITUDE -0.0    /* additional shift for wave amplitude */
-#define VSCALE_AMPLITUDE 1.0     /* additional scaling factor for wave amplitude */
+#define VSCALE_AMPLITUDE 0.5     /* additional scaling factor for wave amplitude */
 #define E_SCALE 300.0       /* scaling factor for energy representation */
 #define LOG_SCALE 0.75     /* scaling factor for energy log representation */
 #define LOG_SHIFT 0.75     /* shift of colors on log scale */
@@ -231,14 +231,14 @@
 #define HUEMEAN 180.0    /* mean value of hue for color scheme C_HUE */
 #define HUEAMP -180.0      /* amplitude of variation of hue for color scheme C_HUE */
-#define DRAW_COLOR_SCHEME 1     /* set to 1 to plot the color scheme */
+#define DRAW_COLOR_SCHEME 0     /* set to 1 to plot the color scheme */
 #define COLORBAR_RANGE 1.7      /* scale of color scheme bar */
 #define COLORBAR_RANGE_B 4.5    /* scale of color scheme bar for 2nd part */
 #define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
 #define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
 #define CIRC_COLORBAR_B 0       /* set to 1 to draw circular color scheme */
-#define DRAW_WAVE_PROFILE 1     /* set to 1 to draw a profile of the wave */
+#define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 #define HORIZONTAL_WAVE_PROFILE 1 /* set to 1 to draw wave profile vertically */
 #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
 #define WAVE_PROFILE_X 1.9      /* value of x to sample wave profile */
@@ -880,11 +880,14 @@ void animation()
        /* add oscillating waves */
        for (source = 0; source < N_SOURCES; source++)
        {
-            wave_source_x[source] = polyarc[2*source].xc;
+            wave_source_x[source] = polyarc[source].xc;
-            wave_source_y[source] = polyarc[2*source].yc;
+            wave_source_y[source] = polyarc[source].yc;
            source_periods[source] = OSCILLATING_SOURCE_PERIOD/(1.0 + 0.5*(double)source);
            source_amp[source] = INITIAL_AMP;
        }
        source_periods[0] = OSCILLATING_SOURCE_PERIOD;
        source_periods[1] = OSCILLATING_SOURCE_PERIOD*3/4;
        source_periods[2] = OSCILLATING_SOURCE_PERIOD*2/3;
        source_periods[3] = OSCILLATING_SOURCE_PERIOD/2;
        for (source = 0; source < N_SOURCES; source++)
        {
            dperiod = source_periods[source];