Add files via upload

global_3d.c

@@ -38,6 +38,8 @@
 #define GF_WING 4           /* wing shape */
 #define GF_COMPUTE_FROM_BC 5    /* compute force field as gradient of bc_field2 */
 #define GF_EARTH 6          /* field depends on altitude on continents */
+#define GF_MARS 7           /* field depends on altitude on Mars */
+#define GF_VENUS 8          /* field depends on altitude on Venus */
 
 /* Choice of water depth for shallow water equation */
 
@@ -45,10 +47,14 @@
 #define SH_CIRCLES 2    /* shallow obstacle specified by CIRCLE_PATTERN */
 #define SH_COAST 3      /* depth varying with x-coordinate */
 #define SH_COAST_MONOTONE 4      /* depth decreasing with x-coordinate */
+#define SH_SPHERE_CUBE 5    /* cube embedded in sphere */
+#define SH_SPHERE_OCTAHEDRON 6   /* octahedron embedded in sphere */
+#define SH_SPHERE_DODECAHEDRON 7 /* dodecahedron embedded in sphere */
+#define SH_SPHERE_ICOSAHEDRON 8  /* icosahedron embedded in sphere */
 
 /* Type of rotating viewpoint */
 
-#define VP_HORIZONTAL 0     /* rotate in a horizontal plane */
+#define VP_HORIZONTAL 0     /* rotate in a horizontal plane (constant latitude) */
 #define VP_ORBIT 1          /* rotate in a plane containing the origin */
 #define VP_ORBIT2 11        /* rotate in a plane specified by max latitude */
 #define VP_POLAR 2          /* polar orbit */
@@ -83,7 +89,7 @@
 #define PLANET ((B_DOMAIN == D_SPHERE_EARTH)||(B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)||(B_DOMAIN == D_SPHERE_VENUS)||(B_DOMAIN == D_SPHERE_MERCURY))
 #define OTHER_PLANET ((B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON)||(B_DOMAIN == D_SPHERE_VENUS)||(B_DOMAIN == D_SPHERE_MERCURY))
 
-#define RDE_PLANET ((ADAPT_STATE_TO_BC)&&(OBSTACLE_GEOMETRY == D_SPHERE_EARTH))
+#define RDE_PLANET ((ADAPT_STATE_TO_BC)&&((OBSTACLE_GEOMETRY == D_SPHERE_EARTH)||(OBSTACLE_GEOMETRY == D_SPHERE_MARS)||(OBSTACLE_GEOMETRY == D_SPHERE_VENUS)))
 
 #define NMAXCIRC_SPHERE 100     /* max number of circles on sphere */
 #define NMAX_TRACER_PTS 20       /* max number of tracer points recorded per cell */

global_ljones.c

@@ -11,14 +11,16 @@
 #define D_CIRCLES 20    /* several circles */
 #define D_CIRCLES_IN_RECT 201   /* several circles in a rectangle */
 
-#define NMAXCIRCLES 100000       /* total number of circles/polygons (must be at least NCX*NCY for square grid) */
+#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 1000  /* max number of obstacles */
 #define NMAXSEGMENTS 1000   /* max number of repelling segments */
 #define NMAXGROUPS 50       /* max number of groups of segments */
 #define NMAXCOLLISIONS 200000   /* max number of collisions */
 #define NMAXPARTNERS 30     /* max number of partners in molecule */
-#define NMAXPARTNERMOLECULES 10 /* max number of partners of a molecule */
+#define NMAXPARTNERMOLECULES 30 /* max number of partners of a molecule */
+#define NMAXPARTINCLUSTER 500   /* max number of particles in cluster */
+#define NMAXBELTS 10        /* max number of conveyor belts */
 
 #define C_SQUARE 0          /* square grid of circles */
 #define C_HEX 1             /* hexagonal/triangular grid of circles */
@@ -76,6 +78,7 @@
 #define S_DAM_BRICKS 17     /* dam made of several bricks */
 #define S_HLINE_HOLE 18    /* horizontal line with a hole in the bottom */
 #define S_HLINE_HOLE_SPOKES 181    /* horizontal line with a hole in the bottom and extra spokes */
+#define S_HLINE_HOLE_SLOPED 182 /* slanted lines with a hole */
 #define S_EXT_CIRCLE_RECT 19    /* particles outside a circle and a rectangle */
 #define S_BIN_OPENING 20        /* bin containing particles opening at deactivation time */
 #define S_BIN_LARGE 201         /* larger bin */
@@ -88,6 +91,10 @@
 #define S_CYLINDER 27           /* walls at top and bottom, for cylindrical b.c. */
 #define S_TREE 28               /* Christmas tree(s) */
 #define S_CONE 29               /* cone */
+#define S_CONVEYOR_BELT 30      /* conveyor belt */
+#define S_TWO_CONVEYOR_BELTS 31 /* two angled conveyor belts */
+#define S_PERIODIC_CONVEYORS 32 /* one wrapping belt, and one short horizontal belt */
+#define S_TEST_CONVEYORS 321    /* test */
 
 /* particle interaction */
 
@@ -108,6 +115,10 @@
 #define I_COULOMB_IMAGINARY 14  /* Coulomb interaction with "imaginary charge" */
 #define I_DNA_CHARGED 15    /* Coulomb-type interaction between end points of DNA nucleotides */
 #define I_DNA_CHARGED_B 151    /* stronger Coulomb-type interaction between end points of DNA nucleotides */
+#define I_SEGMENT 16           /* harmonic interaction between segments */
+#define I_SEGMENT_CHARGED 161  /* harmonic interaction between segments and Coulomb interaction between ends*/
+#define I_POLYGON 17           /* harmonic interaction between regular polygons */
+#define I_POLYGON_ALIGN 171    /* harmonic interaction between polygons with an aligning torque */
 
 /* Boundary conditions */
 
@@ -208,6 +219,9 @@
 #define CHEM_DNA_BASE_SPLIT 254 /* aggregation/splitting of DNA molecules when base pairs don't match */
 #define CHEM_DNA_ENZYME 255 /* aggregation/splitting of DNA molecules in presence of enzymes */
 #define CHEM_DNA_ENZYME_REPAIR 256 /* aggregation/splitting of DNA molecules in presence of enzymes and additional repairing of bad connections */
+#define CHEM_POLYGON_AGGREGATION 26 /* aggregation of polygons */
+#define CHEM_POLYGON_CLUSTER 261    /* clustering of polygons into new clusters */
+#define CHEM_POLYGON_ONECLUSTER 262 /* clustering of polygons, with only one cluster allowed */
 
 /* Initial conditions for chemical reactions */
 
@@ -254,12 +268,16 @@
 #define P_NUMBER 12       /* colors depend on particle number */
 #define P_EMEAN 13        /* averaged kinetic energy (with exponential damping) */
 #define P_LOG_EMEAN 131   /* log of averaged kinetic energy (with exponential damping) */
+#define P_EMEAN_DENSITY 132 /* averaged kinetic energy divided by the cluster size */
 #define P_DIRECT_EMEAN 14 /* averaged version of P_DIRECT_ENERGY */
 #define P_NOPARTICLE 15   /* particles are not drawn (only the links between them) */
 #define P_NPARTNERS 16    /* number of partners */
 #define P_CHARGE 17       /* colors represent charge */
 #define P_MOL_ANGLE 18    /* orientation of molecule defined by partners */
 #define P_CLUSTER 19      /* colors depend on connected component */
+#define P_CLUSTER_SIZE 20 /* colors depend on size of connected component */
+#define P_CLUSTER_SELECTED 21   /* colors show which clusters are slected for growth */
+#define P_COLLISION 22    /* colors depend on number of collision/reaction */
 
 /* Rotation schedules */
 
@@ -275,6 +293,8 @@
 
 #define POLY_STAR 0     /* star-shaped graph (central molecule attracts outer ones) */
 #define POLY_ALL 1      /* all-to-all coupling */
+#define POLY_STAR_CHARGED 11     /* star-shaped graph with charged molecules */
+#define POLY_POLYGON 12 /* polygonal shape */
 #define POLY_WATER 2    /* star-shaped with a 120° separation between anions */
 #define POLY_SOAP 3     /* polymers with all-to-all coupling and polar end */
 #define POLY_SOAP_B 4   /* polymers with pairwise coupling and polar end */
@@ -287,6 +307,9 @@
 #define POLY_DNA_ALT 71     /* simplified model for DNA with different short ends */
 #define POLY_DNA_DOUBLE 72  /* simplified model for DNA with double ends for rigidity */
 #define POLY_DNA_FLEX 73    /* simplified model for DNA with less backbone rigidity (beta) */
+#define POLY_KITE 8         /* kite for kites and darts quasicrystal */
+#define POLY_DART 81        /* dart for kites and darts quasicrystal */
+#define POLY_SEG_POLYGON 9  /* polygon of segments */
 
 /* Background color schemes */
 
@@ -365,18 +388,42 @@ typedef struct
     int partner[NMAXPARTNERS];  /* partner particles for option PAIR_PARTICLES */
     short int npartners;        /* number of partner particles */
     double partner_eqd[NMAXPARTNERS];   /* equilibrium distances between partners */
+    double partner_eqa[NMAXPARTNERS];   /* equilibrium angle between partners */
     int p0, p1;                 /* numbers of two first partners (for P_MOL_ANGLE color scheme) */
 //     short int mol_angle;        /* for color scheme P_MOL_ANGLE */
     int cluster;                /* number of cluster */
+    int cluster_color;          /* color of cluster */
+    int cluster_size;           /* size of cluster */
     int molecule;               /* number of molecule */
     short int tested, cactive;  /* for cluster search */
     short int coulomb;          /* has value 1 if DNA-Coulomb interaction is attractive */
     short int added;            /* has value 1 if particle has been added */
     short int reactive;         /* has value 1 if particle can react */
     short int paired;           /* has value 1 if belongs to base-paired molecule */
+    short int flip;             /* keeps track of which particles in a cluster are flipped by PI */
     int partner_molecule;       /* number of partner molecule */
+    int collision;              /* number of collision */
 } t_particle;
 
+typedef struct
+{
+    short int active;           /* has value 1 if cluster is active */
+    short int thermostat;       /* has value 1 if cluster is coupled to thermostat */
+    short int selected;         /* has value 1 if cluster is selected to be able to grow */
+    double xg, yg;              /* center of gravity */
+    double vx, vy;              /* velocity of center of gravity */
+    double angle;               /* orientation of cluster */
+    double omega;               /* angular velocity of cluster */
+    double mass, mass_inv;      /* mass of cluster and its inverse */
+    double inertia_moment, inertia_moment_inv;   /* moment of inertia */
+    double fx, fy, torque;      /* force and torque */
+    double energy, emean;       /* energy and averaged energy */
+    double dirmean;             /* time-averaged direction */
+    int nparticles;             /* number of particles in cluster */
+    int particle[NMAXPARTINCLUSTER];    /* list of particles in cluster */
+//     int angle_ref;              /* reference particle for orientation */
+} t_cluster;
+
 typedef struct
 {
     int number;                 /* total number of particles in cell */
@@ -433,6 +480,8 @@ typedef struct
     double pressure;            /* pressure acting on segement */
     double avrg_pressure;       /* time-averaged pressure */
     short int inactivate;       /* set to 1 for segment to become inactive at time SEGMENT_DEACTIVATION_TIME */
+    short int conveyor;         /* set to 1 for segment to exert lateral force */
+    double conveyor_speed;      /* speed of conveyor belt */
 } t_segment;
 
 typedef struct
@@ -468,13 +517,24 @@ typedef struct
 typedef struct
 {
     int nparticles;             /* number of particles */
-    int particle[NPARTNERS+1];  /* list of particles */
+    int particle[2*NPARTNERS+1];  /* list of particles */
     int npartners;              /* number of partner molecules */
     int partner[NMAXPARTNERMOLECULES];  /* list of partner molecules */
     int connection_type[NMAXPARTNERMOLECULES];  /* types of particles in connection */
     short int added;            /* has value 1 if molecule has been added */
 } t_molecule;
 
+typedef struct 
+{
+    double x1, y1, x2, y2;      /* positions of extremities */
+    double width;               /* width of belt */
+    double speed;               /* speed of conveyor belt */
+    double position;            /* position of belt (needed for display of rotating parts) */
+    double length;              /* distance between (x1,x2) and (y1,y2) */
+    double angle;               /* angle of (x1,x2) - (y1,y2) */
+    double tx, ty;              /* coordinates of tangent vector */
+} t_belt;
+
 typedef struct
 {
     int nactive;                /* number of active particles */
@@ -496,6 +556,6 @@ typedef struct
 
 
 
-int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0;
+int frame_time = 0, ncircles, nobstacles, nsegments, ngroups = 1, counter = 0, nmolecules = 0, nbelts = 0;
 FILE *lj_log;
 

global_pdes.c

@@ -103,6 +103,9 @@
 #define D_RITCHEY_CHRETIEN_SPHERICAL 75   /* Ritchey-Chrétien telescope with spherical mirrors */
 #define D_RITCHEY_CHRETIEN_HYPERBOLIC 751 /* Ritchey-Chrétien telescope with hyperbolic mirrors */
 #define D_GRADIENT_INDEX_LENS 76    /* gradient index lens (only affects draw_billiard) */ 
+#define D_IMAGE 77              /* Taken from image file */
+#define D_MAGNETRON 78          /* simplified magnetron */
+#define D_MAGNETRON_CATHODE 781 /* simplified magnetron with central cathode */
 
 /* for wave_sphere.c */
 
@@ -141,7 +144,9 @@
 
 #define C_RINGS 20          /* obstacles arranged in concentric rings */
 #define C_RINGS_T 201       /* obstacles arranged in concentric rings, triangular lattice */
-#define C_RINGS_SPIRAL 202  /* obstacles arranged on a "subflower" spiral, similar to C_GOLDEN_SPIRAL */
+#define C_RINGS_SPIRAL 202  /* obstacles arranged on a "sunflower" spiral, similar to C_GOLDEN_SPIRAL */
+#define C_RINGS_POISSONDISC 203 /* obstacles arranged in a Poisson disc pattern */
+#define C_RINGS_LOGSPIRAL 204   /* logarithmic spirals */
 
 #define C_HEX_BOTTOM 101    /* hex/triangular lattice in lower half */
 #define C_HEX_BOTTOM2 102   /* smaller hex/triangular lattice in lower half */

heat.c

@@ -70,6 +70,7 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 300        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.25   /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 1.1	    /* parameter controlling the dimensions of domain */
@@ -158,6 +159,7 @@
 // #define SLOPE 0.1        /* sensitivity of color on wave amplitude */
 #define SLOPE 0.2        /* sensitivity of color on wave amplitude */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
+#define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
 
 #define COLORHUE 260     /* initial hue of water color for scheme C_LUM */
 #define COLORDRIFT 0.0   /* how much the color hue drifts during the whole simulation */
@@ -219,6 +221,7 @@
 #define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
 #define WALL_WIDTH 0.1      /* width of wall separating lenses */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 #define INITIAL_TIME 50      /* time after which to start saving frames */
 #define OSCIL_YMAX 0.35      /* defines oscillation range */
 #define MESSAGE_LDASH 14         /* length of dash for Morse code message */

lennardjones.c

@@ -36,8 +36,8 @@
 #include <omp.h>
 #include <time.h>
 
-#define MOVIE 1         /* set to 1 to generate movie */
+#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 while saving frames */
 #define NO_EXTRA_BUFFER_SWAP 1    /* some OS require one less buffer swap when recording images */
 
@@ -58,10 +58,10 @@
 #define YMIN -1.125
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
-#define INITXMIN -1.95
+#define INITXMIN -1.8
-#define INITXMAX 1.95	/* x interval for initial condition */
+#define INITXMAX 1.8	/* x interval for initial condition */
-#define INITYMIN -1.1
+#define INITYMIN 0.9
-#define INITYMAX 1.1	/* y interval for initial condition */
+#define INITYMAX 2.2	/* y interval for initial condition */
 
 #define THERMOXMIN -1.25
 #define THERMOXMAX 1.25	/* x interval for initial condition */
@@ -78,7 +78,7 @@
 #define BCXMIN -2.0
 #define BCXMAX 2.0	/* x interval for boundary condition */
 #define BCYMIN -1.125
-#define BCYMAX 1.125	/* y interval for boundary condition */
+#define BCYMAX 2.4	/* y interval for boundary condition */
 
 #define OBSXMIN -2.0
 #define OBSXMAX 2.0     /* x interval for motion of obstacle */
@@ -91,14 +91,16 @@
 #define ADD_FIXED_OBSTACLES 0   /* set to 1 do add fixed circular obstacles */
 #define OBSTACLE_PATTERN 71  /* pattern of obstacles, see list in global_ljones.c */
 
-#define ADD_FIXED_SEGMENTS 0    /* set to 1 to add fixed segments as obstacles */
+#define ADD_FIXED_SEGMENTS 1    /* set to 1 to add fixed segments as obstacles */
-#define SEGMENT_PATTERN 29      /* pattern of repelling segments, see list in global_ljones.c */
+#define SEGMENT_PATTERN 182      /* pattern of repelling segments, see list in global_ljones.c */
 #define ROCKET_SHAPE 3        /* shape of rocket combustion chamber, see list in global_ljones.c */
 #define ROCKET_SHAPE_B 3      /* shape of second rocket */
 #define NOZZLE_SHAPE 6        /* shape of nozzle, see list in global_ljones.c */
 #define NOZZLE_SHAPE_B 6      /* shape of nozzle for second rocket, see list in global_ljones.c */
+#define BELT_SPEED1 20.0    /* speed of first conveyor belt */
+#define BELT_SPEED2 -10.0   /* speed of second conveyor belt */
 
-#define TWO_TYPES 1         /* set to 1 to have two types of particles */
+#define TWO_TYPES 0         /* set to 1 to have two types of particles */
 #define TYPE_PROPORTION 0.5 /* proportion of particles of first type */
 #define TWOTYPE_CONFIG 0    /* choice of types, see TTC_ list in global_ljones.c */
 #define SYMMETRIZE_FORCE 1  /* set to 1 to symmetrize two-particle interaction, only needed if particles are not all the same */
@@ -106,31 +108,31 @@
 #define CENTER_PY 0         /* set to 1 to center vertical momentum */
 #define CENTER_PANGLE 0     /* set to 1 to center angular momentum */
 
-#define INTERACTION 12       /* particle interaction, see list in global_ljones.c */
+#define INTERACTION 171     /* particle interaction, see list in global_ljones.c */
-#define INTERACTION_B 12     /* 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 2.0 /* angular frequency of spin-spin interaction */
+#define SPIN_INTER_FREQUENCY 6.0 /* angular frequency of spin-spin interaction */
-#define SPIN_INTER_FREQUENCY_B 2.0 /* angular frequency of spin-spin interaction for second particle type */
+#define SPIN_INTER_FREQUENCY_B 6.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 MOL_ANGLE_FACTOR 6.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 9.6  /* minimal distance in Poisson disc process, controls density of particles */
+#define PDISC_DISTANCE 2.5  /* 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.0087 	    /* parameter controlling radius of particles */
+#define MU 0.025 	    /* parameter controlling radius of particles */
-#define MU_B 0.012          /* parameter controlling radius of particles of second type */
+#define MU_B 0.03           /* parameter controlling radius of particles of second type */
-#define NPOLY 40            /* number of sides of polygon */
+#define NPOLY 6             /* number of sides of polygon */
-#define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
+#define APOLY 0.075           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define AWEDGE 0.5          /* opening angle of wedge, in units of Pi/2 */ 
 #define MDEPTH 4            /* depth of computation of Menger gasket */
 #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 40           /* number of grid point for grid of disks */
+#define NGRIDX 36           /* number of grid point for grid of disks */
-#define NGRIDY 20           /* number of grid point for grid of disks */
+#define NGRIDY 36           /* number of grid point for grid of disks */
 #define EHRENFEST_RADIUS 0.9    /* radius of container for Ehrenfest urn configuration */
 #define EHRENFEST_WIDTH 0.035     /* width of tube for Ehrenfest urn configuration */
 #define TWO_CIRCLES_RADIUS_RATIO 0.8    /* ratio of radii for S_TWO_CIRCLES_EXT segment configuration */
@@ -146,11 +148,10 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 2500      /* number of frames of movie */
+#define NSTEPS 2700      /* number of frames of movie */
-// #define NSTEPS 200      /* number of frames of movie */
+#define NVID 200          /* number of iterations between images displayed on screen */
-#define NVID 50          /* number of iterations between images displayed on screen */
 #define NSEG 25          /* number of segments of boundary of circles */
-#define INITIAL_TIME 5     /* time after which to start saving frames */
+#define INITIAL_TIME 200     /* time after which to start saving frames */
 #define OBSTACLE_INITIAL_TIME 0     /* time after which to start moving obstacle */
 #define BOUNDARY_WIDTH 1    /* width of particle boundary */
 #define LINK_WIDTH 2        /* width of links between particles */
@@ -161,45 +162,49 @@
 #define SLEEP1  1        /* initial sleeping time */
 #define SLEEP2  1   /* final sleeping time */
 #define MID_FRAMES 100    /* number of still frames between parts of two-part movie */
-// #define END_FRAMES 250   /* number of still frames at end of movie */
 #define END_FRAMES 100   /* number of still frames at end of movie */
 
 /* Boundary conditions, see list in global_ljones.c */
 
-#define BOUNDARY_COND 3
+#define BOUNDARY_COND 1
 
 /* Plot type, see list in global_ljones.c  */
 
-#define PLOT 5
+#define PLOT 11
-#define PLOT_B 13        /* plot type for second movie */
+#define PLOT_B 14        /* plot type for second movie */
 
 /* Background color depending on particle properties */
 
-#define COLOR_BACKGROUND 1  /* set to 1 to color background */
+#define COLOR_BACKGROUND 0  /* set to 1 to color background */
 #define BG_COLOR 2          /* type of background coloring, see list in global_ljones.c */
 #define BG_COLOR_B 0        /* type of background coloring, see list in global_ljones.c */
 
-#define DRAW_BONDS 1    /* set to 1 to draw bonds between neighbours */
+#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 */
 #define FILL_TRIANGLES 0    /* set to 1 to fill triangles between neighbours */
+#define DRAW_CLUSTER_LINKS 0    /* set to 1 to draw links between particles in cluster */
 #define ALTITUDE_LINES 0    /* set to 1 to add horizontal lines to show altitude */
 #define COLOR_SEG_GROUPS 0  /* set to 1 to collor segment groups differently */
 #define N_PARTICLE_COLORS 200   /* number of colors for P_NUMBER color scheme */
 #define INITIAL_POS_TYPE 0     /* type of initial position dependence */
 #define ERATIO 0.995          /* ratio for time-averaging in P_EMEAN color scheme */
-#define DRATIO 0.995          /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */
+#define DRATIO 0.999          /* ratio for time-averaging in P_DIRECT_EMEAN color scheme */
 
 /* Color schemes */
 
 #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 10       /* Color palette for angle representation */
+#define COLOR_PALETTE_ANGLE 0        /* Color palette for angle representation */
-#define COLOR_PALETTE_DIRECTION 0    /* Color palette for direction representation */
+#define COLOR_PALETTE_DIRECTION 17   /* 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 */
 #define COLOR_PALETTE_PRESSURE 11    /* Color palette for different neighbours representation */
 #define COLOR_PALETTE_CHARGE 18      /* Color palette for charge representation */
-#define COLOR_PALETTE_CLUSTER 11     /* Color palette for cluster representation */
+#define COLOR_PALETTE_CLUSTER 14     /* Color palette for cluster representation */
+#define COLOR_PALETTE_CLUSTER_SIZE 13 /* Color palette for cluster size representation */
+#define COLOR_PALETTE_CLUSTER_SELECTED 11 /* Color palette for selected cluster representation */
+#define COLOR_HUE_CLUSTER_SELECTED 90.0    /* Color hue for selected cluster */
+#define COLOR_HUE_CLUSTER_NOT_SELECTED 220.0    /* Color hue for selected cluster */
 
 #define BLACK 1          /* background */
 
@@ -234,12 +239,11 @@
 #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_EMAX 50000.0        /* energy of particle with hottest color */
 #define PARTICLE_EMIN 10.0          /* energy of particle with coolest color */
-#define PARTICLE_EMAX 20000.0        /* energy of particle with hottest color */
+#define PARTICLE_EMAX 100.0        /* energy of particle with hottest color */
-#define HUE_TYPE0 300.0      /* hue of particles of type 0 */
+#define HUE_TYPE0 320.0      /* hue of particles of type 0 */
-#define HUE_TYPE1 00.0       /* hue of particles of type 1 */
+#define HUE_TYPE1 60.0       /* hue of particles of type 1 */
-#define HUE_TYPE2 340.0      /* hue of particles of type 2 */
+#define HUE_TYPE2 320.0      /* hue of particles of type 2 */
 #define HUE_TYPE3 260.0      /* hue of particles of type 3 */
 #define HUE_TYPE4 200.0      /* hue of particles of type 4 */
 #define HUE_TYPE5 60.0       /* hue of particles of type 5 */
@@ -249,39 +253,40 @@
 
 #define RANDOM_RADIUS 0     /* set to 1 for random circle radius */
 #define DT_PARTICLE 3.0e-6    /* time step for particle displacement */
-#define KREPEL 15.0          /* constant in repelling force between particles */
+#define KREPEL 50.0          /* constant in repelling force between particles */
-#define EQUILIBRIUM_DIST 5.0    /* Lennard-Jones equilibrium distance */
+#define EQUILIBRIUM_DIST 2.3    /* Lennard-Jones equilibrium distance */
-#define EQUILIBRIUM_DIST_B 5.0  /* Lennard-Jones equilibrium distance for second type of particle */
+#define EQUILIBRIUM_DIST_B 2.5  /* Lennard-Jones equilibrium distance for second type of particle */
 #define REPEL_RADIUS 25.0    /* radius in which repelling force acts (in units of particle radius) */
 #define DAMPING 500.0          /* damping coefficient of particles */
-#define INITIAL_DAMPING 5000.0  /* damping coefficient of particles during initial phase */
+#define INITIAL_DAMPING 1000.0  /* damping coefficient of particles during initial phase */
-#define DAMPING_ROT 1.0e6      /* damping coefficient for rotation of particles */
+#define DAMPING_ROT 1000.0      /* damping coefficient for rotation of particles */
 #define PARTICLE_MASS 2.0    /* mass of particle of radius MU */
-#define PARTICLE_MASS_B 16.0   /* mass of particle of radius MU_B */
+#define PARTICLE_MASS_B 2.0   /* mass of particle of radius MU_B */
-#define PARTICLE_INERTIA_MOMENT 0.2     /* moment of inertia of particle */
+#define PARTICLE_INERTIA_MOMENT 0.5     /* moment of inertia of particle */
-#define PARTICLE_INERTIA_MOMENT_B 0.02     /* moment of inertia of second type of particle */
+#define PARTICLE_INERTIA_MOMENT_B 0.5     /* moment of inertia of second type of particle */
 #define V_INITIAL 50.0        /* initial velocity range */
-#define OMEGA_INITIAL 10.0        /* initial angular velocity range */
+#define OMEGA_INITIAL 50.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 V_INITIAL_TYPE 0    /* type of initial speed distribution (see VI_ in global_ljones.c) */
 
-#define THERMOSTAT 1        /* set to 1 to switch on thermostat */
+#define THERMOSTAT 0        /* set to 1 to switch on thermostat */
 #define VARY_THERMOSTAT 0   /* set to 1 for time-dependent thermostat schedule */
 #define SIGMA 5.0           /* noise intensity in thermostat */
-#define BETA 0.00007          /* initial inverse temperature */
+#define BETA 0.004          /* 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 2.0e11    /* harmonic potential of obstacles */
+#define KSPRING_OBSTACLE 1.0e9    /* harmonic potential of obstacles */
-#define NBH_DIST_FACTOR 5.0        /* radius in which to count neighbours */
+#define NBH_DIST_FACTOR 4.0        /* radius in which to count neighbours */
-#define GRAVITY 0.0            /* gravity acting on all particles */
+#define GRAVITY 4000.0            /* gravity acting on all particles */
 #define GRAVITY_X 0.0          /* horizontal gravity acting on all particles */
 #define CIRCULAR_GRAVITY 0     /* set to 1 to have gravity directed to center */
 #define INCREASE_GRAVITY 0     /* set to 1 to increase gravity during the simulation */
 #define GRAVITY_SCHEDULE 1     /* type of gravity schedule, see list in global_ljones.c */
 #define GRAVITY_FACTOR 10.0     /* factor by which to increase gravity */
-#define GRAVITY_INITIAL_TIME 200    /* time at start of simulation with constant gravity */
+#define GRAVITY_INITIAL_TIME 250    /* time at start of simulation with constant gravity */
 #define GRAVITY_RESTORE_TIME 500    /* time at end of simulation with gravity restored to initial value */
 #define KSPRING_VICSEK 0.2   /* spring constant for I_VICSEK_SPEED interaction */
 #define VICSEK_REPULSION 10.0    /* repulsion between particles in Vicsek model */
@@ -291,9 +296,8 @@
 #define ADD_BFIELD 0     /* set to 1 to add a magnetic field */
 #define BFIELD 2.666666667       /* value of magnetic field */
 #define CHARGE 1.0      /* charge of particles of first type */
-#define CHARGE_B -1.5     /* charge of particles of second type */
+#define CHARGE_B 1.0     /* charge of particles of second type */
 #define INCREASE_E 0     /* set to 1 to increase electric field */
-// #define EFIELD_FACTOR 2500000.0    /* factor by which to increase electric field */
 #define EFIELD_FACTOR 5000000.0    /* factor by which to increase electric field */
 #define INCREASE_B 0     /* set to 1 to increase magnetic field */
 #define BFIELD_FACTOR 20000.0    /* factor by which to increase magnetic field */
@@ -301,25 +305,24 @@
 #define OBSTACLE_CHARGE 3.0     /* charge of obstacles */
 #define KCOULOMB_OBSTACLE 1000.0   /* Coulomb force constant for charged obstacles */
 
-#define ROTATION 0           /* set to 1 to include rotation of particles */
+#define ROTATION 1           /* set to 1 to include rotation of particles */
 #define COUPLE_ANGLE_TO_THERMOSTAT 0    /* set to 1 to couple angular degrees of freedom to thermostat */
 #define DIMENSION_FACTOR 1.0  /* scaling factor taking into account number of degrees of freedom */  
-#define KTORQUE -100.0         /* force constant in angular dynamics */
+#define KTORQUE 5.0e3         /* force constant in angular dynamics */
-#define KTORQUE_BOUNDARY 1.0e6  /* constant in torque from the boundary */
+#define KTORQUE_BOUNDARY 5.0e5  /* constant in torque from the boundary */
 #define KTORQUE_B 10.0        /* force constant in angular dynamics */
-#define KTORQUE_DIFF -150.0    /* force constant in angular dynamics for different particles */
+#define KTORQUE_DIFF 500.0    /* force constant in angular dynamics for different particles */
 #define DRAW_SPIN 0           /* set to 1 to draw spin vectors of particles */
 #define DRAW_SPIN_B 0         /* set to 1 to draw spin vectors of particles */
-#define DRAW_CROSS 1          /* set to 1 to draw cross on particles of second type */
+#define DRAW_CROSS 0          /* set to 1 to draw cross on particles of second type */
-#define DRAW_MINUS 1          /* set to 1 to draw cross on particles of negative charge */
+#define DRAW_MINUS 0          /* set to 1 to draw cross on particles of negative charge */
 #define SPIN_RANGE 10.0       /* range of spin-spin interaction */
 #define SPIN_RANGE_B 10.0     /* range of spin-spin interaction for second type of particle */
 #define QUADRUPOLE_RATIO 0.6  /* anisotropy in quadrupole potential */ 
 
-#define INCREASE_BETA 1  /* set to 1 to increase BETA during simulation */
+#define INCREASE_BETA 0  /* set to 1 to increase BETA during simulation */
 #define BETA_SCHEDULE 3    /* type of temperature schedule, see TS_* in global_ljones */
-#define BETA_FACTOR 10.0    /* factor by which to change BETA during simulation */
+#define BETA_FACTOR 0.06    /* factor by which to change BETA during simulation */
-// #define BETA_FACTOR 0.08    /* 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 */
@@ -380,7 +383,7 @@
 #define SMOOTH_ROTATION 1           /* set to 1 to update segments at each time step (rather than at each movie frame) */
 #define ROTATION_SCHEDULE 0         /* time-dependence of rotation angle, see ROT_* in global_ljones.c */
 #define PERIOD_ROTATE_BOUNDARY 1000  /* period of rotating boundary */
-#define ROTATE_INITIAL_TIME 300       /* initial time without rotation */
+#define ROTATE_INITIAL_TIME 150       /* initial time without rotation */
 #define ROTATE_FINAL_TIME 300       /* final time without rotation */
 #define ROTATE_CHANGE_TIME 0.5     /* relative duration of acceleration/deceleration phases */
 #define OMEGAMAX -2.0*PI              /* maximal rotation speed */
@@ -388,7 +391,7 @@
 #define MOVE_BOUNDARY 0        /* set to 1 to move repelling segments, due to force from particles */
 #define SEGMENTS_MASS 40.0     /* mass of collection of segments */
 #define DEACTIVATE_SEGMENT 1    /* set to 1 to deactivate last segment after a certain time */
-#define SEGMENT_DEACTIVATION_TIME 20   /* time at which to deactivate last segment */
+#define SEGMENT_DEACTIVATION_TIME 200   /* time at which to deactivate last segment */
 #define RELEASE_ROCKET_AT_DEACTIVATION 0    /* set to 1 to limit segments velocity before segment release */
 #define SEGMENTS_X0 1.5        /* initial position of segments */
 #define SEGMENTS_Y0 0.0        /* initial position of segments */
@@ -405,6 +408,7 @@
 #define SEGMENT_GROUP_DAMPING 0.0   /* damping of segment groups */
 #define GROUP_REPULSION 0           /* set to 1 for groups of segments to repel each other */
 #define KSPRING_GROUPS 5.0e11       /* harmonic potential between segment groups */
+#define KSPRING_BELT 1.0e4          /* spring constant from belt */
 #define GROUP_WIDTH 0.05            /* interaction width of groups */
 #define GROUP_G_REPEL 0             /* set to 1 to add repulsion between centers of mass of groups */
 #define GROUP_G_REPEL_RADIUS 1.2    /* radius within which centers of mass of groups repel each other */
@@ -418,20 +422,29 @@
 #define PRINT_ENTROPY 0     /* set to 1 to compute entropy */
 
 #define REACTION_DIFFUSION 0      /* set to 1 to simulate a chemical reaction (particles may change type) */
-#define RD_REACTION 256           /* type of reaction, see list in global_ljones.c */
+#define RD_REACTION 262           /* type of reaction, see list in global_ljones.c */
-#define RD_TYPES 6                /* number of types in reaction-diffusion equation */
+#define RD_TYPES 1                /* number of types in reaction-diffusion equation */
-#define RD_INITIAL_COND 10        /* initial condition of particles */
+#define RD_INITIAL_COND 0         /* initial condition of particles */
-#define REACTION_DIST 4.0         /* maximal distance for reaction to occur */
+#define REACTION_DIST 1.2         /* maximal distance for reaction to occur */
 #define REACTION_PROB 1.0         /* probability controlling reaction term */ 
 #define DISSOCIATION_PROB 0.0     /* probability controlling dissociation reaction */ 
-#define KILLING_PROB 0.0015        /* probability of enzymes being killed */
+#define KILLING_PROB 0.0015       /* probability of enzymes being killed */
+#define DELTAMAX 0.1              /* max orientation difference for pairing polygons */
 #define CENTER_COLLIDED_PARTICLES 0  /* set to 1 to recenter particles upon reaction (may interfere with thermostat) */
 #define EXOTHERMIC 0            /* set to 1 to make reaction exo/endothermic */
 #define DELTA_EKIN 2000.0       /* change of kinetic energy in reaction */
-#define COLLISION_TIME 25      /* time during which collisions are shown */
+#define COLLISION_TIME 25       /* time during which collisions are shown */
-#define DELTAVMAX 1000.0         /* maximal deltav allowed for pairing molecules */
+#define COLLISION_RADIUS 2.0    /* radius of discs showing collisions, in units of MU */
-#define AGREGMAX 11              /* maximal number of partners for CHEM_AGGREGATION reaction */
+#define DELTAVMAX 200.0         /* maximal deltav allowed for pairing molecules */
-#define AGREG_DECOUPLE 12        /* minimal number of partners to decouple from thermostat */
+#define AGREGMAX 6              /* maximal number of partners for CHEM_AGGREGATION reaction */
+#define AGREG_DECOUPLE 12       /* minimal number of partners to decouple from thermostat */
+#define CLUSTER_PARTICLES 0     /* set to 1 for particles to form rigid clusters */
+#define CLUSTER_MAXSIZE 1000     /* max size of clusters */
+#define SMALL_CLUSTER_MAXSIZE 10 /* size limitation on smaller cluster */
+#define SMALL_NP_MAXSIZE 2      /* limitation on number of partners of particle in smaller cluster */
+#define NOTSELECTED_CLUSTER_MAXSIZE 0   /* limit on size of clusters that can merge with non-selected cluster */
+#define REPAIR_CLUSTERS 0       /* set to 1 to repair alignment in clusters */
+#define REPAIR_MIN_DIST 0.75    /* relative distance below which overlapping polygons are inactivated */
 
 #define CHANGE_RADIUS 0         /* set to 1 to change particle radius during simulation */
 #define MU_RATIO 0.666666667    /* ratio by which to increase radius */
@@ -458,35 +471,36 @@
 #define PROP_MIN 0.1        /* min proportion of type 1 particles */
 #define PROP_MAX 0.9        /* max proportion of type 1 particles */
 
-#define PAIR_PARTICLES 1    /* set to 1 to form particles pairs */
+#define PAIR_PARTICLES 0    /* set to 1 to form particles pairs */
 #define RANDOMIZE_ANGLE 0   /* set to 1 for random orientation */
 #define DEACIVATE_CLOSE_PAIRS 1 /* set to 1 to test for closeness to other particles */
 #define PAIR_SAFETY_FACTOR 1.2  /* distance to deactivate divided by sum of radii */
 #define THIRD_TYPE_PROPORTION 1.0   /* proportion of third type pairings, for certain pairing types */
 
-#define KSPRING_PAIRS 1.0e10    /* spring constant for pair interaction */
+#define KSPRING_PAIRS 5.0e10    /* spring constant for pair interaction */
 #define KTORQUE_PAIRS 1.0e10   /* constant for angular coupling in pair interaction */
 #define KTORQUE_PAIR_ANGLE 0.0    /* constant for coupling between orientation in pairs */
-#define NPARTNERS 8        /* number of partners of particles */
+#define NPARTNERS 4         /* number of partners of particles - for DNA, set NPARTNERS_DNA */
-#define NARMS 1              /* number of "arms" for certain paring types */ 
+#define NPARTNERS_DNA 8     /* number of partners of particles, case of DNA, should be at least 8 */
-#define PAIRING_TYPE 42      /* type of pairing, see POLY_ in global_ljones.c */
+#define NARMS 5             /* number of "arms" for certain paring types */ 
+#define PAIRING_TYPE 9      /* type of pairing, see POLY_ in global_ljones.c */
 #define PARTNER_ANGLE 104.45    /* angle (in degrees) between ions for POLY_WATER case */
-#define PAIR_DRATIO 0.9      /* ratio between equilibrium distance and radius (default: 1.0) */
+#define PAIR_DRATIO 1.1      /* ratio between equilibrium distance and radius (default: 1.0) */
-#define MU_C 0.0125            /* radius of partner particle */
+#define MU_C 0.035            /* radius of partner particle */
-#define PARTICLE_MASS_C 8.0  /* mass or partner particle */
+#define PARTICLE_MASS_C 2.0  /* mass or partner particle */
-#define CHARGE_C -1.0         /* charge of partner particle */
+#define CHARGE_C 1.0         /* charge of partner particle */
-#define CLUSTER_COLOR_FACTOR 400  /* factor for initialization of cluster colors */
+#define CLUSTER_COLOR_FACTOR 40  /* factor for initialization of cluster colors */
 #define ALTERNATE_POLY_CHARGE 1   /* set to 1 for alternating charges in molecule */
 #define SECONDARY_PAIRING 0     /* set to 1 to pair with secondary partners, experimental */
 #define DNA_RIGIDITY 0.5     /* controls rigidity for POLY_DNA_DOUBLE pairs, default = 1 */
 
 #define PAIR_TYPEB_PARTICLES 1  /* set to 1 to pair particle of type 1 */
-#define NPARTNERS_B 2         /* number of partners of particles */
+#define NPARTNERS_B 5         /* number of partners of particles */
 #define NARMS_B 1               /* number of "arms" for certain paring types */ 
-#define PAIRING_TYPE_B 2      /* type of pairing, see POLY_ in global_ljones.c */
+#define PAIRING_TYPE_B 81     /* type of pairing, see POLY_ in global_ljones.c */
-#define MU_D 0.008            /* radius of partner particle */
+#define MU_D 0.035            /* radius of partner particle */
-#define PARTICLE_MASS_D 1.0  /* mass or partner particle */
+#define PARTICLE_MASS_D 2.0  /* mass or partner particle */
-#define CHARGE_D 0.75         /* charge of partner particle */
+#define CHARGE_D 1.0         /* charge of partner particle */
 
 #define NXMAZE 12      /* width of maze */
 #define NYMAZE 12      /* height of maze */
@@ -500,8 +514,8 @@
 #define FLOOR_OMEGA 0      /* set to 1 to limit particle momentum to PMAX */
 #define PMAX 1000.0        /* maximal force */
 
-#define HASHX 80     /* size of hashgrid in x direction */
+#define HASHX 100     /* size of hashgrid in x direction */
-#define HASHY 40     /* size of hashgrid in y direction */
+#define HASHY 88      /* 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 */
 
@@ -515,11 +529,12 @@
 #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))
 #define TWO_OBSTACLES ((SEGMENT_PATTERN == S_TWO_CIRCLES_EXT)||(SEGMENT_PATTERN == S_TWO_ROCKETS))
-#define COMPUTE_EMEAN ((PLOT == P_EMEAN)||(PLOT_B == P_EMEAN)||(PLOT == P_LOG_EMEAN)||(PLOT_B == P_LOG_EMEAN)||(PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN))
+#define COMPUTE_EMEAN ((PLOT == P_EMEAN)||(PLOT_B == P_EMEAN)||(PLOT == P_LOG_EMEAN)||(PLOT_B == P_LOG_EMEAN)||(PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN)||(PLOT == P_EMEAN_DENSITY)||(PLOT_B == P_EMEAN_DENSITY))
 #define COMPUTE_DIRMEAN ((PLOT == P_DIRECT_EMEAN)||(PLOT_B == P_DIRECT_EMEAN))
 #define COUNT_PARTNER_TYPE ((RD_REACTION == CHEM_H2O_H_OH)||(RD_REACTION == CHEM_2H2O_H3O_OH))
-#define PAIR_FORCE ((PAIR_PARTICLES)||((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_AGGREGATION)||(RD_REACTION == CHEM_AGGREGATION_CHARGE)||(RD_REACTION == CHEM_AGGREGATION_NNEIGH))))
+#define PAIR_FORCE ((PAIR_PARTICLES)||((REACTION_DIFFUSION)&&((RD_REACTION == CHEM_AGGREGATION)||(RD_REACTION == CHEM_AGGREGATION_CHARGE)||(RD_REACTION == CHEM_AGGREGATION_NNEIGH)||(RD_REACTION == CHEM_POLYGON_AGGREGATION))))
 #define COMPUTE_PAIR_TORQUE (KTORQUE_PAIR_ANGLE != 0.0)
+#define ADD_CONVEYOR_FORCE ((ADD_FIXED_SEGMENTS)&&((SEGMENT_PATTERN == S_CONVEYOR_BELT)||(SEGMENT_PATTERN == S_TWO_CONVEYOR_BELTS)||(SEGMENT_PATTERN == S_PERIODIC_CONVEYORS)||(SEGMENT_PATTERN == S_TEST_CONVEYORS)))
 
 double xshift = 0.0;      /* x shift of shown window */
 double xspeed = 0.0;      /* x speed of obstacle */
@@ -960,6 +975,8 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
     if (initial_phase) damping = INITIAL_DAMPING;
     else damping = DAMPING;
     
+//     printf("Evolving particles\n");
+    
     #pragma omp parallel for private(j,xi,totalenergy,a,move)
     for (j=0; j<ncircles; j++) if (particle[j].active)
     {
@@ -979,7 +996,9 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
         if (COMPUTE_DIRMEAN)
         {
             direction = argument(particle[j].vx, particle[j].vy);
+//             printf("direction = %.3lg\t", direction);
             dmean = particle[j].dirmean;
+//             printf("dirmean = %.3lg\n", particle[j].dirmean);
             if (dmean < direction - PI) dmean += DPI;
             else if (dmean > direction + PI) dmean -= DPI;
             particle[j].dirmean = DRATIO*dmean + (1.0-DRATIO)*direction;
@@ -1092,6 +1111,135 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
     return(totalenergy);
 }
 
+double evolve_clusters(t_particle particle[NMAXCIRCLES], t_cluster cluster[NMAXCIRCLES], 
+                       t_hashgrid hashgrid[HASHX*HASHY], 
+                        double cqx[NMAXCIRCLES], double cqy[NMAXCIRCLES], double cqangle[NMAXCIRCLES],
+                        double cpx[NMAXCIRCLES], double cpy[NMAXCIRCLES], double cpangle[NMAXCIRCLES], 
+                        double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase, int verbose)
+{
+    double a, totalenergy = 0.0, damping, direction, dmean, newx, newy, newangle, deltax, deltay, deltaangle;  
+    static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0;
+    int j, move, ncoup;
+    
+    if (initial_phase) damping = INITIAL_DAMPING;
+    else damping = DAMPING;
+    
+    #pragma omp parallel for private(j,xi,totalenergy,a,move)
+    for (j=0; j<ncircles; j++) if (cluster[j].active)
+    {                
+        cluster[j].vx = cpx[j] + 0.5*DT_PARTICLE*cluster[j].fx;
+        cluster[j].vy = cpy[j] + 0.5*DT_PARTICLE*cluster[j].fy;
+        cluster[j].omega = cpangle[j] + 0.5*DT_PARTICLE*cluster[j].torque;
+        
+        cpx[j] = cluster[j].vx + 0.5*DT_PARTICLE*cluster[j].fx;
+        cpy[j] = cluster[j].vy + 0.5*DT_PARTICLE*cluster[j].fy;
+        cpangle[j] = cluster[j].omega + 0.5*DT_PARTICLE*cluster[j].torque;
+        
+        cluster[j].energy = (cpx[j]*cpx[j] + cpy[j]*cpy[j])*cluster[j].mass_inv;
+                
+        if (COMPUTE_EMEAN)
+            cluster[j].emean = ERATIO*cluster[j].emean + (1.0-ERATIO)*cluster[j].energy;
+        
+        if (COMPUTE_DIRMEAN)
+        {
+            direction = argument(cluster[j].vx, cluster[j].vy);
+            dmean = cluster[j].dirmean;
+            if (dmean < direction - PI) dmean += DPI;
+            else if (dmean > direction + PI) dmean -= DPI;
+            cluster[j].dirmean = DRATIO*dmean + (1.0-DRATIO)*direction;
+            if (cluster[j].dirmean < 0.0) cluster[j].dirmean += DPI;
+            else if (cluster[j].dirmean > DPI) cluster[j].dirmean -= DPI;
+        }
+        
+        if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(cluster[j].thermostat))
+            cluster[j].energy += cpangle[j]*cpangle[j]*cluster[j].inertia_moment_inv;
+        
+        cqx[j] = cluster[j].xg + 0.5*DT_PARTICLE*cpx[j]*cluster[j].mass_inv;
+        cqy[j] = cluster[j].yg + 0.5*DT_PARTICLE*cpy[j]*cluster[j].mass_inv;
+        cqangle[j] = cluster[j].angle + 0.5*DT_PARTICLE*cpangle[j]*cluster[j].inertia_moment_inv;
+                
+        if ((THERMOSTAT_ON)&&(cluster[j].thermostat))
+        {
+            cpx[j] *= exp(- 0.5*DT_PARTICLE*xi);
+            cpy[j] *= exp(- 0.5*DT_PARTICLE*xi);
+        }
+        if ((COUPLE_ANGLE_TO_THERMOSTAT)&&(cluster[j].thermostat)) 
+            cpangle[j] *= exp(- 0.5*DT_PARTICLE*xi);
+    }
+                
+    /* compute kinetic energy */
+//     *nactive = 0;
+    ncoup = 1;
+    for (j=0; j<ncircles; j++)
+        if ((cluster[j].active)&&(cluster[j].thermostat))
+        {
+            totalenergy += cluster[j].energy;
+            ncoup++;
+//             *nactive++;
+        }
+        
+    totalenergy *= DIMENSION_FACTOR;    /* normalize energy to take number of degrees of freedom into account */
+    if (THERMOSTAT_ON)
+    {
+        /* TODO - fix nactive vs ncoupled */
+//         a = DT_PARTICLE*(totalenergy - (double)*nactive/beta)/MU_XI;
+        a = DT_PARTICLE*(totalenergy - (double)ncoup/beta)/MU_XI;
+        a += SIGMA*sqrt(DT_PARTICLE)*gaussian();
+        xi = (xi + a - b*xi)/(1.0 + b);
+    }
+    
+    move = 0;
+    for (j=0; j<ncircles; j++) if (cluster[j].active) 
+    {
+        if ((THERMOSTAT_ON)&&(cluster[j].thermostat))
+        {
+            cpx[j] *= exp(- 0.5*DT_PARTICLE*xi);
+            cpy[j] *= exp(- 0.5*DT_PARTICLE*xi);
+            if (!COUPLE_ANGLE_TO_THERMOSTAT) cpangle[j] *= exp(- DT_PARTICLE*DAMPING_ROT);
+        }
+        else 
+        {
+            cpx[j] *= exp(- DT_PARTICLE*damping);
+            cpy[j] *= exp(- DT_PARTICLE*damping);
+            cpangle[j] *= exp(- DT_PARTICLE*DAMPING_ROT);
+//             printf("Damping cluster angular velocity\n");
+        }
+        if ((THERMOSTAT_ON)&&(COUPLE_ANGLE_TO_THERMOSTAT)&&(cluster[j].thermostat))
+            cpangle[j] *= exp(- 0.5*DT_PARTICLE*xi);
+        
+        newx = cqx[j] + 0.5*DT_PARTICLE*cpx[j]*cluster[j].mass_inv;
+        newy = cqy[j] + 0.5*DT_PARTICLE*cpy[j]*cluster[j].mass_inv;
+        newangle = cqangle[j] + 0.5*DT_PARTICLE*cpangle[j]*cluster[j].inertia_moment_inv;
+        
+        deltax = newx - cluster[j].xg;
+        deltay = newy - cluster[j].yg;
+        deltaangle = newangle - cluster[j].angle;
+                
+//         translate_cluster(j, cluster, particle, deltax, deltay, 0);
+//         rotate_cluster(j, cluster, particle, deltaangle);
+        translate_and_rotate_cluster(j, cluster, particle, deltax, deltay, deltaangle);
+                
+//         cluster[j].vx = cpx[j] + 0.5*DT_PARTICLE*cluster[j].fx;
+//         cluster[j].vy = cpy[j] + 0.5*DT_PARTICLE*cluster[j].fy;
+//         cluster[j].omega = cpangle[j] + 0.5*DT_PARTICLE*cluster[j].torque;
+
+        /* FIXME: adapt to clusters */
+//         else if (!NO_WRAP_BC)
+//         {
+//             move += wrap_particle(&particle[j], &px[j], &py[j]);
+//         }
+//                 if (move > 0) 
+//                 {
+//                     compute_relative_positions(particle, hashgrid);
+//                     update_hashgrid(particle, hashgrid, 0);   /* REDUNDANT ? */
+//                 }
+    }
+    
+//     sleep(1);
+    
+    *ncoupled = ncoup;
+    return(totalenergy);
+}
 
 void evolve_lid(double fboundary)
 {
@@ -1394,16 +1542,18 @@ 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;
+            angle, theta, sum, alpha;
     double *qx, *qy, *px, *py, *qangle, *pangle, *pressure, *obstacle_speeds;
-    int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, 
+    double *cqx, *cqy, *cpx, *cpy, *cqangle, *cpangle;
+    int i, j, k, n, m, s, ij[2], i0, iplus, iminus, j0, jplus, jminus, p, q, p1, q1, p2, q2, total_neighbours = 0, cl, 
         min_nb, max_nb, close, wrapx = 0, wrapy = 0, nadd_particle = 0, nmove = 0, nsuccess = 0, 
         tracer_n[N_TRACER_PARTICLES], traj_position = 0, traj_length = 0, move = 0, old, m0, floor, nthermo, wall = 0,
-        group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1;
+        group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1, np, belt;
     int *particle_numbers;
     static int imin, imax;
     static short int first = 1;
     t_particle *particle;
+    t_cluster *cluster;
     t_obstacle *obstacle;
     t_segment *segment;
     t_group_segments *segment_group;     
@@ -1413,6 +1563,7 @@ void animation()
     t_hashgrid *hashgrid;
     t_molecule *molecule;
     t_lj_parameters params;
+    t_belt *conveyor_belt;
     char message[100];
 
     ratioc = 1.0 - ratio;
@@ -1427,11 +1578,16 @@ void animation()
     params.radius = MU;
     
     particle = (t_particle *)malloc(NMAXCIRCLES*sizeof(t_particle));    /* particles */ 
+    
+    if (CLUSTER_PARTICLES) 
+        cluster = (t_cluster *)malloc(NMAXCIRCLES*sizeof(t_cluster));    /* particle clusters */ 
+    
     if (ADD_FIXED_OBSTACLES) obstacle = (t_obstacle *)malloc(NMAXOBSTACLES*sizeof(t_obstacle));    /* circular obstacles */  
     if (ADD_FIXED_SEGMENTS) 
     {
         segment = (t_segment *)malloc(NMAXSEGMENTS*sizeof(t_segment));         /* linear obstacles */  
         segment_group = (t_group_segments *)malloc(NMAXGROUPS*sizeof(t_group_segments));
+        conveyor_belt = (t_belt *)malloc(NMAXBELTS*sizeof(t_belt));
     }
     
     if (TRACER_PARTICLE) 
@@ -1450,6 +1606,16 @@ void animation()
     pangle = (double *)malloc(NMAXCIRCLES*sizeof(double));
     pressure = (double *)malloc(N_PRESSURES*sizeof(double));
     
+    if (CLUSTER_PARTICLES)
+    {
+        cqx = (double *)malloc(NMAXCIRCLES*sizeof(double));
+        cqy = (double *)malloc(NMAXCIRCLES*sizeof(double));
+        cpx = (double *)malloc(NMAXCIRCLES*sizeof(double));
+        cpy = (double *)malloc(NMAXCIRCLES*sizeof(double));
+        cqangle = (double *)malloc(NMAXCIRCLES*sizeof(double));
+        cpangle = (double *)malloc(NMAXCIRCLES*sizeof(double));
+    }
+    
     if (REACTION_DIFFUSION) 
     {
         collisions = (t_collision *)malloc(2*NMAXCOLLISIONS*sizeof(t_collision));
@@ -1465,7 +1631,7 @@ void animation()
     lj_log = fopen("lj_logfile.txt", "w");
     
     if (ADD_FIXED_OBSTACLES) init_obstacle_config(obstacle);
-    if (ADD_FIXED_SEGMENTS) init_segment_config(segment);
+    if (ADD_FIXED_SEGMENTS) init_segment_config(segment, conveyor_belt);
     
     if ((MOVE_SEGMENT_GROUPS)&&(ADD_FIXED_SEGMENTS)) 
     {
@@ -1504,13 +1670,15 @@ void animation()
     
     printf("%i active particles\n", params.nactive);
     
+    if (CLUSTER_PARTICLES) init_cluster_config(particle, cluster);
+     
 //     xi = 0.0;
     
 //     for (i=0; i<ncircles; i++)
 //     {
 //         printf("Particle %i at (%.3f, %.3f) of energy %.3f\n", i, particle[i].xc, particle[i].yc, particle[i].energy);
 //     }
-    sleep(1);
+//     sleep(1);
         
     update_hashgrid(particle, hashgrid, 1);
     printf("Updated hashgrid\n");
@@ -1695,6 +1863,29 @@ void animation()
                         }
                     }
                 }
+                /* TEST: average angles for clustered polygons */
+//                 if ((REACTION_DIFFUSION)&&(RD_REACTION == CHEM_POLYGON_AGGREGATION))
+//                 {
+//                     np = particle[j].npartners;
+//                     alpha = DPI/(double)NPOLY;
+//                     if (np >= 2)
+//                     {
+//                         angle = particle[j].angle;
+//                         while (angle > alpha) angle -= alpha;
+//                         sum = angle;
+//                         for (p=0; p<np; p++)
+//                         {
+//                             angle =  particle[particle[j].partner[p]].angle;
+//                             while (angle > alpha) angle -= alpha;
+//                             sum += angle;
+//                         }
+//                         sum *= 1.0/(double)(np+1);
+//                         particle[j].angle =  sum;
+//                         for (p=0; p<np; p++)
+//                             particle[particle[j].partner[p]].angle = sum;
+//                     }
+//                     
+//                 }
                 
                 /* add gravity */
                 if (INCREASE_GRAVITY) 
@@ -1742,10 +1933,29 @@ void animation()
                 }
             }
             
-            /* timestep of thermostat algorithm */
+            /* compute force and torque on clusters */
-            totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, params.beta, &params.nactive, &nsuccess, &nmove, &ncoupled, i < INITIAL_TIME);
+            if (CLUSTER_PARTICLES) compute_cluster_force(cluster, particle);
             
             
+            
+            /* timestep of thermostat algorithm */
+            if (CLUSTER_PARTICLES) 
+            {
+                totalenergy = evolve_clusters(particle, cluster, hashgrid, cqx, cqy, cqangle, cpx, cpy, cpangle, params.beta, &params.nactive, &nsuccess, &nmove, &ncoupled, i < INITIAL_TIME, n == 0);
+                
+                /* FIXME: update particle data */
+                for (j=0; j<ncircles; j++)
+                    particle[j].emean = cluster[particle[j].cluster].emean;
+            }
+            else 
+                totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, params.beta, &params.nactive, &nsuccess, &nmove, &ncoupled, i < INITIAL_TIME);
+            
+            /* TEST */
+            /* repair clusters */
+            if ((CLUSTER_PARTICLES)&&(REPAIR_CLUSTERS)) for (cl=0; cl<ncircles; cl++) 
+                if ((cluster[cl].active)&&(cluster[cl].nparticles >= 2))
+                    repair_cluster(cl, particle, cluster, 1000/NVID, 1);
+            
             /* evolution of lid coordinate */
             if (BOUNDARY_COND == BC_RECTANGLE_LID) evolve_lid(params.fboundary);
             if (BOUNDARY_COND == BC_RECTANGLE_WALL) 
@@ -1755,11 +1965,24 @@ void animation()
             }
             if ((MOVE_BOUNDARY)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segments(segment, i);
             
+            if (ADD_CONVEYOR_FORCE) for (belt = 0; belt < nbelts; belt++) 
+                conveyor_belt[belt].position += conveyor_belt[belt].speed*DT_PARTICLE;
+            
             if ((MOVE_SEGMENT_GROUPS)&&(i > INITIAL_TIME + SEGMENT_DEACTIVATION_TIME)) evolve_segment_groups(segment, i, segment_group);
             
 //             if ((MOVE_SEGMENT_GROUPS)&&(i > OBSTACLE_INITIAL_TIME)) evolve_segment_groups(segment, i, segment_group);
         } /* end of for (n=0; n<NVID; n++) */
         
+        /* TEST */
+//         for (j=0; j<ncircles; j++)
+//             {
+//                 if (particle[j].active) 
+//                     printf("Particle %i: (x,y) = (%.3lg, %.3lg), F = (%.3lg, %.3lg)\n", j, particle[j].xc, particle[j].yc, particle[j].fx, particle[j].fy); 
+//                 if (cluster[j].active) 
+//                     printf("Cluster %i: (x,y) = (%.3lg, %.3lg), F = (%.3lg, %.3lg)\n", j, cluster[j].xg, cluster[j].yg, cluster[j].fx, cluster[j].fy); 
+//             }
+//             sleep(1);
+        
         if ((i>INITIAL_TIME)&&(SAVE_TIME_SERIES))
         {
             n_total_active = 0;
@@ -1881,7 +2104,7 @@ void animation()
         /* case of reaction-diffusion equation */
         if ((i > INITIAL_TIME)&&(REACTION_DIFFUSION)) 
         {
-            ncollisions = update_types(particle, molecule, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1,  &delta_energy);
+            ncollisions = update_types(particle, molecule, cluster, collisions, ncollisions, particle_numbers, i - INITIAL_TIME - 1,  &delta_energy);
             if (EXOTHERMIC) params.beta *= 1.0/(1.0 + delta_energy/totalenergy);
             params.nactive = 0;
             for (j=0; j<ncircles; j++) if (particle[j].active)
@@ -1904,8 +2127,8 @@ void animation()
         }
  
         if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-        draw_particles(particle, PLOT, params.beta, collisions, ncollisions, BG_COLOR, hashgrid, params);
+        draw_particles(particle, cluster, PLOT, params.beta, collisions, ncollisions, BG_COLOR, hashgrid, params);
-        draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
+        draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, conveyor_belt, wall);
 
         /* add a particle */
         if ((ADD_PARTICLES)&&(i > ADD_TIME)&&((i - INITIAL_TIME - ADD_TIME)%ADD_PERIOD == 1)&&(i < NSTEPS - FINAL_NOADD_PERIOD))
@@ -1956,8 +2179,8 @@ void animation()
             count_particle_number(particle, particle_numbers, i - INITIAL_TIME);
         
         draw_frame(i, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, 
-        wall, pressure, pleft, pright, particle_numbers, 1, params, particle, 
+        wall, pressure, pleft, pright, particle_numbers, 1, params, particle, cluster, 
-        collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group);
+        collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt);
         
 	if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
 
@@ -1988,8 +2211,8 @@ void animation()
             if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
             {
                 draw_frame(i, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, 
-                wall, pressure, pleft, pright, particle_numbers, 0, params, particle, 
+                wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, 
-                collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group);
+                collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt);
                 glutSwapBuffers();
                 save_frame_lj_counter(NSTEPS + MID_FRAMES + 1 + counter);
                 counter++;
@@ -2029,8 +2252,8 @@ void animation()
         {
             blank();
             draw_frame(NSTEPS, PLOT, BG_COLOR, ncollisions, traj_position, traj_length, 
-            wall, pressure, pleft, pright, particle_numbers, 0, params, particle, 
+            wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, 
-            collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group);
+            collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt);
         }
         if (DOUBLE_MOVIE) for (i=0; i<MID_FRAMES; i++) 
         {
@@ -2041,8 +2264,8 @@ void animation()
         if (DOUBLE_MOVIE) 
         {
             draw_frame(NSTEPS, PLOT_B, BG_COLOR_B, ncollisions, traj_position, traj_length, 
-            wall, pressure, pleft, pright, particle_numbers, 0, params, particle, 
+            wall, pressure, pleft, pright, particle_numbers, 0, params, particle, cluster, 
-            collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group);
+            collisions, hashgrid, trajectory, obstacle, segment, group_speeds, segment_group, conveyor_belt);
             if (!((NO_EXTRA_BUFFER_SWAP)&&(MOVIE))) glutSwapBuffers();
         }
         if ((TIME_LAPSE)&&(!DOUBLE_MOVIE)) 
@@ -2064,6 +2287,8 @@ void animation()
     
 //     printf("1\n");
     free(particle);
+    
+    if (CLUSTER_PARTICLES) free(cluster);
 //     printf("2\n");
     if (ADD_FIXED_OBSTACLES) free(obstacle);
 //     printf("3\n");
@@ -2071,6 +2296,7 @@ void animation()
     {
         free(segment);
         free(segment_group);
+        free(conveyor_belt);
     }
 //     printf("4\n");
     if (MOVE_SEGMENT_GROUPS) free(group_speeds);
@@ -2095,6 +2321,17 @@ void animation()
 //     printf("14\n");
     free(pangle);
 //     printf("15\n");
+    
+    if (CLUSTER_PARTICLES)
+    {
+        free(cqx);
+        free(cqy);
+        free(cpx);
+        free(cpy);
+        free(cqangle);
+        free(cpangle);
+    }
+    
     free(pressure);
 //     printf("16\n");
     if (REACTION_DIFFUSION) free(collisions);

mangrove.c

@@ -71,6 +71,7 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 340        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.25   /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 0.85	    /* parameter controlling the dimensions of domain */
@@ -177,6 +178,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.0   /* additional shift for wave amplitude */
+#define VSCALE_AMPLITUDE 0.5   /* additional scaling factor for wave amplitude */
 #define E_SCALE 2500.0     /* scaling factor for energy representation */
 #define LOG_SCALE 1.0    /* scaling factor for energy log representation */
 #define LOG_SHIFT 0.0     /* shift of colors on log scale */
@@ -247,6 +250,7 @@
 #define POT_MAZE 7
 #define POTENTIAL 0
 #define MAZE_WIDTH 0.02     /* half width of maze walls */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 #define VARIABLE_IOR 0      /* set to 1 for a variable index of refraction */
 #define IOR 7               /* choice of index of refraction, see list in global_pdes.c */
 #define IOR_TOTAL_TURNS 1.5 /* total angle of rotation for IOR_PERIODIC_WELLS_ROTATING */
@@ -283,6 +287,7 @@
 #define SHADE_2D 0       /* set to 1 to add pseudo-3d shading effect */ 
 #define SHADE_SCALE_2D 0.05  /* lower value increases sensitivity of shading */
 #define N_SOURCES 1                     /* number of sources, for option draw_sources */
+#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
 double light[2] = {0.40824829, 0.816496581};   /* location of light source for SHADE_2D option*/
 /* end of constants only used by sub_wave and sub_maze */
 

rde.c

@@ -39,7 +39,7 @@
 #include <omp.h>
 #include <time.h>
 
-#define MOVIE 1         /* set to 1 to generate movie */
+#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 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,13 +48,9 @@
 
 #define WINWIDTH 	1920  /* window width */
 #define WINHEIGHT 	1150  /* window height */
-// #define NX 240          /* number of grid points on x axis */
+#define NX 960          /* number of grid points on x axis */
-// #define NY 120          /* number of grid points on y axis */
+#define NY 575          /* number of grid points on y axis */
-// #define NX 960          /* number of grid points on x axis */
+#define HRES 1          /* factor for high resolution plots */
-// #define NY 500          /* number of grid points on y axis */
-#define NX 780          /* number of grid points on x axis */
-#define NY 400          /* number of grid points on y axis */
-#define HRES 2          /* factor for high resolution plots */
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
@@ -63,31 +59,32 @@
 
 /* Choice of simulated equation */
 
-#define RDE_EQUATION 7  /* choice of reaction term, see list in global_3d.c */
+#define RDE_EQUATION 8  /* 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 */
 
 #define SPHERE 1        /* set to 1 to simulate equation on sphere */
 #define DPOLE 1         /* safety distance to poles */
 #define DSMOOTH 1       /* size of neighbourhood of poles that are smoothed */
-#define SMOOTHPOLE 0.01  /* smoothing coefficient at poles */
+#define SMOOTHPOLE 0.05  /* smoothing coefficient at poles */
-#define SMOOTHCOTPOLE 0.01  /* smoothing coefficient of cotangent at poles */
+#define SMOOTHCOTPOLE 0.05  /* smoothing coefficient of cotangent at poles */
 #define PHISHIFT 0.0    /* shift of phi in 2D plot (in degrees) */
 #define SMOOTHBLOCKS 1  /* set to 1 to use blocks of points near the poles */
 #define BLOCKDIST 64    /* distance to poles where points are blocked */ 
 #define ZERO_MERIDIAN 190.0     /* choice of zero meridian (will be at left/right boundary of 2d plot) */
+#define POLE_NODRAW 10  /* distance around poles where wave is not drawn */
 
 #define ADD_POTENTIAL 0 /* set to 1 to add a potential (for Schrodinger equation) */
 #define ADD_MAGNETIC_FIELD 0    /* set to 1 to add a magnetic field (for Schrodinger equation) - then set POTENTIAL 1 */
-#define ADD_FORCE_FIELD 1   /* set to 1 to add a foce field (for compressible Euler equation) */
+#define ADD_FORCE_FIELD 0   /* set to 1 to add a foce field (for compressible Euler equation) */
 #define POTENTIAL 7         /* type of potential or vector potential, see list in global_3d.c  */
 #define FORCE_FIELD 6       /* type of force field, see list in global_3d.c  */
 #define ADD_CORIOLIS_FORCE 1    /* set to 1 to add Coriolis force (quasigeostrophic Euler equations) */
-#define VARIABLE_DEPTH 0    /* set to 1 for variable depth in shallow water equation */
+#define VARIABLE_DEPTH 1    /* set to 1 for variable depth in shallow water equation */
-#define SWATER_DEPTH 4      /* variable depth in shallow water equation */
+#define SWATER_DEPTH 5      /* variable depth in shallow water equation */
 
 #define ANTISYMMETRIZE_WAVE_FCT 0   /* set tot 1 to make wave function antisymmetric */
-#define ADAPT_STATE_TO_BC 1     /* to smoothly adapt initial state to obstacles */
+#define ADAPT_STATE_TO_BC 0     /* to smoothly adapt initial state to obstacles */
 #define OBSTACLE_GEOMETRY 84     /* geometry of obstacles, as in B_DOMAIN */
 #define BC_STIFFNESS 50.0        /* controls region of boundary condition control */
 #define CHECK_INTEGRAL 1     /* set to 1 to check integral of first field */
@@ -105,6 +102,7 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 300        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.25    /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 1.0	    /* parameter controlling the dimensions of domain */
@@ -120,6 +118,7 @@
 #define NGRIDY 8            /* number of grid point for grid of disks */
 #define REVERSE_TESLA_VALVE 1   /* set to 1 to orient Tesla valve in blocking configuration */
 #define WALL_WIDTH 0.05      /* width of wall separating lenses */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 
 #define X_SHOOTER -0.2
 #define Y_SHOOTER -0.6
@@ -137,17 +136,11 @@
 
 /* Physical parameters of wave equation */
 
-// #define DT 0.0000001
+#define DT 0.00000025
-#define DT 0.00000015
-// #define DT 0.0000002
-// #define DT 0.0000012
-// #define DT 0.000003
-// #define DT 0.0000022
-// #define DT 0.0000012
 
 #define VISCOSITY 0.02
 #define POISSON_STIFFNESS 1.0   /* stiffness of Poisson equation solver for incompressible Euler */
-#define DISSIPATION 0.0
+#define DISSIPATION 1.0e-8
 
 #define RPSA 0.75         /* parameter in Rock-Paper-Scissors-type interaction */
 #define RPSLZB 0.0       /* second parameter in Rock-Paper-Scissors-Lizard-Spock type interaction */
@@ -163,7 +156,7 @@
 #define BZQ 0.0008      /* parameter in BZ equation */
 #define BZF 1.2         /* parameter in BZ equation */
 #define B_FIELD 10.0    /* magnetic field */
-#define G_FIELD 0.03    /* gravity/Coriolis force */
+#define G_FIELD 0.002   /* gravity/Coriolis force */
 #define BC_FIELD 1.0e-5     /* constant in repulsive field from obstacles */
 #define AB_RADIUS 0.2   /* radius of region with magnetic field for Aharonov-Bohm effect */
 #define K_EULER 50.0    /* constant in stream function integration of Euler equation */
@@ -172,8 +165,8 @@
 
 #define SMOOTHEN_VORTICITY 0    /* set to 1 to smoothen vorticity field in Euler equation */
 #define SMOOTHEN_VELOCITY 1     /* set to 1 to smoothen velocity field in Euler equation */
-#define SMOOTHEN_PERIOD 10      /* period between smoothenings */
+#define SMOOTHEN_PERIOD 7       /* period between smoothenings */
-#define SMOOTH_FACTOR 0.05      /* factor by which to smoothen */
+#define SMOOTH_FACTOR 0.15      /* factor by which to smoothen */
 
 #define ADD_OSCILLATING_SOURCE 0        /* set to 1 to add an oscillating wave source */
 #define OSCILLATING_SOURCE_PERIOD 1     /* period of oscillating source */
@@ -181,7 +174,7 @@
 
 #define ADD_TRACERS 1    /* set to 1 to add tracer particles (for Euler equations) */
 #define N_TRACERS 2000    /* number of tracer particles */
-#define TRACERS_STEP 0.005  /* step size in tracer evolution */
+#define TRACERS_STEP 0.1  /* step size in tracer evolution */
 
 #define T_OUT 2.0       /* outside temperature */
 #define T_IN 0.0        /* inside temperature */
@@ -215,8 +208,8 @@
 #define LAMINAR_FLOW_YPERIOD 1.0    /* period of laminar flow in y direction */
 #define PRESSURE_GRADIENT 0.2       /* amplitude of pressure gradient for Euler equation */
 
-#define SWATER_MIN_HEIGHT 0.5      /* min height of initial condition for shallow water equation */
+#define SWATER_MIN_HEIGHT 0.025      /* min height of initial condition for shallow water equation */
-#define DEPTH_FACTOR 0.015        /* proportion of min height in variable depth */
+#define DEPTH_FACTOR 0.075       /* proportion of min height in variable depth */
 #define TANH_FACTOR 1.0         /* steepness of variable depth */
 
 #define EULER_GRADIENT_YSHIFT 0.0    /* y-shift in computation of gradient in Euler equation */
@@ -232,9 +225,8 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 1900           /* number of frames of movie */
+#define NSTEPS 950           /* number of frames of movie */
-// #define NSTEPS 500           /* number of frames of movie */
+#define NVID 85           /* number of iterations between images displayed on screen */
-#define NVID 100          /* 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 */
@@ -256,24 +248,23 @@
 #define PLOT_SPHERE 1   /* draws fields on a sphere */
 
 #define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
-#define ROTATE_ANGLE 360.0  /* total angle of rotation during simulation */
+#define ROTATE_ANGLE 360.0   /* total angle of rotation during simulation */
-// #define ROTATE_ANGLE 90.0  /* total angle of rotation during simulation */
 #define SHADE_3D 1              /* set to 1 to change luminosity according to normal vector */
 #define SHADE_2D 0              /* set to 1 to change luminosity according to normal vector */
-#define VIEWPOINT_TRAJ 0    /* type of viewpoint trajectory */
+#define VIEWPOINT_TRAJ 1    /* type of viewpoint trajectory */
 #define MAX_LATITUDE 45.0   /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
 
 #define DRAW_PERIODICISED 0     /* set to 1 to repeat wave periodically in x and y directions */
 
 /* Plot type - color scheme */
 
-#define CPLOT 62
+#define CPLOT 70
-#define CPLOT_B 64
+#define CPLOT_B 74
 
 /* Plot type - height of 3D plot */
 
-#define ZPLOT 62     /* z coordinate in 3D plot */
+#define ZPLOT 70     /* z coordinate in 3D plot */
-#define ZPLOT_B 64    /* z coordinate in second 3D plot */
+#define ZPLOT_B 71    /* 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 */
@@ -283,7 +274,6 @@
 #define DRAW_OUTSIDE_GRAY 0     /* experimental - draw outside of billiard in gray */
 #define ADD_POTENTIAL_TO_Z 0    /* set to 1 to add the external potential to z-coordinate of plot */
 #define ADD_POT_CONSTANT 0.35   /* constant added to wave height */
-#define ADD_HEIGHT_CONSTANT -0.5   /* constant added to wave height */
 #define DRAW_DEPTH 0            /* set to 1 to draw water depth */
 #define DEPTH_SCALE 0.75         /* vertical scaling of depth plot */
 #define DEPTH_SHIFT -0.015      /* vertical shift of depth plot */
@@ -317,7 +307,7 @@
 /* Color schemes, see list in global_pdes.c  */
 
 #define COLOR_PALETTE 10       /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 17     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 12      /* Color palette, see list in global_pdes.c  */
 
 #define BLACK 1          /* black background */
 #define COLOR_OUT_R 1.0    /* color outside domain */
@@ -326,20 +316,22 @@
 
 #define COLOR_SCHEME 3   /* choice of color scheme */
 
-#define COLOR_PHASE_SHIFT 0.25   /* phase shift of color scheme, in units of Pi */
+#define PHASE_SHIFT -0.25   /* phase shift of color scheme, in units of Pi (formerly COLOR_PHASE_SHIFT) */
 
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
-#define VSCALE_AMPLITUDE 100.0      /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define VSHIFT_AMPLITUDE 0.0   /* additional shift for wave amplitude */
+#define VSCALE_AMPLITUDE 15.0      /* additional scaling factor for color scheme P_3D_AMPLITUDE */
 #define ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
 #define CURL_SCALE 1.0   /* 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) */
-#define SLOPE_SCHROD_LUM 100.0       /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
+#define SLOPE_SCHROD_LUM 10.0      /* sensitivity of luminosity on module, for color scheme Z_ARGUMENT */
 #define MIN_SCHROD_LUM 0.1       /* minimal luminosity in color scheme Z_ARGUMENT*/
 #define VSCALE_PRESSURE 2.0      /* additional scaling factor for color scheme Z_EULER_PRESSURE */
 #define PRESSURE_SHIFT 10.0        /* shift for color scheme Z_EULER_PRESSURE */
 #define PRESSURE_LOG_SHIFT -2.5     /* shift for color scheme Z_EULER_PRESSURE */
-#define VSCALE_WATER_HEIGHT 0.4     /* vertical scaling of water height */
+#define VSCALE_WATER_HEIGHT 30.0     /* vertical scaling of water height */
+#define ADD_HEIGHT_CONSTANT -0.025   /* constant added to wave height */
 #define SHADE_SCALE_2D 0.25     /* controls "depth" of 2D shading */
 
 #define COLORHUE 260     /* initial hue of water color for scheme C_LUM */
@@ -359,10 +351,10 @@
 #define VSCALE_DENSITY 30.0      /* additional scaling factor for color scheme Z_EULER_DENSITY */
 #define VSCALE_VORTICITY 15.0     /* additional scaling factor for color scheme Z_EULERC_VORTICITY */
 #define VORTICITY_SHIFT 0.0     /* vertical shift of vorticity */
-#define ZSCALE_SPEED 0.5        /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
+#define ZSCALE_SPEED 300.0        /* additional scaling factor for z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
 #define ZSHIFT_SPEED 0.0        /* additional shift of z-coord Z_EULER_SPEED and Z_SWATER_SPEED */
 #define ZSCALE_NORMGRADIENT -0.0001  /* vertical scaling for Z_NORM_GRADIENT */
-#define VSCALE_SWATER 250.0        /* additional scaling factor for color scheme Z_EULER_DENSITY */
+#define VSCALE_SWATER 50.0        /* additional scaling factor for color scheme Z_EULER_DENSITY */
 
 #define NXMAZE 7      /* width of maze */
 #define NYMAZE 7      /* height of maze */
@@ -372,7 +364,7 @@
 #define MAZE_WIDTH 0.04     /* half width of maze walls */
 
 #define DRAW_COLOR_SCHEME 1     /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 3.0      /* scale of color scheme bar */
+#define COLORBAR_RANGE 2.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 */
@@ -393,7 +385,6 @@
 #define INITIAL_WAVELENGTH  0.1  /* wavelength of initial condition */
 #define VSCALE_ENERGY 200.0       /* additional scaling factor for color scheme P_3D_ENERGY */
 #define PHASE_FACTOR 20.0       /* factor in computation of phase in color scheme P_3D_PHASE */
-#define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
 #define AVRG_E_FACTOR 0.99   /* controls time window size in P_AVERAGE_ENERGY scheme */
 #define OSCILLATION_SCHEDULE 0  /* oscillation schedule, see list in global_pdes.c */
 #define AMPLITUDE 0.8      /* amplitude of periodic excitation */ 
@@ -436,8 +427,8 @@
 double u_3d[2] = {0.75, -0.45};     /* projections of basis vectors for REP_AXO_3D representation */
 double v_3d[2] = {-0.75, -0.45};
 double w_3d[2] = {0.0, 0.015};
-double light[3] = {-0.40824829, -0.816496581, 0.40824829};      /* vector of "light" direction for P_3D_ANGLE color scheme */
+double light[3] = {0.40824829, -0.816496581, 0.40824829};      /* vector of "light" direction for P_3D_ANGLE color scheme */
-double observer[3] = {-8.0, -4.0, 4.0};    /* location of observer for REP_PROJ_3D representation */ 
+double observer[3] = {8.0, 4.0, 2.5};    /* location of observer for REP_PROJ_3D representation */ 
 int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
 
 /* constants for simulations on planets */
@@ -449,7 +440,7 @@ int reset_view = 0;         /* switch to reset 3D view parameters (for option RO
 #define DEM_SMOOTH_HEIGHT 2.0   /* relative height below which to smoothen */
 #define DEM_MAXHEIGHT 9000.0     /* max height of DEM (estimated from Everest/Olympus Mons) */
 #define DEM_MAXDEPTH -10000     /* max depth of DEM */
-#define PLANET_SEALEVEL 0.0      /* sea level for flooded planet */
+#define PLANET_SEALEVEL 2500.0      /* sea level for flooded planet */
 #define VENUS_NODATA_FACTOR 0.5     /* altitude to assign to DEM points without data (fraction of mean altitude) */
 
 #define Z_SCALING_FACTOR 0.8   /* overall scaling factor of z axis for REP_PROJ_3D representation */
@@ -458,27 +449,27 @@ int reset_view = 0;         /* switch to reset 3D view parameters (for option RO
 #define XSHIFT_3D 0.0          /* overall x shift for REP_PROJ_3D representation */
 #define YSHIFT_3D 0.0          /* overall y shift for REP_PROJ_3D representation */
 #define BORDER_PADDING 0       /* distance from boundary at which to plot points, to avoid boundary effects due to gradient */
-#define DRAW_ARROW 0           /* set to 1 to draw arrow above sphere */
+#define DRAW_ARROW 1           /* set to 1 to draw arrow above sphere */
-
-#define RSCALE 0.01            /* scaling factor of radial component */
-#define RSHIFT -0.01           /* shift in radial component */
-#define RMAX 2.0               /* max value of radial component */
-#define RMIN 0.5               /* min value of radial component */
-// #define COS_VISIBLE -1.1       /* limit on cosine of normal to shown facets */
-#define COS_VISIBLE -0.3        /* limit on cosine of normal to shown facets */
 
+#define RSCALE 0.2               /* scaling factor of radial component */
+#define RSHIFT -0.01             /* shift in radial component */
+#define RMAX 2.0                 /* max value of radial component */
+#define RMIN 0.5                 /* min value of radial component */
+#define COS_VISIBLE -0.75         /* limit on cosine of normal to shown facets */
 
 /* 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 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)))
 #define PRINT_PARAMETERS ((PRINT_TIME)||(PRINT_VISCOSITY)||(PRINT_RPSLZB)||(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)
+#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
 
 int block_sizes[NY];      /* table of block sizes for blocking around poles */
 int block_numbers[NY];    /* table of block numbers for blocking around poles */
@@ -711,83 +702,165 @@ void initialize_gfield(double gfield[2*NX*NY], double bc_field[NX*NY], double bc
     int i, j;
     double dx, dy;
         
-    if (FORCE_FIELD == GF_COMPUTE_FROM_BC)
+    switch (FORCE_FIELD) 
     {
-        dx = (XMAX - XMIN)/(double)NX;
+        case (GF_COMPUTE_FROM_BC):
-        dy = (YMAX - YMIN)/(double)NY;
+        {
+            dx = (XMAX - XMIN)/(double)NX;
+            dy = (YMAX - YMIN)/(double)NY;
 
-        #pragma omp parallel for private(i,j)
+            #pragma omp parallel for private(i,j)
-        for (i=1; i<NX-1; i++){
+            for (i=1; i<NX-1; i++){
-            for (j=1; j<NY-1; j++){
+                for (j=1; j<NY-1; j++){
-                gfield[i*NY+j] = BC_FIELD*(bc_field2[(i+1)*NY+j] - bc_field2[(i-1)*NY+j])/dx;
+                    gfield[i*NY+j] = BC_FIELD*(bc_field2[(i+1)*NY+j] - bc_field2[(i-1)*NY+j])/dx;
-                gfield[NX*NY+i*NY+j] = BC_FIELD*(bc_field2[i*NY+j+1] - bc_field2[i*NY+j-1])/dy;
+                    gfield[NX*NY+i*NY+j] = BC_FIELD*(bc_field2[i*NY+j+1] - bc_field2[i*NY+j-1])/dy;
 //                 printf("gfield at (%i,%i): (%.3lg, %.3lg)\n", i, j, gfield[i*NY+j], gfield[NX*NY+i*NY+j]);
+                }
             }
-        }
         
-        /* boundaries */
+            /* boundaries */
-        for (i=0; i<NX; i++)
+            for (i=0; i<NX; i++)
-        {
+            {
-            gfield[i*NY] = 0.0;
+                gfield[i*NY] = 0.0;
-            gfield[NX*NY+i*NY] = 0.0;
+                gfield[NX*NY+i*NY] = 0.0;
-            gfield[i*NY+NY-1] = 0.0;
+                gfield[i*NY+NY-1] = 0.0;
-            gfield[NX*NY+i*NY+NY-1] = 0.0;
+                gfield[NX*NY+i*NY+NY-1] = 0.0;
-        }
-        for (j=0; j<NY; j++)
-        {
-            gfield[j] = 0.0;
-            gfield[NX*NY+j] = 0.0;
-            gfield[(NX-1)*NY+j] = 0.0;
-            gfield[NX*NY+(NX-1)*NY+j] = 0.0;            
-        }
-    }
-    else if (FORCE_FIELD == GF_EARTH)
-    {
-        dx = (XMAX - XMIN)/(double)NX;
-        dy = (YMAX - YMIN)/(double)NY;
-        init_earth_map_rde(wsphere, 1);
-        init_earth_map_rde(wsphere_hr, HRES);
-        
-        #pragma omp parallel for private(i,j)
-        for (i=1; i<NX-1; i++){
-            for (j=1; j<NY-1; j++){
-                gfield[i*NY+j] = BC_FIELD*(wsphere[(i+1)*NY+j].altitude - wsphere[(i-1)*NY+j].altitude)/dx;
-                gfield[NX*NY+i*NY+j] = BC_FIELD*(wsphere[i*NY+j+1].altitude - wsphere[i*NY+j-1].altitude)/dy;
             }
+            for (j=0; j<NY; j++)
+            {
+                gfield[j] = 0.0;
+                gfield[NX*NY+j] = 0.0;
+                gfield[(NX-1)*NY+j] = 0.0;
+                gfield[NX*NY+(NX-1)*NY+j] = 0.0;            
+            }
+            break;
         }
+        case (GF_EARTH):
+        {
+            dx = (XMAX - XMIN)/(double)NX;
+            dy = (YMAX - YMIN)/(double)NY;
+            init_earth_map_rde(wsphere, 1);
+            init_earth_map_rde(wsphere_hr, HRES);
         
-        /* boundaries TODO */
+            #pragma omp parallel for private(i,j)
-        for (i=0; i<NX; i++)
+            for (i=1; i<NX-1; i++){
-        {
+                for (j=1; j<NY-1; j++){
-            gfield[i*NY] = 0.0;
+                    gfield[i*NY+j] = BC_FIELD*(wsphere[(i+1)*NY+j].altitude - wsphere[(i-1)*NY+j].altitude)/dx;
-            gfield[NX*NY+i*NY] = 0.0;
+                    gfield[NX*NY+i*NY+j] = BC_FIELD*(wsphere[i*NY+j+1].altitude - wsphere[i*NY+j-1].altitude)/dy;
-            gfield[i*NY+NY-1] = 0.0;
+                }
-            gfield[NX*NY+i*NY+NY-1] = 0.0;
+            }
+        
+            /* boundaries TODO */
+            for (i=0; i<NX; i++)
+            {
+                gfield[i*NY] = 0.0;
+                gfield[NX*NY+i*NY] = 0.0;
+                gfield[i*NY+NY-1] = 0.0;
+                gfield[NX*NY+i*NY+NY-1] = 0.0;
+            }
+            for (j=0; j<NY; j++)
+            {
+                gfield[j] = 0.0;
+                gfield[NX*NY+j] = 0.0;
+                gfield[(NX-1)*NY+j] = 0.0;
+                gfield[NX*NY+(NX-1)*NY+j] = 0.0;            
+            }
+            break;
         }
-        for (j=0; j<NY; j++)
+        case (GF_MARS):
         {
-            gfield[j] = 0.0;
+            dx = (XMAX - XMIN)/(double)NX;
-            gfield[NX*NY+j] = 0.0;
+            dy = (YMAX - YMIN)/(double)NY;
-            gfield[(NX-1)*NY+j] = 0.0;
+            init_planet_map_rde(wsphere, D_SPHERE_MARS, 1);
-            gfield[NX*NY+(NX-1)*NY+j] = 0.0;            
+            init_planet_map_rde(wsphere_hr, D_SPHERE_MARS, HRES);
+        
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<NX-1; i++){
+                for (j=1; j<NY-1; j++){
+                    gfield[i*NY+j] = BC_FIELD*(wsphere[(i+1)*NY+j].altitude - wsphere[(i-1)*NY+j].altitude)/dx;
+                    gfield[NX*NY+i*NY+j] = BC_FIELD*(wsphere[i*NY+j+1].altitude - wsphere[i*NY+j-1].altitude)/dy;
+                }
+            }
+        
+            /* boundaries TODO */
+            for (i=0; i<NX; i++)
+            {
+                gfield[i*NY] = 0.0;
+                gfield[NX*NY+i*NY] = 0.0;
+                gfield[i*NY+NY-1] = 0.0;
+                gfield[NX*NY+i*NY+NY-1] = 0.0;
+            }
+            for (j=0; j<NY; j++)
+            {
+                gfield[j] = 0.0;
+                gfield[NX*NY+j] = 0.0;
+                gfield[(NX-1)*NY+j] = 0.0;
+                gfield[NX*NY+(NX-1)*NY+j] = 0.0;            
+            }
+            break;
         }
-    }
+        case (GF_VENUS):
-    else
+        {
-    {
+            dx = (XMAX - XMIN)/(double)NX;
-        #pragma omp parallel for private(i,j)
+            dy = (YMAX - YMIN)/(double)NY;
-        for (i=0; i<NX; i++){
+            init_planet_map_rde(wsphere, D_SPHERE_VENUS, 1);
-            for (j=0; j<NY; j++){
+            init_planet_map_rde(wsphere_hr, D_SPHERE_VENUS, HRES);
-                compute_gfield(i, j, &gfield[i*NY+j], &gfield[NX*NY+i*NY+j]);
+        
+            #pragma omp parallel for private(i,j)
+            for (i=1; i<NX-1; i++){
+                for (j=1; j<NY-1; j++){
+                    gfield[i*NY+j] = BC_FIELD*(wsphere[(i+1)*NY+j].altitude - wsphere[(i-1)*NY+j].altitude)/dx;
+                    gfield[NX*NY+i*NY+j] = BC_FIELD*(wsphere[i*NY+j+1].altitude - wsphere[i*NY+j-1].altitude)/dy;
+                }
+            }
+        
+            /* boundaries TODO */
+            for (i=0; i<NX; i++)
+            {
+                gfield[i*NY] = 0.0;
+                gfield[NX*NY+i*NY] = 0.0;
+                gfield[i*NY+NY-1] = 0.0;
+                gfield[NX*NY+i*NY+NY-1] = 0.0;
+            }
+            for (j=0; j<NY; j++)
+            {
+                gfield[j] = 0.0;
+                gfield[NX*NY+j] = 0.0;
+                gfield[(NX-1)*NY+j] = 0.0;
+                gfield[NX*NY+(NX-1)*NY+j] = 0.0;            
+            }
+            break;
+        }
+        default:
+        {
+            #pragma omp parallel for private(i,j)
+            for (i=0; i<NX; i++){
+                for (j=0; j<NY; j++){
+                    compute_gfield(i, j, &gfield[i*NY+j], &gfield[NX*NY+i*NY+j]);
+                }
             }
         }
     }
 }
 
-void compute_water_depth(int i, int j, t_rde *rde, int ncircles)
+void compute_water_depth(int i, int j, t_rde *rde, t_wave_sphere wsphere[NX*NY], int ncircles)
 /* initialize the vector potential, for Schrodinger equation in a magnetic field */
 {
-    double x, y, xy[2], r0, r, z, h, h0, hh, x1, y1, d;
+    double x, y, xy[2], r0, r, z, h, h0, hh, x1, y1, z1, d, rmax, rmin, phi;
     int n, nx, ny;
+    static double phi1, phi2, phi3, sq3, rico;
+    static int first = 1;
+    
+    if (first)
+    {
+        phi = 0.5*(sqrt(5.0)+1.0);
+        phi1 = phi - 1.0;
+        phi2 = phi - 2.0;
+        phi3 = 2.0*phi - 3.0;
+        sq3 = 0.5/sqrt(3.0);
+        rico = 0.5/(2.0 + phi);
+        first = 0;
+    }
     
     ij_to_xy(i, j, xy);
     x = xy[0];
@@ -843,10 +916,92 @@ void compute_water_depth(int i, int j, t_rde *rde, int ncircles)
             rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
             break;
         }
+        case (SH_SPHERE_CUBE):
+        {
+            rmax = 0.0;
+            /* compute distance from cube to origin, given by 0.5/rmax */
+            if (vabs(wsphere[i*NY+j].x) > rmax) rmax = vabs(wsphere[i*NY+j].x);
+            if (vabs(wsphere[i*NY+j].y) > rmax) rmax = vabs(wsphere[i*NY+j].y);
+            if (vabs(wsphere[i*NY+j].z) > rmax) rmax = vabs(wsphere[i*NY+j].z);
+            h = 1.0 - 0.5/rmax;
+//             printf("h = %.3lg\n", h);
+            if (h < 0.0) h = 0.0;
+            rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
+            break;
+        }
+        case (SH_SPHERE_OCTAHEDRON):
+        {
+            rmax = 0.0;
+            /* compute distance from octahedron to origin, given by 0.5/rmax */
+            rmax = 1.0/(vabs(wsphere[i*NY+j].x) + vabs(wsphere[i*NY+j].y) + vabs(wsphere[i*NY+j].z));
+            h = 1.0 - 0.5/rmax;
+            printf("h = %.3lg\n", h);
+            if (h < 0.0) h = 0.0;
+            rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
+            break;
+        }
+        case (SH_SPHERE_DODECAHEDRON):
+        {
+            rmax = 0.0;
+            x1 = wsphere[i*NY+j].x;
+            y1 = wsphere[i*NY+j].y;
+            z1 = wsphere[i*NY+j].z;
+            /* compute distance from dodecahedron */
+            d = vabs(phi1*y1 - phi2*z1) - sq3;
+            if (d > rmax) rmax = d; 
+            d = vabs(-phi2*x1 + phi1*z1) - sq3;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi1*x1 - phi2*y1) - sq3;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi1*y1 + phi2*z1) - sq3;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi2*x1 + phi1*z1) - sq3;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi1*x1 + phi2*y1) - sq3;
+            if (d > rmax) rmax = d; 
+            h = rmax;
+            printf("h = %.3lg\n", h);
+            if (h < 0.0) h = 0.0;
+            rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
+            break;
+        }
+        case (SH_SPHERE_ICOSAHEDRON):
+        {
+            rmax = 0.0;
+            x1 = wsphere[i*NY+j].x;
+            y1 = wsphere[i*NY+j].y;
+            z1 = wsphere[i*NY+j].z;
+            /* compute distance from dodecahedron */
+            d = vabs(phi2*(x1+y1+z1)) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi2*(-x1+y1+z1)) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi2*(x1-y1+z1)) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi2*(x1+y1-z1)) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi3*x1 + phi1*z1) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi3*z1 + phi1*y1) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi3*y1 + phi1*x1) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi3*x1 - phi1*z1) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi3*z1 - phi1*y1) - rico;
+            if (d > rmax) rmax = d; 
+            d = vabs(phi3*y1 - phi1*x1) - rico;
+            if (d > rmax) rmax = d; 
+            h = rmax;
+            printf("h = %.3lg\n", h);
+            if (h < 0.0) h = 0.0;
+            rde->depth = SWATER_MIN_HEIGHT*DEPTH_FACTOR*h;
+            break;
+        }
     }
 }
 
-double initialize_water_depth(t_rde rde[NX*NY])
+double initialize_water_depth(t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
 {
     int i, j, ncircles;
     double dx, dy, min, max, pscal, norm, vz = 0.01;
@@ -856,7 +1011,7 @@ double initialize_water_depth(t_rde rde[NX*NY])
     #pragma omp parallel for private(i,j)
     for (i=0; i<NX; i++){
         for (j=0; j<NY; j++){
-            compute_water_depth(i, j, &rde[i*NY + j], ncircles);
+            compute_water_depth(i, j, &rde[i*NY + j], wsphere, ncircles);
             if (rde[i*NY + j].depth < 0.0) rde[i*NY + j].depth = 0.0;
         }
     }
@@ -1328,9 +1483,30 @@ double gfield[2*NX*NY], t_rde rde[NX*NY], t_wave_sphere wsphere[NX*NY])
                     vx = rde[i*NY+j].dxv;
                     vy = rde[i*NY+j].dyv;
                     
-                    phi_out[0][i*NY+j] = eta - intstep*(u*etax + v*etay + eta*(ux + vy));
+//                     printf("etax = %.3lg, etay = %.3lg, ux = %.3lg, uy = %.3lg, vx = %.3lg, vy = %.3lg\n", etax, etay, ux, uy, vx, vy);
-                    phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax);
+                    
-                    phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etay);
+                    if (SPHERE)
+                    {
+                        /* TODO */
+                        if ((j > DSMOOTH)&&(j < NY-DSMOOTH))
+                        {
+                            phi_out[0][i*NY+j] = eta - intstep*(u*etax + v*etay + eta*(ux + vy));
+                            phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax);
+                            phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etax);
+                        }   
+                        else
+                        {
+                            phi_out[0][i*NY+j] = SWATER_MIN_HEIGHT;
+                            phi_out[1][i*NY+j] = 0.0;
+                            phi_out[2][i*NY+j] = 0.0;
+                        }
+                    }
+                    else
+                    {
+                        phi_out[0][i*NY+j] = eta - intstep*(u*etax + v*etay + eta*(ux + vy));
+                        phi_out[1][i*NY+j] = u - intstep*(u*ux + v*uy + G_FIELD*etax);
+                        phi_out[2][i*NY+j] = v - intstep*(u*vx + v*vy + G_FIELD*etay);
+                    }
                     
                     if (VISCOSITY > 0.0)
                     {
@@ -1513,6 +1689,12 @@ void evolve_tracers(double *phi[NFIELDS], double tracers[2*N_TRACERS*NSTEPS], t_
                     vy = phi[2][i*NY+j];
                     break;
                 }
+                case (E_SHALLOW_WATER):
+                {
+                    vx = phi[1][i*NY+j];
+                    vy = phi[2][i*NY+j];
+                    break;
+                }
             }
             
 //             v = module2(vx, vy);
@@ -1842,7 +2024,7 @@ void animation()
     if (VARIABLE_DEPTH) 
     {
 //         water_depth = (t_swater_depth *)malloc(NX*NY*sizeof(t_swater_depth));
-        max_depth = initialize_water_depth(rde);
+        max_depth = initialize_water_depth(rde, wsphere);
         printf("Max depth = %.3lg\n", max_depth);
     }
     
@@ -1883,7 +2065,8 @@ void animation()
     /* initialize field */
 //     init_random(0.5, 0.4, phi, xy_in);
 //     init_random(0.5, 0.25, phi, xy_in, wsphere);
-//     init_random_smoothed(0.5, 0.25, phi, xy_in, wsphere);
+//     init_random_smoothed(1.0, 0.1, phi, xy_in, wsphere);
+    
 //     init_gaussian(x, y, mean, amplitude, scalex, phi, xy_in)
 //     init_coherent_state(0.0, 0.0, 10.0, 0.0, 0.1, phi, xy_in);
     
@@ -1914,29 +2097,12 @@ void animation()
 
 //     init_laminar_flow(double amp, double xmodulation, double ymodulation, double xperiod, double yperiod, double yshift, double density_mod, double *phi[NFIELDS], short int xy_in[NX*NY])
     
-    init_laminar_flow_earth(0.04, phi, xy_in);
+//     init_laminar_flow_earth(0.01, phi, xy_in);
+//     
+//     add_random_onefield_smoothed(0, 0.0, 0.03, phi, xy_in, wsphere);
+//     add_random_onefield_smoothed(1, 0.0, 0.07, phi, xy_in, wsphere);
+//     add_random_onefield_smoothed(2, 0.0, 0.07, phi, xy_in, wsphere);
     
-    /* high pressure systems, northern hemisphere */
-    init_vortex_state_sphere_mod(1, 0.5, 2.5*PID, PID + 0.6, 0.03, 0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, 0.6, 3.7*PID, PID + 0.7, 0.035, 0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, 0.2, 2.5*PID, PID + 1.3, 0.01, 0.01, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, 0.5, 1.0*PID, PID + 1.0, 0.03, 0.02, phi, xy_in, wsphere);
-//     init_vortex_state_sphere_mod(1, 0.1, 3.0*PID, 0.9*PI, 0.01, 0.01, phi, xy_in, wsphere);
-    /* low pressure systems, northern hemisphere */
-    init_vortex_state_sphere_mod(1, -0.5, 3.0*PID, PID + 0.6, 0.03, -0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, -0.2, 2.3*PID, PID + 1.1, 0.01, -0.02, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, -0.5, 1.2*PID, PID + 0.4, 0.01, -0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, -0.8, 0.3*PID, PID + 0.6, 0.01, -0.04, phi, xy_in, wsphere);
-    /* high pressure systems, southern hemisphere */
-    init_vortex_state_sphere_mod(1, -0.6, 2.4*PID, PID - 0.7, 0.03, 0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, -0.6, 1.6*PID, PID - 0.7, 0.02, 0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, -0.4, 0.5*PID, PID - 0.7, 0.02, 0.04, phi, xy_in, wsphere);
-   /* low pressure systems, southern hemisphere */
-    init_vortex_state_sphere_mod(1, 0.5, 3.4*PID, PID - 0.7, 0.03, -0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, 0.2, 2.0*PID, PID - 0.1, 0.04, -0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, 0.3, 1.0*PID, PID - 0.6, 0.03, -0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, 0.4, 0.1*PID, PID - 0.5, 0.03, -0.04, phi, xy_in, wsphere);
-    init_vortex_state_sphere_mod(1, 0.2, 0.1*PID, 0.1, 0.1, -0.01, phi, xy_in, wsphere);
 
 //     init_pressure_gradient_flow(flow_speed_schedule(0), 1.0 + PRESSURE_GRADIENT, 1.0 - PRESSURE_GRADIENT, phi, xy_in, bc_field);
 
@@ -1948,13 +2114,12 @@ void animation()
 //     init_vortex_state_sphere(0, amp, phishift, PID + thetashift, 0.15, 0.3, phi, xy_in, wsphere);
     
     
-//     init_gaussian_wave(-1.0, 0.0, 0.005, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
+//     init_gaussian_wave(0.7, 0.0, 0.0225, 0.1, SWATER_MIN_HEIGHT, phi, xy_in);
+//     init_gaussian_wave_sphere(4.0, PID, 0.05, 0.05, SWATER_MIN_HEIGHT, phi, xy_in, wsphere);
 
-//     init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
+    //     init_linear_wave(-1.0, 0.01, 5.0e-8, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
     
-//     init_linear_blob(0.0, -0.75, 0.025, 0.0, 0.001, 0.25, 0.5, SWATER_MIN_HEIGHT, phi, xy_in);
+    init_swater_laminar_flow_sphere(0, -0.75, 0.0075, 6, 1, 0.0025, SWATER_MIN_HEIGHT, phi, xy_in, wsphere);
-        
-//     init_elliptical_vibration(0.0, 0.0, 0.0075, LAMBDA, 1.0, 0.0003, 0.1, 1.0, SWATER_MIN_HEIGHT, phi, xy_in);
         
 //     add_gaussian_wave(-1.6, -0.5, 0.015, 0.25, SWATER_MIN_HEIGHT, phi, xy_in);
     
@@ -2062,7 +2227,9 @@ void animation()
         
         if (ANTISYMMETRIZE_WAVE_FCT) antisymmetrize_wave_function(phi, xy_in);
         
-        for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][NX*NY/2]);
+        for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][NY/2]);
+//         for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][NX*NY/2]);
+//         for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][0]);
         printf("\n");
         
         if (ADD_NOISE == 1) 

schrodinger.c

@@ -66,6 +66,7 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 300        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.25   /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 0.1	    /* parameter controlling the dimensions of domain */
@@ -141,6 +142,7 @@
 #define SCALE 1          /* 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 PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
 #define E_SCALE 150.0     /* scaling factor for energy representation */
 #define LOG_SCALE 1.0    /* scaling factor for energy log representation */
 #define LOG_SHIFT 0.0     /* shift of colors on log scale */
@@ -199,6 +201,7 @@
 #define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
 #define WALL_WIDTH 0.1      /* width of wall separating lenses */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 #define INITIAL_TIME 50      /* time after which to start saving frames */
 #define OSCIL_YMAX 0.35      /* defines oscillation range */
 #define MESSAGE_LDASH 14         /* length of dash for Morse code message */

sub_rde.c

@@ -125,6 +125,64 @@ void init_random_smoothed(double mean, double amplitude, double *phi[NFIELDS], s
     printf("Fields initialised\n");
 }
 
+void add_random_onefield_smoothed(int field, 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) */
+{
+    int i, j, k, in;
+    double xy[2], dist2, module, phase, scale2, random, amp;    
+
+    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 (k=0; k<NFIELDS; k++)
+    k = field;
+    {
+        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=DSMOOTH; j<NY-DSMOOTH; j++)
+            {
+                if (xy_in[i*NY+j]) 
+                    phi[k][i*NY+j] += mean + wsphere[i*NY+j].stheta*amplitude*(2.0*(double)rand()/RAND_MAX - 1.0);
+                else phi[k][i*NY+j] = 0.0;
+            }
+        }
+        random = amplitude*(2.0*(double)rand()/RAND_MAX - 1.0);
+        for (j=0; j<DSMOOTH; j++)
+        {
+            for (i=0; i<NX; i++) 
+            {
+                if (xy_in[i*NY+j]) 
+                    phi[k][i*NY+j] += mean + random*wsphere[i*NY+j].stheta;
+                else phi[k][i*NY+j] = 0.0;
+            }
+        }
+        for (j=NY-DSMOOTH; j<NY; j++)
+        {            
+            for (i=0; i<NX; i++) 
+            {
+                if (xy_in[i*NY+j]) 
+                    phi[k][i*NY+j] += mean + random*wsphere[i*NY+j].stheta;
+                else phi[k][i*NY+j] = 0.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) */
@@ -1173,6 +1231,68 @@ double init_gaussian_wave(double x, double y, double amp, double radius, double
         }
 }
 
+double init_gaussian_wave_sphere(double phi0, double theta0, double amp, double radius, double min, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
+/* initialise gaussian wave height for shallow water equation */
+{
+    int i, j;
+    double xy[2], var, a, pscal, dist2, module;
+    
+    var = radius*radius; 
+    a = amp/(sqrt(DPI)*radius);
+    
+    printf("[init_gaussian_wave_sphere]: a = %.3lg, var = %.3lg\n", a, var); 
+    
+    for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*sin(theta0);
+            pscal += wsphere[i*NY+j].ctheta*cos(theta0);
+            dist2 = acos(pscal);
+            dist2 *= dist2;
+            module = amp*exp(-dist2/var);
+            
+//             printf("[init_gaussian_wave_sphere]: dist2 = %.3lg, module = %.3lg\n", dist2, module);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] = min + a*module*wsphere[i*NY+j].stheta;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+            else
+            {
+                phi[0][i*NY+j] = min;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+            
+//             printf("[init_gaussian_wave_sphere]:");
+//             for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][0]);
+//             printf("\n");
+        }
+        for (j=0; j<DPOLE; j++)
+        {
+            phi[0][i*NY+j] = min;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+        for (j=NY-DPOLE; j<NY; j++)
+        {
+            phi[0][i*NY+j] = min;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+    }
+        
+    printf("[init_gaussian_wave_sphere]:");
+    for (j=0; j<NFIELDS; j++) printf("field[%i] = %.3lg\t", j, phi[j][0]);
+    printf("\n"); 
+}
+
 
 double init_linear_wave(double x, double vx, double amp, double radius, double min, double *phi[NFIELDS], short int xy_in[NX*NY])
 /* initialise gaussian wave height for shallow water equation */
@@ -1237,6 +1357,83 @@ double init_linear_blob(double x, double y, double vx, double vy, double amp, do
         }
 }
 
+double init_linear_blob_sphere(int add, double phi0, double theta0, double vx, double vy, double amp, double radiusx, double radiusy, double min, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
+/* initialise gaussian wave height for shallow water equation */
+{
+    int i, j;
+    double xy[2], varx, vary, a, pscal, dist2, module;
+    
+    varx = radiusx*radiusx; 
+    vary = radiusy*radiusy; 
+    a = amp/(sqrt(DPI*radiusx*radiusy));
+    
+    printf("[init_gaussian_wave_sphere]: a = %.3lg, var = %.3lg\n", a, varx); 
+    
+    if (add) for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*sin(theta0);
+            pscal += wsphere[i*NY+j].ctheta*cos(theta0);
+            dist2 = acos(pscal);
+            dist2 *= dist2;
+            module = amp*exp(-dist2/varx);
+            
+//             printf("[init_gaussian_wave_sphere]: dist2 = %.3lg, module = %.3lg\n", dist2, module);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] += a*module*wsphere[i*NY+j].stheta;
+                phi[1][i*NY+j] += vx*module*wsphere[i*NY+j].stheta;
+                phi[2][i*NY+j] += vy*module*wsphere[i*NY+j].stheta;
+            }
+        }
+    }
+    else for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            pscal = cos(wsphere[i*NY+j].phi - phi0)*wsphere[i*NY+j].stheta*sin(theta0);
+            pscal += wsphere[i*NY+j].ctheta*cos(theta0);
+            dist2 = acos(pscal);
+            dist2 *= dist2;
+            module = amp*exp(-dist2/varx);
+            
+//             printf("[init_gaussian_wave_sphere]: dist2 = %.3lg, module = %.3lg\n", dist2, module);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] = min + a*module*wsphere[i*NY+j].stheta;
+                phi[1][i*NY+j] = vx*module*wsphere[i*NY+j].stheta;
+                phi[2][i*NY+j] = vy*module*wsphere[i*NY+j].stheta;
+            }
+            else
+            {
+                phi[0][i*NY+j] = min;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+        }
+        for (j=0; j<DPOLE; j++)
+        {
+            phi[0][i*NY+j] = min;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+        for (j=NY-DPOLE; j<NY; j++)
+        {
+            phi[0][i*NY+j] = min;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+    }
+}
 
 double init_circular_vibration(double x, double y, double v, double radius, double amp, double var, double omega, double min, double *phi[NFIELDS], short int xy_in[NX*NY])
 /* initialise gaussian wave height for shallow water equation */
@@ -1350,6 +1547,76 @@ double add_gaussian_wave(double x, double y, double amp, double radius, double m
         }
 }
 
+double init_swater_laminar_flow_sphere(int add, double vx, double amp, int nx, int ny, double deltavy, double min, double *phi[NFIELDS], short int xy_in[NX*NY], t_wave_sphere wsphere[NX*NY])
+/* initialise gaussian wave height for shallow water equation */
+{
+    int i, j;
+    double xy[2], ct, snt, module;
+    
+//     printf("[init_gaussian_wave_sphere]: a = %.3lg, var = %.3lg\n", a, varx); 
+    
+    if (add) for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            snt = 0.5*sin((double)ny*wsphere[i*NY+j].theta);
+            ct = wsphere[i*NY+j].ctheta;
+            module = snt*exp(-ct*ct);
+//             st = wsphere[i*NY+j].stheta;    
+//             module = st*ct*exp(-ct*ct);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] += 0.0;
+                phi[1][i*NY+j] += vx*amp*module;
+                phi[2][i*NY+j] += deltavy*module*cos((double)nx*wsphere[i*NY+j].phi);
+            }
+        }
+    }
+    else for (i=0; i<NX; i++)
+    {
+        for (j=DPOLE; j<NY-DPOLE; j++)
+        {
+            ij_to_xy(i, j, xy);
+            xy_in[i*NY+j] = xy_in_billiard(xy[0],xy[1]);
+            
+            snt = 0.5*sin((double)ny*wsphere[i*NY+j].theta);
+            ct = wsphere[i*NY+j].ctheta;
+            module = snt*exp(-ct*ct);
+//             st = wsphere[i*NY+j].stheta;
+//             module = st*ct*exp(-ct*ct);
+            
+            if (xy_in[i*NY+j])
+            {
+                phi[0][i*NY+j] = min;
+                phi[1][i*NY+j] = vx*amp*module;
+                phi[2][i*NY+j] = deltavy*module*cos((double)nx*wsphere[i*NY+j].phi);
+            }
+            else
+            {
+                phi[0][i*NY+j] = min;
+                phi[1][i*NY+j] = 0.0;
+                phi[2][i*NY+j] = 0.0;
+            }
+        }
+        for (j=0; j<DPOLE; j++)
+        {
+            phi[0][i*NY+j] = min;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+        for (j=NY-DPOLE; j<NY; j++)
+        {
+            phi[0][i*NY+j] = min;
+            phi[1][i*NY+j] = 0.0;
+            phi[2][i*NY+j] = 0.0;
+        }
+    }
+}
+
 double distance_to_segment(double x, double y, double x1, double y1, double x2, double y2)
 /* distance of (x,y) to segment from (x1,y1) to (x2,y2) */
 {
@@ -1637,7 +1904,38 @@ void initialize_bcfield(double bc_field[NX*NY], double bc_field2[NX*NY], t_recta
         }
         case (D_SPHERE_EARTH):
         {
-            /* do nothing? */
+            for (i=0; i<NX; i++)
+                for (j=0; j<NY; j++)
+                {
+                    /* set dummy values */
+                    height = wsphere[i*NY+j+1].altitude;
+                    if (height == 0.0) bc_field[i*NY+j] = 1.0;
+                    else 
+                    {   
+                        f = tanh(BC_STIFFNESS*height);
+                        bc_field[i*NY+j] = 1.0 - f;
+                    }
+                 }
+            break;
+        }
+        case (D_SPHERE_MARS):
+        {
+            for (i=0; i<NX; i++)
+                for (j=0; j<NY; j++)
+                {
+                    /* set dummy values */
+                    height = wsphere[i*NY+j+1].altitude;
+                    if (height == 0.0) bc_field[i*NY+j] = 1.0;
+                    else 
+                    {   
+                        f = tanh(BC_STIFFNESS*height);
+                        bc_field[i*NY+j] = 1.0 - f;
+                    }
+                 }
+            break;
+        }
+        case (D_SPHERE_VENUS):
+        {
             for (i=0; i<NX; i++)
                 for (j=0; j<NY; j++)
                 {
@@ -2802,7 +3100,7 @@ void compute_gradient_polar(t_rde rde[NX*NY], double factor)
         for (j=0; j<NY; j++)
         {
             rde[i*NY+j].field_norm = factor*module2(rde[i*NY+j].nablax, rde[i*NY+j].nablay);
-            angle = argument(rde[i*NY+j].nablax, rde[i*NY+j].nablay) + PI*COLOR_PHASE_SHIFT;
+            angle = argument(rde[i*NY+j].nablax, rde[i*NY+j].nablay) + PI*PHASE_SHIFT;
             if (angle < 0.0) angle += DPI;
             if (angle >= DPI) angle -= DPI;
             rde[i*NY+j].field_arg = angle;
@@ -2832,7 +3130,7 @@ void compute_field_argument(double *phi[NFIELDS], t_rde rde[NX*NY])
     for (i=0; i<NX; i++)
         for (j=0; j<NY; j++)
         {
-            arg = argument(phi[0][i*NY+j], phi[1][i*NY+j]) + COLOR_PHASE_SHIFT*PI;
+            arg = argument(phi[0][i*NY+j], phi[1][i*NY+j]) + PHASE_SHIFT*PI;
             while (arg < 0.0) arg += DPI;
             while (arg >= DPI) arg -= DPI;
             rde[i*NY+j].field_arg = arg;
@@ -2953,6 +3251,7 @@ void adjust_height(double *phi[NFIELDS], t_rde rde[NX*NY])
         {
             value = VSCALE_WATER_HEIGHT*(phi[0][i*NY+j] + ADD_HEIGHT_CONSTANT);
             rde[i*NY+j].height = value;
+//             printf("[adjust_height] value = %.3lg\n", value);
         }
 }
 
@@ -2961,6 +3260,14 @@ void compute_direction(double *phi[NFIELDS], t_rde rde[NX*NY])
 {
     int i, j; 
     double value;
+//     static double phaseshift;
+//     static int first = 1;
+//     
+//     if (first)
+//     {
+//         phaseshift = PHASE_SHIFT*PID/90.0;
+//         first = 0;
+//     }
     
     #pragma omp parallel for private(i,j,value)
     for (i=0; i<NX; i++)
@@ -3934,7 +4241,7 @@ void compute_probabilities(t_rde rde[NX*NY], short int xy_in[NX*NY], double prob
 //     for (i=0; i<NX; i++)
 //         for (j=0; j<NY; j++)
 //         {
-//             angle = argument(phi_x[i*NY+j],phi_y[i*NY+j]) + PI*COLOR_PHASE_SHIFT;
+//             angle = argument(phi_x[i*NY+j],phi_y[i*NY+j]) + PI*PHASE_SHIFT;
 //             if (angle < 0.0) angle += DPI;
 //             if (angle >= DPI) angle -= DPI;
 //             phi_arg[i*NY+j] = 360.0*angle/DPI;
@@ -4852,7 +5159,7 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3]
         case (Z_POLAR): 
         {
 //             hsl_to_rgb_palette(360.0*value/DPI, 0.9, 0.5, rgb, palette);
-            value += COLOR_PHASE_SHIFT*PI;
+            value += PHASE_SHIFT*PI;
             if (value > DPI) value -= DPI;
             color_scheme_palette(C_ONEDIM_LINEAR, palette, value/DPI, 1.0, 1, rgb);
             break;
@@ -5011,7 +5318,7 @@ void compute_field_color_rde(double value, int cplot, int palette, double rgb[3]
         case (Z_SWATER_DIRECTION_SPEED): 
         {
             value *= 360.0/DPI;
-            value += 180.0*COLOR_PHASE_SHIFT;
+            value += 180.0*PHASE_SHIFT;
             if (value > 360.0) value -= 360.0; 
             hsl_to_rgb_palette(value, 0.9, 0.5, rgb, palette);
             break;
@@ -5163,6 +5470,7 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot,
 {
     int i, j;
     
+//     printf("[compute_rde_fields] zplot = %i\n", zplot); 
     switch (RDE_EQUATION) {
         case (E_RPS):
         {
@@ -5233,11 +5541,13 @@ void compute_rde_fields(double *phi[NFIELDS], short int xy_in[NX*NY], int zplot,
         }
         case (E_SHALLOW_WATER):
         {
+//             printf("[compute_rde_fields]\n"); 
             if (zplot == Z_SWATER_HEIGHT) adjust_height(phi, rde);
             if ((zplot == Z_SWATER_SPEED)||(cplot == Z_SWATER_SPEED)||(zplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED))
                 compute_speed(phi, rde, wsphere);
             if ((zplot == Z_SWATER_DIRECTION_SPEED)||(cplot == Z_SWATER_DIRECTION_SPEED))
                 compute_direction(phi, rde);
+            break;
         }
         default : break;
     }
@@ -5605,6 +5915,7 @@ void init_cfield_rde(double *phi[NFIELDS], short int xy_in[NX*NY], int cplot, t_
         {
             #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] = &phi[0][i*NY+j];
+//             for (i=0; i<NX; i++) for (j=0; j<NY; j++) rde[i*NY+j].p_cfield[movie] = &rde[i*NY+j].height;
             break;
         }
         case (Z_SWATER_SPEED):
@@ -5727,7 +6038,8 @@ void draw_wave_sphere_rde_ij(int i, int iplus, int j, int jplus, int jcolor, int
     if (NON_DIRICHLET_BC) 
         draw_bc = (xy_in[i*NY+j])&&(xy_in[iplus*NY+j])&&(xy_in[i*NY+jplus])&&(xy_in[iplus*NY+jplus]);
     
-    draw = wsphere[i*NY+j].indomain;
+//     draw = wsphere[i*NY+j].indomain;
+    draw = 1;
     for (p=0; p<HRES; p++)
         for (q=0; q<HRES; q++)
             draw_hr[p*HRES+q] = wsphere_hr[(HRES*i+p)*HRES*NY+HRES*j+q].indomain;
@@ -5762,7 +6074,9 @@ void draw_wave_sphere_rde_ij(int i, int iplus, int j, int jplus, int jcolor, int
         {
             glBegin(GL_TRIANGLE_FAN);
             glColor3f(rgb[0], rgb[1], rgb[2]);
+//             printf("[draw_wave_sphere_rde_ij] draw = %i, zfield = %.3lg\n", draw, *rde[i*NY+j].p_zfield[movie]);
             drawij = ij_to_sphere(i, j, *rde[i*NY+j].p_zfield[movie], wsphere, xyz, draw);
+            
             if (drawij)
                 draw_vertex_sphere(xyz);
             if (ij_to_sphere(iplus, j, *rde[iplus*NY+j].p_zfield[movie], wsphere, xyz, draw))
@@ -5871,7 +6185,7 @@ void draw_arrow(double z, int fade, double fade_value)
 
 void draw_wave_sphere_3d_rde(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)
 {
-    int i, j, imax, imin, jmax, imid, jmid;
+    int i, j, imax, imin, jmin, jmax, imid, jmid;
     double observer_angle, angle2, observer_latitude, xyz[3];
     
     blank();
@@ -5897,7 +6211,8 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX
     if (imin >= NX-1) imin = NX-2;
     if (imax >= NX-1) imax = NX-2;
     
-    jmax = NY-3;
+    jmax = NY-POLE_NODRAW;
+    jmin = POLE_NODRAW;
 //     jmax = NY - 50;
     jmid = (int)((double)NY*(observer_latitude + PID)/PI);
     imid = (imin + imax)/2;
@@ -5910,12 +6225,12 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX
         {
             for (i=imax; i>imid; i--)
             {
-                for (j=0; j<=jmax; j++)
+                for (j=jmin; j<=jmax; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
             for (i=imid; i>imin; i--)
             {
-                for (j=0; j<=jmid; j++)
+                for (j=jmin; j<=jmid; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
                 for (j=jmax; j>=jmid; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
@@ -5923,16 +6238,16 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX
             
             for (i=imax+1; i<NX-1; i++)
             {
-                for (j=0; j<=jmax; j++)
+                for (j=jmin; j<=jmax; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
         
-            for (j=0; j<=jmax; j++)
+            for (j=jmin; j<=jmax; j++)
                 draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
         
             for (i=0; i<=imin; i++)
             {
-                for (j=0; j<=jmid; j++)
+                for (j=jmin; j<=jmid; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
                 for (j=jmax; j>=jmid; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
@@ -5950,13 +6265,13 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX
         {
             for (i=imax; i<imid; i++)
             {
-                for (j=0; j<=jmax; j++)
+                for (j=jmin; j<=jmax; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
         
             for (i=imid; i<imin; i++)
             {
-                for (j=0; j<=jmid; j++)
+                for (j=jmin; j<=jmid; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
                 for (j=jmax; j>=jmid; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
@@ -5964,16 +6279,16 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX
         
             for (i=imax-1; i>=0; i--)
             {
-                for (j=0; j<=jmax; j++)
+                for (j=jmin; j<=jmax; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
             
-            for (j=0; j<=jmax; j++)
+            for (j=jmin; j<=jmax; j++)
                 draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
         
             for (i=NX-2; i>=imin; i--)
             {
-                for (j=0; j<=jmid; j++)
+                for (j=jmin; j<=jmid; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
                 for (j=jmax; j>=jmid; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
@@ -6006,7 +6321,7 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX
         {
             for (i=imax; i>imid; i--)
             {
-                for (j=jmax; j>=0; j--)
+                for (j=jmax; j>=jmin; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
             
@@ -6014,59 +6329,102 @@ void draw_wave_sphere_3d_rde(int movie, double *phi[NFIELDS], short int xy_in[NX
             {
                 for (j=jmax; j>=jmid; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
-                for (j=0; j<=jmid; j++)
+                for (j=jmin; j<=jmid; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
             
             for (i=imax+1; i<NX-1; i++)
             {
-                for (j=jmax; j>=0; j--)
+                for (j=jmax; j>=jmin; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
             
-            for (j=jmax; j>=0; j--)
+            for (j=jmax; j>=jmin; j--)
                 draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
         
             for (i=0; i<=imin; i++)
             {
                 for (j=jmax; j>=jmid; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
-                for (j=0; j<=jmid; j++)
+                for (j=jmin; j<=jmid; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
+            
+            /* TEST */
+            for (j=jmin; j<=jmid; j++)
+                for (i=-1; i<2; i++) if ((imin+i >= 0)&&(imin+i+1 < NX))
+                    draw_wave_sphere_rde_ij(imin+i, imin+i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+            for (j=jmax; j>=jmid; j--)
+                for (i=-1; i<2; i++) if ((imin+i >= 0)&&(imin+i+1 < NX))
+                    draw_wave_sphere_rde_ij(imin+i, imin+i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+            
+//             for (i=0; i<2; i++) if (imin >= i)
+//             {
+//                 for (j=2*jmax/3; j>=jmid; j--)
+//                     draw_wave_sphere_rde_ij(imin-i, imin-i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+//                 for (j=jmin; j<=jmid; j++)
+//                     draw_wave_sphere_rde_ij(imin-i, imin-i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+//             }
+
         }
         else
         {
             for (i=imax; i<imid; i++)
             {
-                for (j=jmax; j>=0; j--)
+                for (j=jmax; j>=jmin; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
             
-            for (i=imid; i<imin; i++)
+            for (i=imid; i<imin-1; i++)
             {
                 for (j=jmax; j>=jmid; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
-                for (j=0; j<=jmid; j++)
+                for (j=jmin; j<=jmid; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
             
             for (i=imax-1; i>=0; i--)
             {
-                for (j=jmax; j>=0; j--)
+                for (j=jmax; j>=jmin; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
             
-            for (j=jmax; j>=0; j--)
+            for (j=jmax; j>=jmin; j--)
                 draw_wave_sphere_rde_ij(NX-1, 0, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
         
-            for (i=NX-2; i>=imin-1; i--)
+            for (i=NX-2; i>=imin; i--)
             {
                 for (j=jmax; j>=jmid; j--)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
-                for (j=0; j<=jmid; j++)
+                for (j=jmin; j<=jmid; j++)
                     draw_wave_sphere_rde_ij(i, i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
             }
+            
+            for (i=0; i<2; i++) if (imin >= i)
+            {
+                for (j=2*jmax/3; j>=jmid; j--)
+                    draw_wave_sphere_rde_ij(imin-i, imin-i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+                for (j=jmin; j<=jmid; j++)
+                    draw_wave_sphere_rde_ij(imin-i, imin-i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+            }
+            
+            /* TEST */
+            for (j=jmin; j<=jmid; j++)
+                for (i=-1; i<2; i++) if ((imin+i >= 0)&&(imin+i+1 < NX))
+                    draw_wave_sphere_rde_ij(imin+i, imin+i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+            for (j=jmax; j>=jmid; j--)
+                for (i=-1; i<2; i++) if ((imin+i >= 0)&&(imin+i+1 < NX))
+                    draw_wave_sphere_rde_ij(imin+i, imin+i+1, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+            
+//             if (imin >= 1)
+//             {
+//                 /* experimental */
+//                 for (j=2*jmax/3; j>=jmid; j--)
+//                     draw_wave_sphere_rde_ij(imin-1, imin, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+//                 for (j=0; j<=jmid; j++)
+//                     draw_wave_sphere_rde_ij(imin-1, imin, j, j+1, j+1, movie, phi, xy_in, rde, wsphere, wsphere_hr, zplot, cplot, palette, fade, fade_value);
+//             }
+
         }
     
         /* South pole */
@@ -6329,9 +6687,9 @@ void draw_wave_rde(int movie, double *phi[NFIELDS], short int xy_in[NX*NY], t_rd
 
     if (refresh)
     {
-//         printf("Computing fields\n");
+        printf("Computing fields\n");
         compute_rde_fields(phi, xy_in, zplot, cplot, rde, wsphere);      
-//         printf("Computed fields\n");
+        printf("Computed fields\n");
         if ((PLOT_3D)&&(SHADE_3D)) 
         {
             if (SPHERE) 

sub_wave.c

@@ -685,12 +685,74 @@ double michelson_schedule(int i)
 }   
 
 
+int generate_poisson_discs(t_circle circles[NMAXCIRCLES], double xmin, double xmax, double ymin, double ymax, double dpoisson)
+/* generate a Poisson disc sample in a given rectangle */
+{
+    int i, j, k, n_p_active, ncandidates=5000, naccepted; 
+    double r, phi, x, y;
+    short int active_poisson[NMAXCIRCLES], far;
+    
+    printf("Generating Poisson disc sample\n");
+    /* generate first circle */
+    circles[0].xc = (xmax - xmin)*(double)rand()/RAND_MAX + xmin;
+    circles[0].yc = (ymax - ymin)*(double)rand()/RAND_MAX + ymin;
+    active_poisson[0] = 1;
+    n_p_active = 1;
+    ncircles = 1;
+            
+    while ((n_p_active > 0)&&(ncircles < NMAXCIRCLES))
+    {
+        /* randomly select an active circle */
+        i = rand()%(ncircles);
+        while (!active_poisson[i]) i = rand()%(ncircles);                 
+        /* 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 = circles[i].xc + r*cos(phi);
+            y = circles[i].yc + r*sin(phi);
+            far = 1;
+            for (k=0; k<ncircles; k++) if ((k!=i))
+            {
+                /* new circle is far away from circle k */
+                far = far*((x - circles[k].xc)*(x - circles[k].xc) + (y - circles[k].yc)*(y - circles[k].yc) >= dpoisson*dpoisson);
+                /* new circle is in domain */
+                far = far*(x < xmax)*(x > xmin)*(y < ymax)*(y > ymin);
+            }
+            if (far)    /* accept new circle */
+            {
+                printf("New circle at (%.3f,%.3f) accepted\n", x, y);
+                circles[ncircles].xc = x;
+                circles[ncircles].yc = y;
+                circles[ncircles].radius = MU;
+                circles[ncircles].active = 1;
+                active_poisson[ncircles] = 1;
+                ncircles++;
+                n_p_active++;
+                naccepted++;
+            }
+        }
+        if (naccepted == 0)    /* inactivate circle i */ 
+        {
+            active_poisson[i] = 0;
+            n_p_active--;
+        }
+        printf("%i active circles\n", n_p_active);
+    }
+            
+    printf("Generated %i circles\n", ncircles);
+    return(ncircles);
+}
+
+
 int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern)
 /* initialise the arrays circlex, circley, circlerad and circleactive */
 /* for billiard shape D_CIRCLES */
 {
     int i, j, k, n, ncirc0, n_p_active, ncandidates=5000, naccepted; 
-    double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = 3.25*MU, xx[4], yy[4], dr, dphi;
+    double dx, dy, p, phi, r, r0, ra[5], sa[5], height, x, y = 0.0, gamma, dpoisson = PDISC_FACTOR*MU, xx[4], yy[4], dr, dphi;
     short int active_poisson[NMAXCIRCLES], far;
     
     switch (circle_pattern) {
@@ -737,7 +799,8 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
                 {
                     n = NGRIDY*i + j;
                     circles[n].xc = ((double)(i-NGRIDX/2) + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy;
-                    circles[n].yc = YMIN + 0.5 + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy;
+                    circles[n].yc = YMIN + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy;
+//                     circles[n].yc = YMIN + 0.5 + ((double)j + 0.5 + 0.5*((double)rand()/RAND_MAX - 0.5))*dy;
                     circles[n].radius = MU*sqrt(1.0 + 0.8*((double)rand()/RAND_MAX - 0.5));
                     circles[n].active = 1;
                 }
@@ -955,7 +1018,7 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
             r0 = 2.0*MU;
             r = r0 + MU;
             
-            for (i=0; i<1000; i++) 
+            for (i=0; i<2000; i++) 
             {
                 x = r*cos(phi);
                 y = r*sin(phi);
@@ -1003,13 +1066,13 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
         {
             ncircles = NGRIDX*NGRIDY;
             dphi = DPI/((double)NGRIDX);
-            dr = 0.5*LAMBDA/(double)NGRIDY;
+            dr = (1.0 - RADIUS_FACTOR)*LAMBDA/(double)NGRIDY;
             for (i = 0; i < NGRIDX; i++)
                 for (j = 0; j < NGRIDY; j++)
                 {
                     n = NGRIDY*i + j;
                     phi = (double)i*dphi;
-                    r = 0.5*LAMBDA + (double)j*dr;
+                    r = RADIUS_FACTOR*LAMBDA + (double)j*dr;
                     circles[n].xc = r*cos(phi);
                     circles[n].yc = r*sin(phi);
                     circles[n].radius = MU;
@@ -1023,14 +1086,14 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
         {
             ncircles = NGRIDX*NGRIDY;
             dphi = DPI/((double)NGRIDX);
-            dr = 0.5*LAMBDA/(double)NGRIDY;
+            dr = (1.0 - RADIUS_FACTOR)*LAMBDA/(double)NGRIDY;
             for (i = 0; i < NGRIDX; i++)
                 for (j = 0; j < NGRIDY; j++)
                 {
                     n = NGRIDY*i + j;
                     phi = (double)i*dphi;
                     phi += 0.5*(double)j*dphi;
-                    r = 0.5*LAMBDA + (double)j*dr;
+                    r = RADIUS_FACTOR*LAMBDA + (double)j*dr;
                     circles[n].xc = r*cos(phi);
                     circles[n].yc = r*sin(phi);
                     circles[n].radius = MU;
@@ -1048,10 +1111,11 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
             
             gamma = (sqrt(5.0) - 1.0)*PI;    /* golden mean times 2Pi */
             phi = 0.0;
-            r0 = 0.5*LAMBDA;
+//             r0 = 0.5*LAMBDA;
+            r0 = RADIUS_FACTOR*LAMBDA;
             r = r0 + MU;
             
-            for (i=0; i<1000; i++) 
+            for (i=0; i<(int)(1000.0/RADIUS_FACTOR); i++) 
             {
                 x = r*cos(phi);
                 y = r*sin(phi);
@@ -1075,6 +1139,40 @@ int init_circle_config_pattern(t_circle circles[NMAXCIRCLES], int circle_pattern
             }
             break;
         }
+        case (C_RINGS_POISSONDISC):
+        {
+            ncircles = generate_poisson_discs(circles, YMIN, YMAX, YMIN, YMAX, PDISC_FACTOR*MU);
+            for (i=0; i<ncircles; i++)
+            {
+                circles[i].radius = MU;
+                /* inactivate circles outside the domain */
+                r = module2(circles[i].xc, circles[i].yc);
+                if ((r < LAMBDA)&&(r > RADIUS_FACTOR*LAMBDA)) circles[i].active = 1;
+                else circles[i].active = 0;
+            }
+            break;
+        }
+        case (C_RINGS_LOGSPIRAL):
+        {
+            ncircles = NGRIDX*NGRIDY;
+            dr = (1.0 - RADIUS_FACTOR)*LAMBDA/(double)NGRIDY;
+            dphi = DPI/((double)NGRIDX);
+            for (i = 0; i < NGRIDX; i++)
+                for (j = 0; j < NGRIDY; j++)
+                {
+                    n = NGRIDY*i + j;
+                    r = RADIUS_FACTOR*LAMBDA + (double)j*dr;
+                    phi = (double)i*dphi;
+                    phi += log(r/(RADIUS_FACTOR*LAMBDA));
+                    circles[n].xc = r*cos(phi);
+                    circles[n].yc = r*sin(phi);
+                    circles[n].radius = MU;
+                    /* activate only circles that intersect the domain */
+                    if ((circles[n].yc < YMAX + MU)&&(circles[n].yc > YMIN - MU)) circles[n].active = 1;
+                    else circles[n].active = 0;
+                }
+            break;
+        }
         case (C_HEX_BOTTOM):
         {
             ncircles = NGRIDX*(NGRIDY+1);
@@ -2545,6 +2643,60 @@ int rc_hyp(double x, double y)
     }
 }
 
+
+int init_xyin_from_image(short int * xy_in[NX])
+/* initialize table xy_in from an image file */
+{
+    FILE *image_file;
+    int nx, ny, maxrgb, i, j, k, ii, jj, nmaxpixels = 2000000, rgbtot, scan, rgbval;
+    int *rgb_values;
+    double scalex, scaley;
+
+    image_file = fopen("PHOTONSband.ppm", "r");
+    scan = fscanf(image_file,"%i %i\n", &nx, &ny);
+    scan = fscanf(image_file,"%i\n", &maxrgb);
+    
+    rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int));
+    
+    scalex = (double)nx/(double)NX;
+    scaley = (double)ny/(double)NY;
+    
+    if (nx*ny > nmaxpixels)
+    {
+        printf("Image too large, increase nmaxpixels in init_xyin_from_image()\n");
+        exit(0);
+    }
+    
+    /* read rgb values */
+    for (j=0; j<ny; j++)
+        for (i=0; i<nx; i++)
+            for (k=0; k<3; k++)
+                {
+                    scan = fscanf(image_file,"%i\n", &rgbval);
+                    rgb_values[3*(j*nx+i)+k] = rgbval;
+                }
+                
+    for (i=0; i<NX; i++)
+        for (j=0; j<NY; j++)
+        {
+            ii = nx/2 + (int)((double)(i-NX/2)*scalex);
+            jj = ny/2 - (int)((double)(j-NY/2)*scalex);
+            if (ii > nx-1) ii = nx-1;
+            if (ii < 0) ii = 0;
+            if (jj > ny-1) jj = ny-1;
+            if (jj < 0) jj = 0;
+            k = 3*(jj*nx+ii);
+            rgbtot = rgb_values[k] + rgb_values[k+1] + rgb_values[k+2];
+            if (rgbtot > 3*maxrgb/2) xy_in[i][j] = 1;
+            else xy_in[i][j] = 0;
+        }   
+    
+    fclose(image_file);
+    free(rgb_values);
+    return(1);
+}
+
+
 int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_circle *circles)
 /* returns 1 if (x,y) represents a point in the billiard */
 {
@@ -3547,6 +3699,52 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
             /* use IOR_GRADIENT_INDEX_LENS for index of refraction control */
             return(1);
         }
+        case (D_IMAGE):
+        {
+            /* Not needed, uses option XYIN_INITIALISED */
+            return(1);
+        }
+        case (D_MAGNETRON):
+        {
+            r = module2(x,y);
+            if (r > LAMBDA) return(1);
+            if (r < 0.33*LAMBDA) return(1);
+            if ((vabs(y) < WALL_WIDTH)&&(x > -0.66*LAMBDA)) return(1);
+            if ((vabs(x) < WALL_WIDTH)&&(vabs(y) < 0.66*LAMBDA)) return(1);
+            if ((vabs(x-y) < WALL_WIDTH*sqrt(2.0))&&(r < 0.66*LAMBDA)) return(1);
+            if ((vabs(x+y) < WALL_WIDTH*sqrt(2.0))&&(r < 0.66*LAMBDA)) return(1);
+            for (k=0; k<8; k++)
+            {
+                x1 = 0.66*cos((double)k*0.5*PID);
+                y1 = 0.66*sin((double)k*0.5*PID);
+                r = module2(x - x1, y - y1);
+                if (r < MU) return(1);
+            }
+            
+            
+            return(0);
+        }
+        case (D_MAGNETRON_CATHODE):
+        {
+            r = module2(x,y);
+            if (r < WALL_WIDTH) return(0);
+            if (r > LAMBDA) return(1);
+            if (r < 0.33*LAMBDA) return(1);
+            if ((vabs(y) < WALL_WIDTH)&&(x > -0.66*LAMBDA)) return(1);
+            if ((vabs(x) < WALL_WIDTH)&&(vabs(y) < 0.66*LAMBDA)) return(1);
+            if ((vabs(x-y) < WALL_WIDTH*sqrt(2.0))&&(r < 0.66*LAMBDA)) return(1);
+            if ((vabs(x+y) < WALL_WIDTH*sqrt(2.0))&&(r < 0.66*LAMBDA)) return(1);
+            for (k=0; k<8; k++)
+            {
+                x1 = 0.66*cos((double)k*0.5*PID);
+                y1 = 0.66*sin((double)k*0.5*PID);
+                r = module2(x - x1, y - y1);
+                if (r < MU) return(1);
+            }
+            
+            
+            return(0);
+        }
         case (D_MENGER):       
         {
             x1 = 0.5*(x+1.0);
@@ -5630,6 +5828,21 @@ void draw_billiard(int fade, double fade_value)      /* draws the billiard bound
             draw_line(WAVE_PROFILE_X, YMIN, WAVE_PROFILE_X, YMAX);
             break;
         }
+        case (D_IMAGE):
+        {
+            /* do nothing */
+            break;
+        }
+        case (D_MAGNETRON):
+        {
+            /* TODO */
+            break;
+        }
+        case (D_MAGNETRON_CATHODE):
+        {
+            /* TODO */
+            break;
+        }
         case (D_MENGER):
         {
             glLineWidth(3);
@@ -6337,7 +6550,7 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2,
 void draw_circular_color_scheme_palette_fade(double x1, double y1, double radius, int plot, double min, double max, int palette, int fade, double fade_value)
 {
     int j, k;
-    double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2], zscale = 0.9;
+    double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2], zscale = 0.85;
     
 //     glBegin(GL_TRIANGLE_FAN);
     xy_to_pos(x1, y1, xy);
@@ -6465,6 +6678,13 @@ void draw_circular_color_scheme_palette_fade(double x1, double y1, double radius
                 color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
                 break;
             }
+            case (Z_SWATER_DIRECTION_SPEED):
+            {
+                value = dy_phase*(double)(j) - 0.5*PHASE_SHIFT + 0.5;
+                if (value > 1.0) value -= 1.0;
+                color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
+                break;
+            }
             default:
             {
                 value = min + 1.0*dy*(double)(j);

sub_wave_3d_rde.c

@@ -604,6 +604,145 @@ void init_earth_map_rde(t_wave_sphere *wsphere, int res)
         }
 }
 
+
+void init_planet_map_rde(t_wave_sphere wsphere[NX*NY], int planet, int res)
+/* init file from planetary map */
+{
+    int i, j, ii, jj, k, nx, ny, maxrgb, nmaxpixels, scan, rgbval, diff, sshift, nshift, ishift, dem_number;
+    int *rgb_values;
+    double cratio, rx, ry, cy, vshift;
+    FILE *image_file;
+    
+    switch (planet){
+        case (D_SPHERE_MARS): 
+        {
+            printf("Reading Mars map\n");
+            nmaxpixels = 8388608;
+            image_file = fopen("Mars_Viking_ClrMosaic_global_925m_scaled.ppm", "r");
+            dem_number = DEM_MARS;
+            break;
+        }
+        case (D_SPHERE_MOON):
+        {
+            printf("Reading Moon map\n");
+            nmaxpixels = 2048*1024;
+            image_file = fopen("Moon_photo_map.ppm", "r");
+            dem_number = DEM_MOON;
+            break;
+        }
+        case (D_SPHERE_VENUS):
+        {
+            printf("Reading Venus map\n");
+            nmaxpixels = 1440*720;
+            image_file = fopen("Venus_map_NASA_JPL_Magellan-Venera-Pioneer.ppm", "r");
+            dem_number = DEM_VENUS;
+            break;
+        }
+        case (D_SPHERE_MERCURY):
+        {
+            printf("Reading Mercury map\n");
+            nmaxpixels = 2304*1152;
+            image_file = fopen("Mercury_color_photo.ppm", "r");
+            dem_number = DEM_MERCURY;
+            break;
+        }
+    }
+    
+    scan = fscanf(image_file,"%i %i\n", &nx, &ny);
+    scan = fscanf(image_file,"%i\n", &maxrgb);
+    printf("nx*ny = %i\n", nx*ny);
+    rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int));
+    
+    if (nx*ny > nmaxpixels)
+    {
+        printf("Image too large, increase nmaxpixels in init_planet_map()\n");
+        exit(0);
+    }
+    
+    /* shift due to min/max latitudes of image */
+    sshift = 0 + DPOLE;
+    nshift = 0 + DPOLE; 
+    
+    /* choice of zero meridian */
+    ishift = (int)(nx*ZERO_MERIDIAN/360.0);
+    
+    /* read rgb values */
+    for (j=0; j<ny; j++)
+        for (i=0; i<nx; i++)
+            for (k=0; k<3; k++)
+                {
+                    scan = fscanf(image_file,"%i\n", &rgbval);
+                    rgb_values[3*(j*nx+i)+k] = rgbval;
+                }
+                    
+    cratio = 1.0/(double)maxrgb;
+    rx = (double)nx/(double)(res*NX);
+    ry = (double)ny/(double)(res*NY - sshift - nshift);
+    
+    printf("cratio = %.3lg, rx = %.3lg, ry = %.3lg\n", cratio, rx, ry);
+    
+//     cy = rx*(double)(NY - nshift);
+    
+    /* build wave table */
+    for (i=0; i<res*NX; i++)
+        for (j=0; j<res*NY; j++)
+        {
+            ii = (int)(rx*(double)(res*NX-1 - i)) + nx/2 + ishift;
+            if (ii > nx-1) ii -= nx;
+//             jj = (int)(-ry*(double)j + cy);
+//             jj = (int)(ry*(double)(NY-nshift - j)) + sshift;
+            jj = (int)(ry*(double)(res*NY-nshift - j));
+            if (jj > ny-1) jj = ny-1;
+            if (jj < 0) jj = 0;
+            wsphere[i*res*NY+j].r = (double)rgb_values[3*(jj*nx+ii)]*cratio;
+            wsphere[i*res*NY+j].g = (double)rgb_values[3*(jj*nx+ii)+1]*cratio;
+            wsphere[i*res*NY+j].b = (double)rgb_values[3*(jj*nx+ii)+2]*cratio;
+            
+//             printf("RGB at (%i, %i) = (%.3lg, %3.lg, %.3lg)\n", i, j, wsphere[i*NY+j].r, wsphere[i*NY+j].g, wsphere[i*NY+j].b);
+            
+            /* decide which points are in the Sea */
+            wsphere[i*NY+j].indomain = 1;
+            wsphere[i*NY+j].draw_wave = 1;
+        }
+
+    
+    /* smooth colors in case of high resolution */
+    if ((res*NX > nx)||(res*NY > ny))
+        for (i=1; i<res*NX-1; i++)
+            for (j=1; j<res*NY-1; j++)
+            {
+                wsphere[i*res*NY+j].r *= 0.2; 
+                wsphere[i*res*NY+j].r += 0.2*wsphere[(i+1)*res*NY+j].r;
+                wsphere[i*res*NY+j].r += 0.2*wsphere[(i-1)*res*NY+j].r;
+                wsphere[i*res*NY+j].r += 0.2*wsphere[i*res*NY+j-1].r;
+                wsphere[i*res*NY+j].r += 0.2*wsphere[i*res*NY+j+1].r;
+
+                wsphere[i*res*NY+j].g *= 0.2; 
+                wsphere[i*res*NY+j].g += 0.2*wsphere[(i+1)*res*NY+j].g;
+                wsphere[i*res*NY+j].g += 0.2*wsphere[(i-1)*res*NY+j].g;
+                wsphere[i*res*NY+j].g += 0.2*wsphere[i*res*NY+j-1].g;
+                wsphere[i*res*NY+j].g += 0.2*wsphere[i*res*NY+j+1].g;
+
+                wsphere[i*res*NY+j].b *= 0.2; 
+                wsphere[i*res*NY+j].b += 0.2*wsphere[(i+1)*res*NY+j].b;
+                wsphere[i*res*NY+j].b += 0.2*wsphere[(i-1)*res*NY+j].b;
+                wsphere[i*res*NY+j].b += 0.2*wsphere[i*res*NY+j-1].b;
+                wsphere[i*res*NY+j].b += 0.2*wsphere[i*res*NY+j+1].b;
+            }
+
+    
+    free(rgb_values);
+    fclose(image_file);
+    
+    if (ADD_DEM) init_dem_rde(wsphere, dem_number, res);
+    
+    vshift = PLANET_SEALEVEL/(DEM_MAXHEIGHT - DEM_MAXDEPTH);
+    for (i=0; i<res*NX; i++)
+        for (j=0; j<res*NY; j++)
+            wsphere[i*res*NY+j].indomain = (wsphere[i*res*NY+j].altitude < vshift);
+}
+
+
 int ij_to_sphere(int i, int j, double r, t_wave_sphere wsphere[NX*NY], double xyz[3], int use_wave_radius)
 /* convert spherical to rectangular coordinates */
 {
@@ -2348,7 +2487,7 @@ void draw_color_scheme_palette_3d(double x1, double y1, double x2, double y2, in
 void draw_circular_color_scheme_palette_3d(double x1, double y1, double radius, int plot, double min, double max, int palette, int fade, double fade_value)
 {
     int j, k, ij_center[2], ij_right[2], ic, jc, ir;
-    double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2], zscale = 0.9;
+    double x, y, dy, dy_e, dy_phase, rgb[3], value, lum, amp, dphi, pos[2], phi, xy[2], zscale = 0.95;
     
 //     printf("Drawing color bar\n");
     
@@ -2508,7 +2647,9 @@ void draw_circular_color_scheme_palette_3d(double x1, double y1, double radius,
             }
             case (Z_SWATER_DIRECTION_SPEED):
             {
-                value = dy_phase*(double)(j);
+                value = dy_phase*(double)(j) + 0.5 - 0.5*PHASE_SHIFT;
+                if (value > 1.0) value -= 1.0;
+                if (value < 0.0) value += 1.0;
                 color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
                 break;
             }

wave_3d.c

@@ -43,7 +43,7 @@
 #include <omp.h>
 #include <time.h>
 
-#define MOVIE 0         /* set to 1 to generate movie */
+#define MOVIE 1         /* set to 1 to generate movie */
 #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 */
@@ -83,6 +83,7 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 1000        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.25   /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 1.0	    /* parameter controlling the dimensions of domain */
@@ -97,6 +98,7 @@
 #define NGRIDX 30           /* number of grid point for grid of disks */
 #define NGRIDY 18           /* number of grid point for grid of disks */
 #define WALL_WIDTH 0.1      /* width of wall separating lenses */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 
 #define X_SHOOTER -0.2
 #define Y_SHOOTER -0.6
@@ -227,6 +229,7 @@
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
 #define VSCALE_AMPLITUDE 2.0   /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define VSHIFT_AMPLITUDE 0.0   /* additional shift for wave amplitude */
 #define VSCALE_ENERGY 3.0      /* additional scaling factor for color scheme P_3D_ENERGY */
 #define PHASE_FACTOR 20.0      /* factor in computation of phase in color scheme P_3D_PHASE */
 #define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
@@ -287,6 +290,7 @@
 #define SHADE_2D 0       /* set to 1 to add pseudo-3d shading effect */ 
 #define SHADE_SCALE_2D 0.05  /* lower value increases sensitivity of shading */
 #define N_SOURCES 1                     /* number of sources, for option draw_sources */
+#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
 /* end of constants only used by sub_wave and sub_maze */
 
 

wave_billiard.c

@@ -57,8 +57,6 @@
 
 #define WINWIDTH 	1920  /* window width */
 #define WINHEIGHT 	1150  /* window height */
-// #define NX 1920          /* number of grid points on x axis */
-// #define NY 1150          /* number of grid points on y axis */
 #define NX 3840          /* number of grid points on x axis */
 #define NY 2300          /* number of grid points on y axis */
 
@@ -74,20 +72,21 @@
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 20        /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 6        /* choice of domain shape, see list in global_pdes.c */
 
-#define CIRCLE_PATTERN 1   /* pattern of circles or polygons, see list in global_pdes.c */
+#define CIRCLE_PATTERN 202   /* pattern of circles or polygons, see list in global_pdes.c */
 
 #define COMPARISON 0        /* set to 1 to compare two different patterns (beta) */
-#define B_DOMAIN_B 20       /* second domain shape, for comparisons */
+#define B_DOMAIN_B 20      /* second domain shape, for comparisons */
 #define CIRCLE_PATTERN_B 0  /* second pattern of circles or polygons */
 
-#define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
+#define P_PERCOL 0.15       /* probability of having a circle in C_RAND_PERCOL arrangement */
-#define NPOISSON 1000        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define NPOISSON 1000       /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.5    /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
-#define LAMBDA 0.2	    /* parameter controlling the dimensions of domain */
+#define LAMBDA 0.3	    /* parameter controlling the dimensions of domain */
-#define MU 0.012            /* parameter controlling the dimensions of domain */
+#define MU 1.0              /* parameter controlling the dimensions of domain */
 #define MU_B 1.0            /* parameter controlling the dimensions of domain */
 #define NPOLY 6             /* number of sides of polygon */
 #define APOLY -0.666666666666          /* angle by which to turn polygon, in units of Pi/2 */ 
@@ -96,9 +95,10 @@
 #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 24            /* number of grid point for grid of disks */
+#define NGRIDX 60           /* number of grid point for grid of disks */
-#define NGRIDY 40           /* number of grid point for grid of disks */
+#define NGRIDY 25           /* number of grid point for grid of disks */
-#define WALL_WIDTH 0.1      /* width of wall separating lenses */
+#define WALL_WIDTH 0.05      /* width of wall separating lenses */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 
 #define X_SHOOTER -0.2
 #define Y_SHOOTER -0.6
@@ -116,7 +116,7 @@
 
 /* Physical parameters of wave equation */
 
-#define TWOSPEEDS 0         /* set to 1 to replace hardcore boundary by medium with different speed */
+#define TWOSPEEDS 1         /* set to 1 to replace hardcore boundary by medium with different speed */
 #define OSCILLATE_LEFT 0     /* set to 1 to add oscilating boundary condition on the left */
 #define OSCILLATE_TOPBOT 0   /* set to 1 to enforce a planar wave on top and bottom boundary */
 #define OSCILLATION_SCHEDULE 61  /* oscillation schedule, see list in global_pdes.c */
@@ -128,8 +128,10 @@
 #define AMPLITUDE 0.5      /* amplitude of periodic excitation */ 
 #define ACHIRP 0.25        /* acceleration coefficient in chirp */
 #define DAMPING 0.0        /* damping of periodic excitation */
-#define COURANT 0.2       /* Courant number */
+#define COURANT 0.12       /* Courant number */
-#define COURANTB 0.0       /* Courant number in medium B */
+// #define COURANT 0.07       /* Courant number */
+#define COURANTB 0.08       /* Courant number in medium B */
+// #define GAMMA 1.0e-7          /* damping factor in wave equation */
 #define GAMMA 0.0          /* damping factor in wave equation */
 #define GAMMAB 0.0         /* damping factor in wave equation */
 #define GAMMA_SIDES 1.0e-4      /* damping factor on boundary */
@@ -144,16 +146,16 @@
 /* For similar wave forms, COURANT^2*GAMMA should be kept constant */
 
 #define ADD_OSCILLATING_SOURCE 1        /* set to 1 to add an oscillating wave source */
-#define OSCILLATING_SOURCE_PERIOD 8     /* period of oscillating source */
+#define OSCILLATING_SOURCE_PERIOD 3     /* period of oscillating source */
 #define ALTERNATE_OSCILLATING_SOURCE 1  /* set to 1 to alternate sign of oscillating source */
-#define N_SOURCES 2                     /* number of sources, for option draw_sources */
+#define N_SOURCES 1                     /* number of sources, for option draw_sources */
 
 #define ADD_WAVE_PACKET_SOURCES 0       /* set to 1 to add several sources emitting wave packets */
 #define WAVE_PACKET_SOURCE_TYPE 3       /* type of wave packet sources */
 #define N_WAVE_PACKETS 5                /* number of wave packets */
 #define WAVE_PACKET_RADIUS 50           /* radius of wave packets */
 
-#define USE_INPUT_TIMESERIES 1          /* set to 1 to use a time series (Morse code) as input * /
+#define USE_INPUT_TIMESERIES 0          /* set to 1 to use a time series (Morse code) as input * /
 
 /* Boundary conditions, see list in global_pdes.c  */
 
@@ -161,9 +163,10 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 1600       /* number of frames of movie */
+#define NSTEPS 1000      /* number of frames of movie */
-#define NVID 8            /* number of iterations between images displayed on screen */
+// #define NSTEPS 2000         /* number of frames of movie */
-#define NSEG 1000         /* number of segments of boundary */
+#define NVID 10             /* number of iterations between images displayed on screen */
+#define NSEG 1000           /* number of segments of boundary */
 #define INITIAL_TIME 0      /* time after which to start saving frames */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
 #define PRINT_SPEED 0       /* print speed of moving source */
@@ -179,7 +182,7 @@
 
 /* Parameters of initial condition */
 
-#define INITIAL_AMP 1.5             /* amplitude of initial condition */
+#define INITIAL_AMP 0.75             /* amplitude of initial condition */
 #define INITIAL_VARIANCE 0.00001    /* variance of initial condition */
 #define INITIAL_WAVELENGTH  0.025   /* wavelength of initial condition */
 
@@ -192,7 +195,7 @@
 /* Color schemes */
 
 #define COLOR_PALETTE 11     /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 14     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE_B 13   /* Color palette, see list in global_pdes.c  */
 
 #define BLACK 1          /* background */
 
@@ -202,15 +205,19 @@
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
 #define PHASE_FACTOR 1.0       /* factor in computation of phase in color scheme P_3D_PHASE */
 #define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
+// #define COLOR_PHASE_SHIFT 0.0   /* phase shift of color scheme, in units of Pi */
 #define ATTENUATION 0.0   /* exponential attenuation coefficient of contrast with time */
-#define E_SCALE 55.0      /* scaling factor for energy representation */
+#define VSHIFT_AMPLITUDE 0.0   /* additional shift for wave amplitude */
-#define LOG_SCALE 1.0     /* scaling factor for energy log representation */
+#define VSCALE_AMPLITUDE 0.5   /* additional scaling factor for wave amplitude */
-#define LOG_SHIFT 0.5     /* shift of colors on log scale */
+#define E_SCALE 50.0       /* scaling factor for energy representation */
-#define FLUX_SCALE 5.0e3    /* scaling factor for energy flux represtnation */
+#define LOG_SCALE 0.75     /* scaling factor for energy log representation */
-#define AVRG_E_FACTOR 0.95   /* controls time window size in P_AVERAGE_ENERGY scheme */
+#define LOG_SHIFT 0.75     /* shift of colors on log scale */
+#define FLUX_SCALE 250.0    /* scaling factor for energy flux represtnation */
+#define AVRG_E_FACTOR 0.75   /* controls time window size in P_AVERAGE_ENERGY scheme */
 #define RESCALE_COLOR_IN_CENTER 0   /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
-#define FADE_IN_OBSTACLE 0      /* set to 1 to fade color inside obstacles */
+#define FADE_IN_OBSTACLE 1      /* set to 1 to fade color inside obstacles */
 #define SHADE_2D 1       /* set to 1 to add pseudo-3d shading effect */ 
+// #define SHADE_SCALE_2D 0.25  /* lower value increases sensitivity of shading */
 #define SHADE_SCALE_2D 0.05  /* lower value increases sensitivity of shading */
 
 #define COLORHUE 260     /* initial hue of water color for scheme C_LUM */
@@ -220,24 +227,24 @@
 #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 0    /* set to 1 to plot the color scheme */
+#define DRAW_COLOR_SCHEME 1    /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 2.0     /* scale of color scheme bar */
+#define COLORBAR_RANGE 1.4     /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 0.8   /* scale of color scheme bar for 2nd part */
+#define COLORBAR_RANGE_B 1.2   /* scale of color scheme bar for 2nd part */
 #define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
 #define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
 #define CIRC_COLORBAR_B 0       /* set to 1 to draw circular color scheme */
 
-#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 0 /* set to 1 to draw wave profile vertically */
 #define VERTICAL_WAVE_PROFILE 1 /* set to 1 to draw wave profile vertically */
 #define WAVE_PROFILE_X 1.9      /* value of x to sample wave profile */
-#define WAVE_PROFILE_Y 0.075    /* value of y to sample wave profile */
+#define WAVE_PROFILE_Y 0.12     /* value of y to sample wave profile */
 #define PROFILE_AT_BOTTOM 1     /* draw wave profile at bottom instead of top */
-#define AVERAGE_WAVE_PROFILE 1  /* set to 1 to draw time-average of wave profile squared*/
+#define AVERAGE_WAVE_PROFILE 0  /* set to 1 to draw time-average of wave profile squared*/
 #define DRAW_WAVE_TIMESERIES 0  /* set to 1 to draw a time series of the wave, 2 to also draw it at the top */
 #define TIMESERIES_NVALUES 400  /* number of values plotted in time series */
 #define SAVE_TIME_SERIES 0      /* set to 1 to save wave time series at a point */
-#define DRAW_WAVE_SOURCE 1      /* set to 1 to draw source of wave at (wave_source_x, wave_source_y), set to 2 to also draw focus */
+#define DRAW_WAVE_SOURCE 0      /* set to 1 to draw source of wave at (wave_source_x, wave_source_y), set to 2 to also draw focus */
 
 #define MESSAGE_LDASH 14         /* length of dash for Morse code message */
 #define MESSAGE_LDOT 8           /* length of dot for Morse code message */
@@ -265,6 +272,7 @@
 
 #define MEAN_FLUX (PLOT == P_TOTAL_ENERGY_FLUX)||(PLOT_B == P_TOTAL_ENERGY_FLUX)
 #define REFRESH_IOR ((IOR == IOR_PERIODIC_WELLS_ROTATING)||(IOR == IOR_PERIODIC_WELLS_ROTATING_LARGE))
+#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
 
 double light[2] = {0.40824829, 0.816496581};   /* location of light source for SHADE_2D option*/
 
@@ -588,10 +596,10 @@ void draw_color_bar_palette(int plot, double range, int palette, int circular, i
 
 void animation()
 {
-    double time, scale, ratio, startleft[2], startright[2], sign = 1.0, r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle = 0.0, x1, sign1 = 1.0, ior_angle = 0.0, omega, phase_shift, vshift, dsource, finv; 
+    double time, scale, ratio, startleft[2], startright[2], sign[N_SOURCES], r2, xy[2], fade_value, yshift, speed = 0.0, a, b, c, x, y, angle = 0.0, x1, ior_angle = 0.0, omega, phase_shift, vshift, dsource, finv, source_amp; 
     double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *average_energy[NX], *color_scale[NX], *total_flux, *tcc_table[NX], *tgamma_table[NX], *fade_table;
     short int *xy_in[NX];
-    int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q, first_source = 1, imin, imax, ij[2];
+    int i, j, k, s, sample_left[2], sample_right[2], period = 0, fade, source_counter = 0, p, q, first_source = 1, imin, imax, ij[2], source, source_period, source_shift;
 //     static int image_counter = 0;
     int image_counter = 0;
     long int wave_value;
@@ -628,6 +636,8 @@ void animation()
         init_wave_packets(packet, WAVE_PACKET_RADIUS);
     }
     
+    for (i=0; i<N_SOURCES; i++) sign[i] = 1.0;
+    
     /* initialise positions and radii of circles */
     if ((B_DOMAIN == D_CIRCLES)||(B_DOMAIN == D_CIRCLES_IN_RECT)) ncircles = init_circle_config(circles);
     else if (B_DOMAIN == D_POLYGONS) ncircles = init_polygon_config(polygons);
@@ -695,6 +705,8 @@ void animation()
         printf("imin = %i, imax = %i\n", imin, imax);
     }
     
+    if (XYIN_INITIALISED) init_xyin_from_image(xy_in);
+    
 //     isospectral_initial_point(0.2, 0.0, startleft, startright);    /* for isospectral billiards */
 //     homophonic_initial_point(0.5, -0.25, 1.5, -0.25, startleft, startright);
 //     homophonic_initial_point(0.5, -0.25, 1.5, -0.25, startleft, startright);
@@ -709,7 +721,7 @@ void animation()
     init_ior_2d(xy_in, tcc_table, tgamma_table, ior_angle);
     if (FADE_IN_OBSTACLE) init_fade_table(tcc_table, fade_table); 
 
-//     init_circular_wave(-1.5, 0.0, phi, psi, xy_in);
+//     init_circular_wave(-LAMBDA, 0.0, phi, psi, xy_in);
 //     x = XMIN + (XMAX - XMIN)*rand()/RAND_MAX;
 //     y = YMIN + (YMAX - YMIN)*rand()/RAND_MAX;
 //     init_circular_wave(0.0, -0.8, phi, psi, xy_in);
@@ -760,7 +772,7 @@ void animation()
     blank();
     glColor3f(0.0, 0.0, 0.0);
 //     draw_wave(phi, psi, xy_in, 1.0, 0, PLOT);
-    if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 1, 0, 1.0);
+    if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 1, 0, 1.0);
     else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, 1.0, 0, PLOT, COLOR_PALETTE, 0, 1.0);
 
 //     draw_billiard(0, 1.0);
@@ -797,7 +809,7 @@ void animation()
         else scale = 1.0;
 
 //         draw_wave(phi, psi, xy_in, scale, i, PLOT);
-        if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 0, 1.0);
+        if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 0, 1.0);
         else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
         for (j=0; j<NVID; j++) 
         {
@@ -811,7 +823,6 @@ void animation()
                 if ((j == 0)&&(i%10 == 0)) printf("Frame %i of %i\n", i, NSTEPS);
 //                 fprintf(time_series_right, "%.15f\n", phi[sample_right[0]][sample_right[1]]);
             }
-//             if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
         }
         
 //         draw_billiard(0, 1.0);
@@ -819,21 +830,16 @@ void animation()
         if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, fade, fade_value);
         
         /* add oscillating waves */
-        wave_source_x[0] = -1.6;
+        wave_source_x[0] = XMIN + 2.0*(XMAX - XMIN)*(double)i/(double)NSTEPS;
-        wave_source_y[0] = 0.5;
+        wave_source_y[0] = 0.3;
-        wave_source_x[1] = -1.6;
-        wave_source_y[1] = -0.5;
         if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
         {
-            if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
+            if (ALTERNATE_OSCILLATING_SOURCE) sign[0] = -sign[0];
-            add_circular_wave(sign*INITIAL_AMP, wave_source_x[0], wave_source_y[0], phi, psi, xy_in);
+            for (source = 0; source < N_SOURCES; source++)
-        }
+                add_circular_wave(-sign[0]*INITIAL_AMP, wave_source_x[source], wave_source_y[source], phi, psi, xy_in);
-        if ((ADD_OSCILLATING_SOURCE)&&(i%(OSCILLATING_SOURCE_PERIOD/3) == 1))
-        {
-            if (ALTERNATE_OSCILLATING_SOURCE) sign1 = -sign1;
-            add_circular_wave(sign1*INITIAL_AMP/sqrt(3.0), wave_source_x[1], wave_source_y[1], phi, psi, xy_in);
         }
 
+        
         if (ADD_WAVE_PACKET_SOURCES) add_wave_packets(phi, psi, xy_in, packet, i, WAVE_PACKET_RADIUS, 1, 4, 1);
         if (PRINT_SPEED) print_speed(speed, 0, 1.0);
         if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0);
@@ -862,7 +868,7 @@ void animation()
             {
 //                 draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
                 if (HIGHRES) 
-                    draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 1, 0, 1.0);
+                    draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 1, 0, 1.0);
                 else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
 //                 draw_billiard(0, 1.0);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0);  
@@ -893,7 +899,7 @@ void animation()
         if (DOUBLE_MOVIE) 
         {
 //             draw_wave(phi, psi, xy_in, scale, i, PLOT);
-            if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 0, 1.0);
+            if (HIGHRES) draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 0, 1.0);
             else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
 //             draw_billiard(0, 1.0);
             if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 0, 1.0);   
@@ -906,7 +912,7 @@ void animation()
         {
             fade_value = 1.0 - (double)i/(double)MID_FRAMES;
             if (HIGHRES) 
-                draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1, fade_value);
+                draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1, fade_value);
             else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
 //             draw_billiard(1, fade_value);
             if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, CIRC_COLORBAR, 1, fade_value); 
@@ -919,7 +925,7 @@ void animation()
         {
 //             draw_wave(phi, psi, xy_in, scale, i, PLOT_B);
             if (HIGHRES) 
-                draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, 0, 1.0);
+                draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, 0, 1.0);
             else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
 //             draw_billiard(0, 1.0);
             if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 0, 1.0); 
@@ -932,7 +938,7 @@ void animation()
             {
                 fade_value = 1.0 - (double)i/(double)END_FRAMES;
                 if (HIGHRES) 
-                    draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, 1, fade_value);
+                    draw_wave_highres_palette(2, phi, psi, total_energy, average_energy, fade_table, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, 1, fade_value);
                 else draw_wave_epalette(phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
 //                 draw_billiard(1, fade_value);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, CIRC_COLORBAR_B, 1, fade_value); 

wave_common.c

@@ -297,14 +297,22 @@ void init_wave_flat( double *phi[NX], double *psi[NX], short int * xy_in[NX])
     int i, j;
     double xy[2], dist2;
 
+    if (!XYIN_INITIALISED) 
+        for (i=0; i<NX; i++) 
+        {
+            if (i%100 == 0) printf("Table xy_in - Initialising column %i of %i\n", i, NX);
+            for (j=0; j<NY; j++)
+            {
+                ij_to_xy(i, j, xy);
+                xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
+            }
+        }
+
     for (i=0; i<NX; i++) {
-        if (i%100 == 0) printf("Wave and table xy_in - Initialising column %i of %i\n", i, NX);
+        if (i%100 == 0) printf("Wave - Initialising column %i of %i\n", i, NX);
         for (j=0; j<NY; j++)
         {
-            ij_to_xy(i, j, xy);
+            phi[i][j] = 0.0;
-	    xy_in[i][j] = xy_in_billiard(xy[0],xy[1]);
-//             if ((i%10 == 0)&&(j == NY/2)) printf ("xy_in[%i][%i] = %i\n", i, j, xy_in[i][j]);
-	    phi[i][j] = 0.0;
             psi[i][j] = 0.0;
         }
     }
@@ -1042,7 +1050,7 @@ void draw_wave_epalette(double *phi[NX], double *psi[NX], double *total_energy[N
 }
 
 
-double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, double rgb[3])
+double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *total_flux, double *color_scale[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, double rgb[3])
 /* compute value of wave and color */
 {
     int k;
@@ -1052,7 +1060,7 @@ double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_
         case (P_AMPLITUDE):
         {
             value = phi[i][j];
-            color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
+            color_scheme_palette(COLOR_SCHEME, palette, VSHIFT_AMPLITUDE + VSCALE_AMPLITUDE*value, scale, time, rgb);
             break;
         }
         case (P_ENERGY):
@@ -1125,7 +1133,8 @@ double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_
             if (arg < 0.0) arg += DPI;
             value = arg/DPI;
             color_scheme_palette(C_ONEDIM_LINEAR, palette, value, 1.0, 1, rgb);
-            flux_factor = tanh(mod*E_SCALE);
+//             flux_factor = tanh(mod*E_SCALE);
+            flux_factor = tanh(mod*FLUX_SCALE);
             for (k=0; k<3; k++) rgb[k] *= flux_factor;
             break;
         }
@@ -1150,6 +1159,13 @@ double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_
         }
     }
     
+    if (RESCALE_COLOR_IN_CENTER)
+    {
+        value *= color_scale[i][j];
+        for (k=0; k<3; k++) rgb[k] = 1.0 + color_scale[i][j]*(rgb[k]-1.0);
+//         for (k=0; k<3; k++) rgb[k] = 0.5 + color_scale[i][j]*(rgb[k]-0.5);
+    }
+    
     return(value);
 }
 
@@ -1160,7 +1176,7 @@ void draw_wave_profile_horizontal(double *values, int size, int average, int pnb
     int i, k, ij[2];
     double deltav, a, b, y;
     static int imin, imax, jmin, jmax, jmid, jval, d, first = 1;
-    static double vmin[2], vmax[2], deltaj, ymin, average_vals[NX];
+    static double vmin[2], vmax[2], deltaj, ymin, average_vals[2*NX];
     
     if (first)
     {
@@ -1184,7 +1200,7 @@ void draw_wave_profile_horizontal(double *values, int size, int average, int pnb
         d = (jmax - jmin)/10 + 1;
         deltaj = (double)(jmax - jmin - 2*d); 
         ymin = (double)(jmin + d);
-        for (i=0; i<NX; i++) average_vals[i] = 0.0;
+        for (i=0; i<2*NX; i++) average_vals[i] = 0.0;
         for (k=0; k<2; k++)
         {
             vmin[k] = 1.0e10;
@@ -1197,12 +1213,12 @@ void draw_wave_profile_horizontal(double *values, int size, int average, int pnb
     if ((average)&&(!fade))
     {
         for (i=imin; i<imax; i+=size)
-            average_vals[i] = 0.99*average_vals[i] + 0.01*values[i*NY+jval]*values[i*NY+jval];
+            average_vals[pnb*NX+i] = 0.99*average_vals[pnb*NX+i] + 0.01*values[i*NY+jval]*values[i*NY+jval];
         for (i=imin; i<imax; i+=size)
-            if (average_vals[i] > vmax[pnb]) vmax[pnb] = average_vals[i];
+            if (average_vals[pnb*NX+i] > vmax[pnb]) vmax[pnb] = average_vals[pnb*NX+i];
         
     }
-    else /*if (!fade)*/ 
+    else if (!fade) 
     {
         for (i=imin; i<imax; i+=size)
         {
@@ -1213,11 +1229,11 @@ void draw_wave_profile_horizontal(double *values, int size, int average, int pnb
         }
     }
 //     if (vmax <= vmin) vmax = vmin + 0.01;
-    if ((vmin[pnb] < 0.0)&&(vmax[pnb] > 0.0))
+//     if ((vmin[pnb] < 0.0)&&(vmax[pnb] > 0.0))
-    {
+//     {
-        if (vmax[pnb] > -vmin[pnb]) vmin[pnb] = -vmax[pnb];
+//         if (vmax[pnb] > -vmin[pnb]) vmin[pnb] = -vmax[pnb];
-        else if (vmax[pnb] < -vmin[pnb]) vmax[pnb] = -vmin[pnb];
+//         else if (vmax[pnb] < -vmin[pnb]) vmax[pnb] = -vmin[pnb];
-    }
+//     }
 //     printf("vmin = %.3lg, vmax = %.3lg\n", vmin, vmax);
     deltav = vmax[pnb]-vmin[pnb];
     if (deltav <= 0.0) deltav = 0.01;
@@ -1232,7 +1248,7 @@ void draw_wave_profile_horizontal(double *values, int size, int average, int pnb
     for (i=imin; i<imax; i+=size)
     {
 //         y = a*values[i*NY+jmin] + b;
-        if (average) y = a*average_vals[i] + b;
+        if (average) y = a*average_vals[pnb*NX+i] + b;
         else y = a*values[i*NY+jval] + b;
         if (y > ymin) glVertex2d((double)i, y);
     }
@@ -1250,6 +1266,8 @@ void draw_wave_profile_horizontal(double *values, int size, int average, int pnb
     {
         vmax[pnb] *= 0.99;
         vmin[pnb] *= 0.99;
+//         vmax[pnb] *= 0.95;
+//         vmin[pnb] *= 0.95;
     }
 }
 
@@ -1260,7 +1278,7 @@ void draw_wave_profile_vertical(double *values, int size, int average, int pnb,
     int j, k, ij[2];
     double deltav, a, b, x;
     static int imin, imax, ival, jmin, jmax, imid, d, first = 1;
-    static double vmin[2], vmax[2], deltai, xmin, average_vals[NY];
+    static double vmin[2], vmax[2], deltai, xmin, average_vals[2*NY];
     
     if (first)
     {        
@@ -1276,7 +1294,7 @@ void draw_wave_profile_vertical(double *values, int size, int average, int pnb,
         d = (imax - imin)/10 + 1;
         deltai = (double)(imax - imin - 2*d); 
         xmin = (double)(imin + d);
-        for (j=0; j<NY; j++) average_vals[j] = 0.0;
+        for (j=0; j<2*NY; j++) average_vals[j] = 0.0;
         for (k=0; k<2; k++)
         {
             vmin[k] = 1.0e10;
@@ -1298,9 +1316,9 @@ void draw_wave_profile_vertical(double *values, int size, int average, int pnb,
     if ((average)&&(!fade))
     {
         for (j=jmin; j<jmax; j+=size)
-            average_vals[j] = 0.99*average_vals[j] + 0.01*values[ival*NY+j]*values[ival*NY+j];
+            average_vals[pnb*NY+j] = 0.99*average_vals[pnb*NY+j] + 0.01*values[ival*NY+j]*values[ival*NY+j];
         for (j=jmin; j<jmax; j+=size)
-            if (average_vals[j] > vmax[pnb]) vmax[pnb] = average_vals[j];
+            if (average_vals[pnb*NY+j] > vmax[pnb]) vmax[pnb] = average_vals[pnb*NY+j];
         
     }
     
@@ -1345,7 +1363,7 @@ void draw_wave_profile_vertical(double *values, int size, int average, int pnb,
     glBegin(GL_LINE_STRIP);
     for (j=jmin; j<jmax; j+=size)
     {
-        if (average) x = a*average_vals[j] + b;
+        if (average) x = a*average_vals[pnb*NY+j] + b;
         else x = a*values[ival*NY+j] + b;
 //         else x = a*values[imin*NY+j] + b;
         if (x >= xmin) glVertex2d(x, (double)j);
@@ -1722,25 +1740,33 @@ void init_fade_table(double *tcc_table[NX], double *fade_table)
     
     for (i=0; i<NX; i++)
         for (j=0; j<NY; j++)
-            fade_table[i*NY + j] = tcc_table[i][j]/courant2;
+            fade_table[i*NY + j] = pow(tcc_table[i][j]/courant2, 0.25);
+//             fade_table[i*NY + j] = tcc_table[i][j]/courant2;
 }
 
-void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *fade_table, double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, int pnumber, int fade, double fade_value)
+void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *average_energy[NX], double *fade_table, double *total_flux, double *color_scale[NX], short int *xy_in[NX], double scale, int time, int plot, int palette, int pnumber, int fade, double fade_value)
 /* same as draw_wave_highres, but with color scheme option */
 {
     int i, j, k, iplus, iminus, jplus, jminus;
     double value, rgb[3], xy[2], x1, y1, x2, y2, velocity, energy, gradientx, gradienty, arg, mod, flux_factor, gx, gy, mgx, mgy, factor, norm, pscal, ca, vscale2;
 //     double vmin, vmax, deltav;
-    static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
+    static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX), col_ratio;
+    static int first = 1;
     double *values, *shade, *rgbvals;
         
+    if (first)
+    {
+        col_ratio = 1.0/1.01;
+        first = 0;
+    }
+    
     values = (double *)malloc(NX*NY*sizeof(double));
     rgbvals = (double *)malloc(3*NX*NY*sizeof(double));
     #pragma omp parallel for private(i,j,rgb)
     for (i=0; i<NX; i+=size)
         for (j=0; j<NY; j+=size)
         {
-            values[i*NY+j] = wave_value(i, j, phi, psi, total_energy, average_energy, total_flux, xy_in, scale, time, plot, palette, rgb);
+            values[i*NY+j] = wave_value(i, j, phi, psi, total_energy, average_energy, total_flux, color_scale, xy_in, scale, time, plot, palette, rgb);
             rgbvals[i*NY+j] = rgb[0];
             rgbvals[NX*NY+i*NY+j] = rgb[1];
             rgbvals[2*NX*NY+i*NY+j] = rgb[2];
@@ -1761,7 +1787,10 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
                 pscal = -gradientx*light[0] - gradienty*light[1] + SHADE_SCALE_2D;
                 ca = pscal/norm;
                 ca = (ca + 1.0)*0.4 + 0.2;
-                for (k=0; k<3; k++) rgbvals[k*NX*NY+i*NY+j] *= ca;
+                if (RESCALE_COLOR_IN_CENTER) for (k=0; k<3; k++) 
+                    rgbvals[k*NX*NY+i*NY+j] = 1.0 + (ca*rgbvals[k*NX*NY+i*NY+j] - 1.0)*(color_scale[i][j]+0.01)*col_ratio;
+//                     rgbvals[k*NX*NY+i*NY+j] = 0.5 + (ca*rgbvals[k*NX*NY+i*NY+j] - 0.5)*(color_scale[i][j]+0.01)*col_ratio;
+                else for (k=0; k<3; k++) rgbvals[k*NX*NY+i*NY+j] *= ca;
             }
         }
     }

wave_comparison.c

@@ -93,6 +93,7 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 300        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.25   /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 0 /* set to 1 to randomize angle of polygons */
 #define RANDOM_POLY_ANGLE_B 0 /* set to 1 to randomize angle of polygons */
 
@@ -207,6 +208,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.0   /* additional shift for wave amplitude */
+#define VSCALE_AMPLITUDE 0.5   /* additional scaling factor for wave amplitude */
 #define E_SCALE 100.0     /* scaling factor for energy representation */
 #define LOG_SCALE 1.0     /* scaling factor for energy log representation */
 #define LOG_SHIFT 2.0     /* shift of colors on log scale */
@@ -255,6 +258,7 @@
 #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
 #define DRAW_WAVE_TIMESERIES 0  /* set to 1 to draw a time series of the wave */
 #define WALL_WIDTH 0.1      /* width of wall separating lenses */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 #define OSCIL_YMAX 0.35      /* defines oscillation range */
 #define MESSAGE_LDASH 14         /* length of dash for Morse code message */
 #define MESSAGE_LDOT 8           /* length of dot for Morse code message */
@@ -277,6 +281,7 @@
 #define SHADE_2D 0       /* set to 1 to add pseudo-3d shading effect */ 
 #define SHADE_SCALE_2D 0.05  /* lower value increases sensitivity of shading */
 #define N_SOURCES 1                     /* number of sources, for option draw_sources */
+#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
 double light[2] = {0.40824829, 0.816496581};   /* location of light source for SHADE_2D option*/
 /* end of constants only used by sub_wave and sub_maze */
 

wave_energy.c

@@ -73,6 +73,7 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 300        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.25   /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 #define RANDOM_POLY_ANGLE_B 0 /* set to 1 to randomize angle of polygons */
 
@@ -176,6 +177,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.0   /* additional shift for wave amplitude */
+#define VSCALE_AMPLITUDE 0.5   /* additional scaling factor for wave amplitude */
 #define E_SCALE 500.0     /* scaling factor for energy representation */
 #define LOG_SCALE 1.5     /* scaling factor for energy log representation */
 #define LOG_SHIFT 1.0     /* shift of colors on log scale */
@@ -229,6 +232,7 @@
 #define VERTICAL_WAVE_PROFILE 0 /* set to 1 to draw wave profile vertically */
 #define DRAW_WAVE_TIMESERIES 0  /* set to 1 to draw a time series of the wave */
 #define WALL_WIDTH 0.1      /* width of wall separating lenses */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 #define OSCIL_YMAX 0.35      /* defines oscillation range */
 #define MESSAGE_LDASH 14         /* length of dash for Morse code message */
 #define MESSAGE_LDOT 8           /* length of dot for Morse code message */
@@ -251,6 +255,7 @@
 #define SHADE_2D 0       /* set to 1 to add pseudo-3d shading effect */ 
 #define SHADE_SCALE_2D 0.05  /* lower value increases sensitivity of shading */
 #define N_SOURCES 1                     /* number of sources, for option draw_sources */
+#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
 double light[2] = {0.40824829, 0.816496581};   /* location of light source for SHADE_2D option*/
 /* end of constants only used by sub_wave and sub_maze */
 

wave_sphere.c

@@ -90,6 +90,7 @@
 
 #define P_PERCOL 0.25       /* probability of having a circle in C_RAND_PERCOL arrangement */
 #define NPOISSON 1000        /* number of points for Poisson C_RAND_POISSON arrangement */
+#define PDISC_FACTOR 3.25   /* controls density of Poisson disc process (default: 3.25) */
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 0.75	    /* parameter controlling the dimensions of domain */
@@ -105,6 +106,7 @@
 #define NGRIDX 30            /* number of grid point for grid of disks */
 #define NGRIDY 18            /* number of grid point for grid of disks */
 #define WALL_WIDTH 0.6      /* width of wall separating lenses */
+#define RADIUS_FACTOR 0.3   /* controls inner radius for C_RING arrangements */
 
 #define X_SHOOTER -0.2
 #define Y_SHOOTER -0.6
@@ -260,6 +262,7 @@
 #define SCALE 0          /* set to 1 to adjust color scheme to variance of field */
 #define SLOPE 1.0        /* sensitivity of color on wave amplitude */
 #define VSCALE_AMPLITUDE 50.0  /* additional scaling factor for color scheme P_3D_AMPLITUDE */
+#define VSHIFT_AMPLITUDE 0.0   /* additional shift for wave amplitude */
 #define VSCALE_ENERGY 5.0     /* additional scaling factor for color scheme P_3D_ENERGY */
 #define PHASE_FACTOR 20.0      /* factor in computation of phase in color scheme P_3D_PHASE */
 #define PHASE_SHIFT 0.0      /* shift of phase in color scheme P_3D_PHASE */
@@ -325,6 +328,7 @@
 #define WAVE_PACKET_SHIFT 200.0     /* time shift between wave packets (in oscillation periods) */
 #define FADE_IN_OBSTACLE 0      /* set to 1 to fade color inside obstacles */
 #define N_SOURCES 1                     /* number of sources, for option draw_sources */
+#define XYIN_INITIALISED (B_DOMAIN == D_IMAGE)
 /* end of constants only used by sub_wave and sub_maze */
 
 /* For debugging purposes only */