Add files via upload

global_3d.c

@@ -1,21 +1,5 @@
 /* global variables and definitions used by sub_wave_3d.c */
 
-/* plot types used by wave_3d */
-
-#define P_3D_AMPLITUDE  101     /* height/color depends on amplitude - DEPRECATED, instead use set SHADE_3D to 0 */
-#define P_3D_ANGLE 102          /* height/color depends on angle with fixed direction - TODO */
-#define P_3D_AMP_ANGLE 103      /* height/color depends on amplitude, luminosity depends on angle */
-#define P_3D_ENERGY 104         /* height/color depends on energy, luminosity depends on angle */
-#define P_3D_LOG_ENERGY 105     /* height/color depends on logarithm of energy, luminosity depends on angle */
-#define P_3D_TOTAL_ENERGY 106   /* height/color depends on total energy over time, luminosity depends on angle */
-#define P_3D_LOG_TOTAL_ENERGY 107 /* height/color depends on log on total energy over time, luminosity depends on angle */
-#define P_3D_MEAN_ENERGY 108      /* height/color depends on energy averaged over time, luminosity depends on angle */
-#define P_3D_LOG_MEAN_ENERGY 109  /* height/color depends on log on energy averaged over time, luminosity depends on angle */
-
-#define P_3D_PHASE 111          /* phase of wave */
-#define P_3D_FLUX_INTENSITY 112    /* energy flux intensity */
-#define P_3D_FLUX_DIRECTION 113    /* energy flux direction */
-
 /* Choice of simulated reaction-diffusion equation in rde.c */
 
 #define E_HEAT 0            /* heat equation */
@@ -73,6 +57,8 @@
 #define DEM_EARTH 0         /* DEM of Earth */
 #define DEM_MARS 1          /* DEM of Mars */
 #define DEM_MOON 2          /* DEM of the Moon */
+#define DEM_VENUS 3         /* DEM of Venus */
+#define DEM_MERCURY 4       /* DEM of Mercury */
 
 /* macros to avoid unnecessary computations in 3D plots */
 
@@ -93,8 +79,8 @@
 
 #define COMPUTE_TOTAL_ENERGY ((ZPLOT == P_3D_TOTAL_ENERGY)||(CPLOT == P_3D_TOTAL_ENERGY)||(ZPLOT == P_3D_LOG_TOTAL_ENERGY)||(CPLOT == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT == P_3D_MEAN_ENERGY)||(CPLOT == P_3D_MEAN_ENERGY)||(ZPLOT == P_3D_LOG_MEAN_ENERGY)||(CPLOT == P_3D_LOG_MEAN_ENERGY)||(ZPLOT_B == P_3D_TOTAL_ENERGY)||(CPLOT_B == P_3D_TOTAL_ENERGY)||(ZPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(CPLOT_B == P_3D_LOG_TOTAL_ENERGY)||(ZPLOT_B == P_3D_MEAN_ENERGY)||(CPLOT_B == P_3D_MEAN_ENERGY)||(ZPLOT_B == P_3D_LOG_MEAN_ENERGY)||(CPLOT_B == P_3D_LOG_MEAN_ENERGY))
 
-#define PLANET ((B_DOMAIN == D_SPHERE_EARTH)||(B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON))
-#define OTHER_PLANET ((B_DOMAIN == D_SPHERE_MARS)||(B_DOMAIN == D_SPHERE_MOON))
+#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 NMAXCIRC_SPHERE 100     /* max number of circles on sphere */
 

global_ljones.c

@@ -13,10 +13,11 @@
 
 #define NMAXCIRCLES 100000       /* 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 100   /* max number of obstacles */
+#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 10     /* max number of partners in molecule */
 
 #define C_SQUARE 0          /* square grid of circles */
 #define C_HEX 1             /* hexagonal/triangular grid of circles */
@@ -46,6 +47,7 @@
 #define O_HLINE_HOLE_SPOKES 181    /* tips of spokes for S_HLINE_HOLE_SPOKES segment pattern */
 #define O_CIRCLE 4          /* one circle at the origin */
 #define O_FOUR_CIRCLES 5    /* four circles  */
+#define O_HEX 6             /* hexagonal lattice */
 
 /* pattern of additional repelling segments */
 #define S_RECTANGLE 0       /* segments forming a rectangle */
@@ -79,6 +81,8 @@
 #define S_AIRFOIL 24            /* exterior of an air foil */
 #define S_COANDA 25             /* wall for Coanda effect */
 #define S_COANDA_SHORT 26       /* shorter wall for Coanda effect */
+#define S_CYLINDER 27           /* walls at top and bottom, for cylindrical b.c. */
+#define S_TREE 28               /* Christmas tree(s) */
 
 /* particle interaction */
 
@@ -121,7 +125,7 @@
 #define TH_VERTICAL 0       /* only particles at the right of x = PARTIAL_THERMO_SHIFT are coupled */
 #define TH_INSEGMENT 1      /* only particles in region defined by segments are coupled */
 #define TH_INBOX 2          /* only particles in a given box are coupled */
-#define TH_LAYER 3          /* only particles above -LAMBDA are coupled */
+#define TH_LAYER 3          /* only particles above PARTIAL_THERMO_HEIGHT are coupled */
 #define TH_LAYER_TYPE2 4    /* only particles above highest type 2 particle are coupled */
 
 /* Gravity schedules */
@@ -221,6 +225,20 @@
 #define IP_X 0      /* color depends on x coordinate of initial position */
 #define IP_Y 1      /* color depends on y coordinate of initial position */
 
+/* Interaction types for polyatomic molecules */
+
+#define POLY_STAR 0     /* star-shaped graph (central molecule attracts outer ones) */
+#define POLY_ALL 1      /* all-to-all coupling */
+#define POLY_WATER 2    /* star-shaped with a 120° separation between anions */
+
+/* Background color schemes */
+
+#define BG_NONE 0       /* no background color */
+#define BG_DENSITY 1    /* background color depends on number of particles */
+#define BG_CHARGE 2     /* background color depends on charge density */
+#define BG_EKIN 3       /* background color depends on kinetic energy */
+#define BG_FORCE 4      /* background color depends on total force */
+
 /* Color schemes */
 
 #define C_LUM 0          /* color scheme modifies luminosity (with slow drift of hue) */
@@ -282,6 +300,9 @@ typedef struct
     double spin_freq;           /* angular frequency of spin-spin interaction */
     double color_hue;           /* color hue of particle, for P_INITIAL_POS plot type */
     int color_rgb[3];           /* RGB colors code of particle, for use in ljones_movie.c */
+    int partner[NMAXPARTNERS];  /* partners particle for option PAIR_PARTICLES */
+    short int npartners;        /* number of partner particles */
+    double partner_eqd[NMAXPARTNERS];   /* equilibrium distances between partners */
 } t_particle;
 
 typedef struct
@@ -290,6 +311,9 @@ typedef struct
     int particles[HASHMAX];     /* numbers of particles in cell */
     int nneighb;                /* number of neighbouring cells */
     int neighbour[9];           /* numbers of neighbouring cells */
+    double x1, y1, x2, y2;      /* coordinates of hashcell corners */
+    double hue1, hue2;          /* color hues */
+    double charge;              /* charge of fixed obstacles */
 } t_hashgrid;
 
 typedef struct
@@ -314,6 +338,7 @@ typedef struct
 {
     double xc, yc, radius;      /* center and radius of circle */
     short int active;           /* circle is active */
+    double charge;              /* charge of obstacle, for EM simulations */
 } t_obstacle;
 
 typedef struct
@@ -384,6 +409,7 @@ typedef struct
     double omega;               /* angular speed of obstacle */
     double bdry_fx, bdry_fy;    /* components of boundary force */
     double efield, bfield;      /* electric and magnetic field */
+    double prop;                /* proportion of types */
 } t_lj_parameters;
 
 int ncircles, nobstacles, nsegments, ngroups = 1, counter = 0;

global_pdes.c

@@ -66,6 +66,9 @@
 #define D_LSHAPE 50             /* L-shaped billiard (surface of genus 2) */
 #define D_WAVEGUIDE 51          /* wave guide */
 #define D_WAVEGUIDE_W 52        /* W-shaped wave guide */
+#define D_WAVEGUIDES_W 521      /* two W-shaped wave guides */
+#define D_WAVEGUIDES_COUPLED 522    /* two coupled wave guides */
+#define D_WAVEGUIDE_S 523       /* S-shaped wave guide */
 #define D_MAZE 53               /* maze */
 #define D_MAZE_CLOSED 54        /* closed maze */
 #define D_MAZE_CHANNELS 541     /* maze with two channels attached */
@@ -78,11 +81,21 @@
 #define D_LENSES_RING 59        /* several lenses forming a ring */
 #define D_MAZE_CIRCULAR 60      /* circular maze */
 #define D_LENS 61               /* symmetric lens made of circular faces */
+#define D_LENS_WALL 62          /* symmetric lens made of circular faces with separating wall (to use with IOR_LENS_WALL) */
+#define D_TWO_LENSES_WALL 63    /* two lenses separated by a wall with a hole (to use with IOR_LENS_WALL) */
+#define D_TWO_LENSES_OBSTACLE 64    /* two lenses with an obstacle in between (to use with IOR_LENS_OBSTACLE) */
+#define D_FRESNEL_ZONE_PLATE 65 /* Fresnel zone plate */
+#define D_FRESNEL_ZONE_PLATE_INV 66 /* Fresnel zone plate, with central hole */
+#define D_LENS_ROTATED 67       /* rotated lens */
+#define D_LENS_CONCAVE 68       /* biconcave lens (to use with IOR_LENS_CONCAVE) */
+#define D_LENS_CONVEX_CONCAVE 69    /* a convex and a biconcave lens (to use with IOR_LENS_CONVEX_CONCAVE) */
 
 #define D_WING 70               /* complement of wing-shaped domain */
 #define D_TESLA 71              /* Tesla valve */
 #define D_TESLA_FOUR 72         /* four Tesla valves */
 
+#define D_TREE 73               /* Christmas tree, to use with IOR_TREE */
+
 /* for wave_sphere.c */
 
 #define D_LATITUDE 80           /* strip between two latitudes */
@@ -93,6 +106,11 @@
 #define D_SPHERE_JULIA_CUBIC 85 /* Julia set for cubic polynomial on Riemann sphere */
 #define D_SPHERE_MARS 86        /* map of Mars */
 #define D_SPHERE_MOON 87        /* map of the Moon */
+#define D_SPHERE_VENUS 88       /* map of Venus */
+#define D_SPHERE_MERCURY 89     /* map of Mercury */
+#define D_SPHERE_MAZE 100       /* circular maze on the sphere */
+#define D_SPHERE_MAZE_SPIRAL 101 /* circular maze on the sphere with slanted walls */
+#define D_SPHERE_MAZE_WAVE 102  /* circular maze on the sphere with wavy walls */
 
 #define NMAXCIRCLES 10000       /* total number of circles/polygons (must be at least NCX*NCY for square grid) */
 #define NMAXPOLY 50000          /* maximal number of vertices of polygonal lines (for von Koch et al) */
@@ -152,8 +170,14 @@
 #define IOR_PERIODIC_WELLS_ROTATING_LARGE 7   /* periodic superposition rotating in time, larger area */
 #define IOR_POISSON_WELLS 8     /* wells located on a random Poisson disc process */
 #define IOR_PPP_WELLS 9         /* wells located on a Poisson point process */
+#define IOR_LENS_WALL 10        /* lens with separating wall (to use with D_LENS_WALL) */
+#define IOR_LENS_OBSTACLE 11    /* lens with separating wall (to use with D_TWO_LENSES_OBSTACLE) */
+#define IOR_LENS_CONCAVE 12     /* lens with separating wall (to use with D_LENS_CONCAVE) */
+#define IOR_LENS_CONVEX_CONCAVE 13     /* lens with separating wall (to use with D_LENS_CONVEX_CONCAVE) */
+#define IOR_TREE 14             /* Christmas tree, to use with D_TREE */
+#define IOR_WAVE_GUIDES_COUPLED 15  /* coupled wave guides */
 
-#define IOR_EARTH_DEM 10        /* digital elevation model (for waves on sphere) */
+#define IOR_EARTH_DEM 20        /* digital elevation model (for waves on sphere) */
 
 /* Boundary conditions */
 
@@ -207,6 +231,22 @@
 #define P_REAL 12          /* plot real part */
 #define P_IMAGINARY 13     /* plot imaginary part */
 
+/* plot types used by wave_3d */
+
+#define P_3D_AMPLITUDE  101     /* height/color depends on amplitude - DEPRECATED, instead use set SHADE_3D to 0 */
+#define P_3D_ANGLE 102          /* height/color depends on angle with fixed direction - TODO */
+#define P_3D_AMP_ANGLE 103      /* height/color depends on amplitude, luminosity depends on angle */
+#define P_3D_ENERGY 104         /* height/color depends on energy, luminosity depends on angle */
+#define P_3D_LOG_ENERGY 105     /* height/color depends on logarithm of energy, luminosity depends on angle */
+#define P_3D_TOTAL_ENERGY 106   /* height/color depends on total energy over time, luminosity depends on angle */
+#define P_3D_LOG_TOTAL_ENERGY 107 /* height/color depends on log on total energy over time, luminosity depends on angle */
+#define P_3D_MEAN_ENERGY 108      /* height/color depends on energy averaged over time, luminosity depends on angle */
+#define P_3D_LOG_MEAN_ENERGY 109  /* height/color depends on log on energy averaged over time, luminosity depends on angle */
+
+#define P_3D_PHASE 111          /* phase of wave */
+#define P_3D_FLUX_INTENSITY 112    /* energy flux intensity */
+#define P_3D_FLUX_DIRECTION 113    /* energy flux direction */
+
 
 /* Color schemes */
 
@@ -366,6 +406,8 @@ int npolyrect = NMAXPOLY;           /* actual number of polyrect */
 int npolyrect_rot = NMAXPOLY;       /* actual number of rotated polyrect */
 int npolyarc = NMAXPOLY;            /* actual number of arcs */
 
+short int input_signal[NSTEPS];     /* time-dependent source signal */
+
 t_circle circles[NMAXCIRCLES];      /* circular scatterers */
 t_polygon polygons[NMAXCIRCLES];    /* polygonal scatterers */
 t_vertex polyline[NMAXPOLY];        /* vertices of polygonal line */

heat.c

@@ -213,6 +213,11 @@
 #define N_WAVE_PACKETS 15               /* number of wave packets */
 #define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 #define OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
+#define MU_B 1.0           /* parameter controlling the dimensions of domain */
+#define GAMMA 0.0          /* damping factor in wave equation */
+#define GAMMAB 0.0         /* damping factor in wave equation */
+#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 */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"

lennardjones.c

@@ -59,15 +59,14 @@
 #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
 
 #define INITXMIN -2.0
-#define INITXMAX 2.03	/* x interval for initial condition */
-#define INITYMIN -1.1
-#define INITYMAX 1.1	/* y interval for initial condition */
-
-#define ADDXMIN -1.97
-#define ADDXMAX -0.8	/* x interval for adding particles */
-#define ADDYMIN -1.125
-#define ADDYMAX 1.125	/* y interval for adding particles */
+#define INITXMAX 2.0	/* x interval for initial condition */
+#define INITYMIN -1.125
+#define INITYMAX 1.125	/* y interval for initial condition */
 
+#define ADDXMIN 1.9
+#define ADDXMAX 2.0	/* x interval for adding particles */
+#define ADDYMIN -0.9
+#define ADDYMAX 0.9	/* y interval for adding particles */
 
 #define BCXMIN -2.0
 #define BCXMAX 2.0	/* x interval for boundary condition */
@@ -83,25 +82,25 @@
 #define CIRCLE_PATTERN_B 1  /* pattern of circles for additional particles */
 
 #define ADD_FIXED_OBSTACLES 0   /* set to 1 do add fixed circular obstacles */
-#define OBSTACLE_PATTERN 4  /* pattern of obstacles, see list in global_ljones.c */
+#define OBSTACLE_PATTERN 6  /* pattern of obstacles, see list in global_ljones.c */
 
-#define ADD_FIXED_SEGMENTS 1    /* set to 1 to add fixed segments as obstacles */
-#define SEGMENT_PATTERN 25      /* pattern of repelling segments, see list in global_ljones.c */
+#define ADD_FIXED_SEGMENTS 0    /* set to 1 to add fixed segments as obstacles */
+#define SEGMENT_PATTERN 27      /* 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 TWO_TYPES 1         /* set to 1 to have two types of particles */
-#define TYPE_PROPORTION 0.5 /* proportion of particles of first type */
-#define TWOTYPE_CONFIG 2    /* choice of types, see TTC_ list in global_ljones.c */
+#define TYPE_PROPORTION 1.0 /* 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 */
 #define CENTER_PX 0         /* set to 1 to center horizontal momentum */
 #define CENTER_PY 0         /* set to 1 to center vertical momentum */
 #define CENTER_PANGLE 0     /* set to 1 to center angular momentum */
 
-#define INTERACTION 1       /* particle interaction, see list in global_ljones.c */
-#define INTERACTION_B 1     /* particle interaction for second type of particle, see list in global_ljones.c */
+#define INTERACTION 12       /* 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 SPIN_INTER_FREQUENCY 4.0 /* angular frequency of spin-spin interaction */
 #define SPIN_INTER_FREQUENCY_B 4.0 /* angular frequency of spin-spin interaction for second particle type */
 
@@ -111,9 +110,9 @@
 #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.3	    /* parameter controlling the dimensions of domain */
-#define MU 0.014 	    /* parameter controlling radius of particles */
-#define MU_B 0.01           /* parameter controlling radius of particles of second type */
+#define LAMBDA 0.2	    /* parameter controlling the dimensions of domain */
+#define MU 0.01 	    /* parameter controlling radius of particles */
+#define MU_B 0.012          /* parameter controlling radius of particles of second type */
 #define NPOLY 40            /* number of sides of polygon */
 #define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define AWEDGE 0.5          /* opening angle of wedge, in units of Pi/2 */ 
@@ -122,12 +121,15 @@
 #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 101          /* number of grid point for grid of disks */
-#define NGRIDY 47           /* number of grid point for grid of disks */
+#define NGRIDX 60           /* number of grid point for grid of disks */
+#define NGRIDY 30           /* 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 */
 #define DAM_WIDTH 0.05       /* width of dam for S_DAM segment configuration */
+#define NOBSX 30
+#define NOBSY 20            /* obstacles for O_HEX obstacle pattern */
+#define NTREES 15           /* number of trees in S_TREES */
 
 #define X_SHOOTER -0.2
 #define Y_SHOOTER -0.6
@@ -136,14 +138,14 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 3200      /* number of frames of movie */
-#define NVID 15         /* number of iterations between images displayed on screen */
+#define NSTEPS 2400      /* number of frames of movie */
+#define NVID 30         /* number of iterations between images displayed on screen */
 #define NSEG 25         /* number of segments of boundary of circles */
 #define INITIAL_TIME 0     /* time after which to start saving frames */
 #define OBSTACLE_INITIAL_TIME 150     /* time after which to start moving obstacle */
 #define BOUNDARY_WIDTH 1    /* width of particle boundary */
 #define LINK_WIDTH 2        /* width of links between particles */
-#define CONTAINER_WIDTH 4   /* width of container boundary */
+#define CONTAINER_WIDTH 2   /* width of container boundary */
 
 #define PAUSE 1000         /* number of frames after which to pause */
 #define PSLEEP 1         /* sleep time during pause */
@@ -159,7 +161,13 @@
 /* Plot type, see list in global_ljones.c  */
 
 #define PLOT 5
-#define PLOT_B 11        /* plot type for second movie */
+#define PLOT_B 13        /* plot type for second movie */
+
+/* Background color depending on particle properties */
+
+#define COLOR_BACKGROUND 1  /* 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 COLOR_BONDS 1   /* set to 1 to color bonds according to length */
@@ -176,8 +184,10 @@
 #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_DIRECTION 17  /* Color palette for direction representation */
+#define COLOR_PALETTE_DIRECTION 17   /* Color palette for direction representation */
 #define COLOR_PALETTE_INITIAL_POS 0  /* 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 BLACK 1          /* background */
 
@@ -193,15 +203,17 @@
 #define LUMAMP 0.3       /* amplitude of luminosity variation for scheme C_LUM */
 #define HUEMEAN 220.0    /* mean value of hue for color scheme C_HUE */
 #define HUEAMP -50.0      /* amplitude of variation of hue for color scheme C_HUE */
-#define COLOR_HUESHIFT 0.0     /* shift in color hue (for some cyclic palettes) */
+#define COLOR_HUESHIFT -0.5     /* shift in color hue (for some cyclic palettes) */
 
-#define PRINT_PARAMETERS 0  /* set to 1 to print certain parameters */
+#define PRINT_PARAMETERS 1  /* set to 1 to print certain parameters */
 #define PRINT_TEMPERATURE 0 /* set to 1 to print current temperature */
 #define PRINT_ANGLE 0               /* set to 1 to print obstacle orientation */
 #define PRINT_OMEGA 0               /* set to 1 to print angular speed */
 #define PRINT_PARTICLE_SPEEDS 0     /* set to 1 to print average speeds/momenta of particles */
 #define PRINT_SEGMENTS_SPEEDS 0     /* set to 1 to print velocity of moving segments */
 #define PRINT_SEGMENTS_FORCE 0      /* set to 1 to print force on segments */
+#define PRINT_NPARTICLES 0          /* print number of active particles */
+#define PRINT_TYPE_PROP 0           /* print type proportion */
 #define FORCE_FACTOR 0.1            /* factor controlling length of force vector */
 
 /* particle properties */
@@ -210,41 +222,41 @@
 #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 1000.0        /* energy of particle with hottest color */
-#define HUE_TYPE0 280.0     /* hue of particles of type 0 */
-#define HUE_TYPE1 60.0      /* hue of particles of type 1 */
+#define PARTICLE_EMAX 10000.0        /* energy of particle with hottest color */
+#define HUE_TYPE0 60.0     /* hue of particles of type 0 */
+#define HUE_TYPE1 280.0      /* hue of particles of type 1 */
 #define HUE_TYPE2 140.0      /* hue of particles of type 2 */
 #define HUE_TYPE3 200.0     /* hue of particles of type 3 */
+#define BG_FORCE_SLOPE 7.5e-8   /* contant in BG_FORCE backgound color scheme*/
 
 #define RANDOM_RADIUS 0     /* set to 1 for random circle radius */
 #define DT_PARTICLE 3.0e-6    /* time step for particle displacement */
-#define KREPEL 1.0          /* constant in repelling force between particles */
-#define EQUILIBRIUM_DIST 7.0   /* Lennard-Jones equilibrium distance */
-#define EQUILIBRIUM_DIST_B 6.5  /* Lennard-Jones equilibrium distance for second type of particle */
-#define REPEL_RADIUS 15.0    /* radius in which repelling force acts (in units of particle radius) */
-#define DAMPING 0.02          /* damping coefficient of particles */
-#define OMEGA_INITIAL 10.0        /* initial angular velocity range */
-#define INITIAL_DAMPING 5.0  /* damping coefficient of particles during initial phase */
+#define KREPEL 50.0          /* constant in repelling force between particles */
+#define EQUILIBRIUM_DIST 5.0    /* Lennard-Jones equilibrium distance */
+#define EQUILIBRIUM_DIST_B 5.0  /* Lennard-Jones equilibrium distance for second type of particle */
+#define REPEL_RADIUS 20.0    /* radius in which repelling force acts (in units of particle radius) */
+#define DAMPING 3000.0          /* damping coefficient of particles */
+#define INITIAL_DAMPING 5000.0  /* damping coefficient of particles during initial phase */
 #define DAMPING_ROT 100.0      /* dampint coefficient for rotation of particles */
-#define PARTICLE_MASS 1.0    /* mass of particle of radius MU */
-#define PARTICLE_MASS_B 6.5   /* mass of particle of radius MU */
+#define PARTICLE_MASS 2.0    /* mass of particle of radius MU */
+#define PARTICLE_MASS_B 1.0   /* mass of particle of radius MU */
 #define PARTICLE_INERTIA_MOMENT 0.2     /* moment of inertia of particle */
 #define PARTICLE_INERTIA_MOMENT_B 0.02     /* moment of inertia of second type of particle */
-#define V_INITIAL 1000.0        /* initial velocity range */
+#define V_INITIAL -50.0        /* initial velocity range */
 #define OMEGA_INITIAL 10.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 V_INITIAL_TYPE 1    /* type of initial speed distribution (see VI_ in global_ljones.c) */
+#define V_INITIAL_TYPE 0    /* type of initial speed distribution (see VI_ in global_ljones.c) */
 
-#define THERMOSTAT 0        /* set to 1 to switch on thermostat */
+#define THERMOSTAT 1        /* set to 1 to switch on thermostat */
 #define VARY_THERMOSTAT 0   /* set to 1 for time-dependent thermostat schedule */
 #define SIGMA 5.0           /* noise intensity in thermostat */
-#define BETA 0.0002          /* initial inverse temperature */
+#define BETA 0.00004          /* initial inverse temperature */
 #define MU_XI 0.01           /* friction constant in thermostat */
-#define KSPRING_BOUNDARY 1.0e7    /* confining harmonic potential outside simulation region */
-#define KSPRING_OBSTACLE 1.0e11    /* harmonic potential of obstacles */
-#define NBH_DIST_FACTOR 3.8       /* radius in which to count neighbours */
+#define KSPRING_BOUNDARY 2.0e11    /* confining harmonic potential outside simulation region */
+#define KSPRING_OBSTACLE 2.0e11    /* harmonic potential of obstacles */
+#define NBH_DIST_FACTOR 3.2       /* radius in which to count neighbours */
 #define GRAVITY 0.0            /* gravity acting on all particles */
 #define GRAVITY_X 0.0          /* horizontal gravity acting on all particles */
 #define CIRCULAR_GRAVITY 0     /* set to 1 to have gravity directed to center */
@@ -256,14 +268,19 @@
 #define KSPRING_VICSEK 0.2   /* spring constant for I_VICSEK_SPEED interaction */
 #define VICSEK_REPULSION 10.0    /* repulsion between particles in Vicsek model */
 
-#define ADD_EFIELD 1     /* set to 1 to add an electric field */
-#define EFIELD 10000.0  /* value of electric field */
+#define ADD_EFIELD 0     /* set to 1 to add an electric field */
+#define EFIELD 200000.0   /* value of electric field */
 #define ADD_BFIELD 0     /* set to 1 to add a magnetic field */
-#define BFIELD 1000.0    /* value of magnetic field */
-#define CHARGE 0.0       /* charge of particles of first type */
-#define CHARGE_B 1.0     /* charge of particles of second type */
+#define BFIELD 10000.0       /* value of magnetic field */
+#define CHARGE 1.5      /* charge of particles of first type */
+#define CHARGE_B -1.0     /* charge of particles of second type */
 #define INCREASE_E 0     /* set to 1 to increase electric field */
-#define EFIELD_FACTOR 1000.0    /* factor by which to increase electric field */
+#define EFIELD_FACTOR 1000000.0    /* factor by which to increase electric field */
+#define INCREASE_B 0     /* set to 1 to increase magnetic field */
+#define BFIELD_FACTOR 10000.0    /* factor by which to increase magnetic field */
+#define CHARGE_OBSTACLES 1      /* set to 1 for obstacles to be charged */
+#define OBSTACLE_CHARGE 3.0     /* charge of obstacles */
+#define KCOULOMB_OBSTACLE 1000.0   /* Coulomb force constant for charged obstacles */
 
 #define ROTATION 0           /* set to 1 to include rotation of particles */
 #define COUPLE_ANGLE_TO_THERMOSTAT 1    /* set to 1 to couple angular degrees of freedom to thermostat */
@@ -274,30 +291,30 @@
 #define KTORQUE_DIFF 150.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 0          /* set to 1 to draw cross on particles of second type */
+#define DRAW_CROSS 1          /* set to 1 to draw cross on particles of second type */
 #define SPIN_RANGE 10.0       /* range of spin-spin interaction */
 #define SPIN_RANGE_B 5.0     /* range of spin-spin interaction for second type of particle */
 #define QUADRUPOLE_RATIO 0.6  /* anisotropy in quadrupole potential */ 
 
-#define INCREASE_BETA 0  /* set to 1 to increase BETA during simulation */
-#define BETA_FACTOR 0.02   /* factor by which to change BETA during simulation */
-#define N_TOSCILLATIONS 1.5   /* number of temperature oscillations in BETA schedule */
-#define NO_OSCILLATION 1        /* set to 1 to have exponential BETA change only */
-#define MIDDLE_CONSTANT_PHASE 2000   /* final phase in which temperature is constant */
-#define FINAL_DECREASE_PHASE 1300    /* final phase in which temperature decreases */ 
+#define INCREASE_BETA 1  /* set to 1 to increase BETA during simulation */
+#define BETA_FACTOR 5.0   /* factor by which to change BETA during simulation */
+#define N_TOSCILLATIONS 2.5   /* number of temperature oscillations in BETA schedule */
+#define NO_OSCILLATION 0        /* set to 1 to have exponential BETA change only */
+#define MIDDLE_CONSTANT_PHASE 0   /* final phase in which temperature is constant */
+#define FINAL_DECREASE_PHASE 0    /* final phase in which temperature decreases */ 
 #define FINAL_CONSTANT_PHASE -1     /* final phase in which temperature is constant */
 
 #define DECREASE_CONTAINER_SIZE 0   /* set to 1 to decrease size of container */
 #define SYMMETRIC_DECREASE 0        /* set tp 1 to decrease container symmetrically */
 #define COMPRESSION_RATIO 0.3       /* final size of container */
 #define RESTORE_CONTAINER_SIZE 1    /* set to 1 to restore container to initial size at end of simulation */
-#define RESTORE_TIME 700            /* time before end of sim at which to restore size */
+#define RESTORE_TIME 1200            /* time before end of sim at which to restore size */
 
 #define MOVE_OBSTACLE 0     /* set to 1 to have a moving obstacle */
 #define CENTER_VIEW_ON_OBSTACLE 0   /* set to 1 to center display on moving obstacle */
 #define RESAMPLE_Y 0         /* set to 1 to resample y coordinate of moved particles (for shock waves) */
 #define NTRIALS 2000         /* number of trials when resampling */
-#define OBSTACLE_RADIUS 0.15  /* radius of obstacle for circle boundary conditions */
+#define OBSTACLE_RADIUS 0.015  /* radius of obstacle for circle boundary conditions */
 #define FUNNEL_WIDTH  0.25   /* funnel width for funnel boundary conditions */
 #define OBSTACLE_XMIN 0.0    /* initial position of obstacle */
 #define OBSTACLE_XMAX 3.0    /* final position of obstacle */
@@ -321,8 +338,8 @@
 
 #define ADD_PARTICLES 0   /* set to 1 to add particles */
 #define ADD_TIME 0       /* time at which to add first particle */
-#define ADD_PERIOD 10       /* time interval between adding further particles */
-#define N_ADD_PARTICLES 2  /* number of particles to add */
+#define ADD_PERIOD 4       /* time interval between adding further particles */
+#define N_ADD_PARTICLES 3  /* number of particles to add */
 #define FINAL_NOADD_PERIOD 0  /* final period where no particles are added */
 #define SAFETY_FACTOR 2.0  /* no particles are added at distance less than MU*SAFETY_FACTOR of other particles */
 
@@ -346,12 +363,12 @@
 #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 RELEASE_ROCKET_AT_DEACTIVATION 0    /* set to 1 to limit segments velocity before segment release */
-#define SEGMENTS_X0 0.0        /* initial position of segments */
-#define SEGMENTS_Y0 -0.6      /* initial position of segments */
+#define SEGMENTS_X0 1.5        /* initial position of segments */
+#define SEGMENTS_Y0 0.0        /* initial position of segments */
 #define SEGMENTS_VX0 0.0       /* initial velocity of segments */
 #define SEGMENTS_VY0 0.0      /* initial velocity of segments */
 #define DAMP_SEGS_AT_NEGATIVE_Y 0   /* set to 1 to dampen segments when y coordinate is negative */
-#define SHOW_SEGMENTS_PRESSURE 0    /* set to 1 to show (averaged) pressure acting on segments */
+#define SHOW_SEGMENTS_PRESSURE 1    /* set to 1 to show (averaged) pressure acting on segments */
 #define SEGMENT_PMAX 7.5e7        /* pressure of segment with hottest color */
 #define P_AVRG_FACTOR 0.02      /* factor in computation of mean pressure */
 
@@ -406,6 +423,16 @@
 #define WALL_VMAX 100.0     /* max speed of wall */
 #define WALL_TIME 0         /* time during which to keep wall */
 
+#define CHANGE_TYPES 0      /* set to 1 to change type proportion in course of simulation */
+#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 KSPRING_PAIRS 1.0e10    /* spring constant for pair interaction */
+#define NPARTNERS 2         /* number of partners of particles */
+#define PAIRING_TYPE 1      /* type of pairing, see POLY_ in global_ljones.c */
+#define PARTNER_ANGLE 104.45    /* angle (in degrees) between anions for POLY_WATER case */
+
 #define NXMAZE 12      /* width of maze */
 #define NYMAZE 12      /* height of maze */
 #define MAZE_MAX_NGBH 4     /* max number of neighbours of maze cell */
@@ -418,8 +445,8 @@
 #define FLOOR_OMEGA 0      /* set to 1 to limit particle momentum to PMAX */
 #define PMAX 1000.0        /* maximal force */
 
-#define HASHX 60    /* size of hashgrid in x direction */
-#define HASHY 30    /* size of hashgrid in y direction */
+#define HASHX 100   /* size of hashgrid in x direction */
+#define HASHY 50    /* 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 */
 
@@ -457,8 +484,8 @@ double sinangle[NSEG+1];      /* precomputed trig functions of angles to draw ci
 
 #include "global_ljones.c"
 #include "sub_maze.c"
-#include "sub_lj.c"
 #include "sub_hashgrid.c"
+#include "sub_lj.c"
 
 FILE *lj_time_series, *lj_final_position;
 
@@ -501,6 +528,23 @@ double efield_schedule(int i)
 }
 
 
+double bfield_schedule(int i)
+{
+    static double bfactor;
+    static int first = 1;
+    double bfield;
+    
+    if (first) 
+    {
+        bfactor = BFIELD_FACTOR/(double)(INITIAL_TIME + NSTEPS);
+        first = 0;
+    }
+    bfield = BFIELD*(double)i*bfactor;
+    printf("B = %.3lg\n", bfield);
+    return(bfield);
+}
+
+
 double temperature_schedule(int i)
 {
     static double bexponent, omega, bexp2;
@@ -728,11 +772,11 @@ double radius_schedule(int i)
 double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], 
                         double qx[NMAXCIRCLES], double qy[NMAXCIRCLES], double qangle[NMAXCIRCLES],
                         double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], 
-                        double beta, int *nactive, int *nsuccess, int *nmove, int initial_phase)
+                        double beta, int *nactive, int *nsuccess, int *nmove, int *ncoupled, int initial_phase)
 {
     double a, totalenergy = 0.0, damping, direction, dmean;  
     static double b = 0.25*SIGMA*SIGMA*DT_PARTICLE/MU_XI, xi = 0.0;
-    int j, move;
+    int j, move, ncoup;
     
     if (initial_phase) damping = INITIAL_DAMPING;
     else damping = DAMPING;
@@ -782,16 +826,20 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
                 
     /* compute kinetic energy */
 //     *nactive = 0;
+    ncoup = 1;
     for (j=0; j<ncircles; j++)
         if ((particle[j].active)&&(particle[j].thermostat))
         {
             totalenergy += particle[j].energy;
+            ncoup++;
 //             *nactive++;
         }
     totalenergy *= DIMENSION_FACTOR;    /* normalize energy to take number of degrees of freedom into account */
     if (THERMOSTAT_ON)
     {
-        a = DT_PARTICLE*(totalenergy - (double)*nactive/beta)/MU_XI;
+        /* 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);
     }
@@ -860,6 +908,7 @@ double evolve_particles(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HA
 //                 }
     }
     
+    *ncoupled = ncoup;
     return(totalenergy);
 }
 
@@ -1164,12 +1213,12 @@ void evolve_segment_groups(t_segment segment[NMAXSEGMENTS], int time, t_group_se
 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, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio;
+            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, ymin, ymax, delta_energy, speed, ratio = 1.0, ratioc, cum_etot = 0.0, emean = 0.0, radius_ratio, t;
     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, 
         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;
+        group, gshift, n_total_active = 0, ncollisions = 0, ncoupled = 1;
     int *particle_numbers;
     static int imin, imax;
     static short int first = 1;
@@ -1261,6 +1310,8 @@ void animation()
     printf("Initializing configuration\n");
 
     params.nactive = initialize_configuration(particle, hashgrid, obstacle, px, py, pangle, tracer_n, segment);
+    
+    printf("%i active particles\n", params.nactive);
      
 //     xi = 0.0;
     
@@ -1288,6 +1339,7 @@ void animation()
         if (INCREASE_KREPEL) params.krepel = repel_schedule(i);
         if (INCREASE_BETA) params.beta = temperature_schedule(i);  
         if (INCREASE_E) params.efield = efield_schedule(i);
+        if (INCREASE_B) params.bfield = bfield_schedule(i);
         if (DECREASE_CONTAINER_SIZE) 
         {
             params.xmincontainer = container_size_schedule(i);
@@ -1368,7 +1420,7 @@ void animation()
                 compute_particle_force(j, params.krepel, particle, hashgrid);
                 
                 /* take care of boundary conditions */
-                params.fboundary += compute_boundary_force(j, particle, obstacle, segment, params.xmincontainer, params.xmaxcontainer, &pleft, &pright, pressure, wall);
+                params.fboundary += compute_boundary_force(j, particle, obstacle, segment, params.xmincontainer, params.xmaxcontainer, &pleft, &pright, pressure, wall, params.krepel);
 
                 /* align velocities in case of Vicsek models */
 //                 if (VICSEK_INT) 
@@ -1414,10 +1466,8 @@ void animation()
                 /* add magnetic force */
                 if (ADD_BFIELD)
                 {
-//                     particle[j].fx += BFIELD*particle[j].charge*particle[j].vy;
-//                     particle[j].fy -= BFIELD*particle[j].charge*particle[j].vx;
-                    particle[j].fx += BFIELD*particle[j].charge*particle[j].vy*particle[j].mass_inv;
-                    particle[j].fy -= BFIELD*particle[j].charge*particle[j].vx*particle[j].mass_inv;
+                    particle[j].fx += params.bfield*particle[j].charge*particle[j].vy*particle[j].mass_inv;
+                    particle[j].fy -= params.bfield*particle[j].charge*particle[j].vx*particle[j].mass_inv;
                 }
                 
                 if (FLOOR_FORCE)
@@ -1432,7 +1482,7 @@ void animation()
             }
             
             /* timestep of thermostat algorithm */
-            totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, params.beta, &params.nactive, &nsuccess, &nmove, i < INITIAL_TIME);
+            totalenergy = evolve_particles(particle, hashgrid, qx, qy, qangle, px, py, pangle, params.beta, &params.nactive, &nsuccess, &nmove, &ncoupled, i < INITIAL_TIME);
             
             
             /* evolution of lid coordinate */
@@ -1529,6 +1579,7 @@ void animation()
         }
                 
         printf("Mean kinetic energy: %.3f\n", totalenergy/(double)ncircles); 
+        printf("Kinetic energy by coupled particle: %.3f\n",  totalenergy/(double)ncoupled);
         
         if ((!THERMOSTAT)&&(LIMIT_ENERGY))
         {
@@ -1582,9 +1633,17 @@ void animation()
             }
 //             draw_collisions(collisions, ncollisions);
         }
+        
+        /* case of varying type proportion */
+        if (CHANGE_TYPES) 
+        {
+            t = (double)i/(double)NSTEPS;
+            params.prop = PROP_MIN*(1.0-t) + PROP_MAX*t;
+            change_type_proportion(particle, params.prop);
+        }
  
         if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-        draw_particles(particle, PLOT, params.beta, collisions, ncollisions);
+        draw_particles(particle, PLOT, params.beta, collisions, ncollisions, BG_COLOR, hashgrid);
         draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
 
         /* add a particle */
@@ -1592,6 +1651,9 @@ void animation()
         {
             for (k=0; k<N_ADD_PARTICLES; k++)
             nadd_particle = add_particles(particle, px, py, nadd_particle);
+//             params.nactive = nadd_particle;
+            params.nactive = 0;
+            for (j=0; j<ncircles; j++) if (particle[j].active) params.nactive++;
         }
         
         /* change particle radius */
@@ -1670,7 +1732,7 @@ void animation()
             if ((i >= INITIAL_TIME)&&(DOUBLE_MOVIE))
             {
                 if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-                draw_particles(particle, PLOT_B, params.beta, collisions, ncollisions);
+                draw_particles(particle, PLOT_B, params.beta, collisions, ncollisions, BG_COLOR_B, hashgrid);
                 draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
                 if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 0);
                 if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i);
@@ -1722,7 +1784,7 @@ void animation()
         {
             blank();
             if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-            draw_particles(particle, PLOT, params.beta, collisions, ncollisions); 
+            draw_particles(particle, PLOT, params.beta, collisions, ncollisions, BG_COLOR, hashgrid); 
             draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
             if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 1);
             if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i);
@@ -1748,7 +1810,7 @@ void animation()
         if (DOUBLE_MOVIE) 
         {
             if (TRACER_PARTICLE) draw_trajectory(trajectory, traj_position, traj_length);
-            draw_particles(particle, PLOT_B, params.beta, collisions, ncollisions);
+            draw_particles(particle, PLOT_B, params.beta, collisions, ncollisions, BG_COLOR_B, hashgrid);
             draw_container(params.xmincontainer, params.xmaxcontainer, obstacle, segment, wall);
             if (PRINT_PARAMETERS) print_parameters(params, PRINT_LEFT, pressure, 0);
             if (PLOT_SPEEDS) draw_speed_plot(group_speeds, i);

mangrove.c

@@ -256,6 +256,10 @@
 #define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 #define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 #define OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
+#define MU_B 1.0           /* parameter controlling the dimensions of domain */
+#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 DRAW_WAVE_TIMESERIES 0  /* set to 1 to draw a time series of the wave */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"

rde.c

@@ -396,6 +396,11 @@
 #define WAVE_PACKET_RADIUS 20            /* radius of wave packets */
 #define OSCIL_LEFT_YSHIFT 25.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 #define DRAW_WAVE_PROFILE 1     /* set to 1 to draw a profile of the wave */
+#define MU_B 1.0           /* parameter controlling the dimensions of domain */
+#define GAMMA 0.0          /* damping factor in wave equation */
+#define GAMMAB 0.0         /* damping factor in wave equation */
+#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 */
 /* end of constants added only for compatibility with wave_common.c */
 
 

schrodinger.c

@@ -193,6 +193,11 @@
 #define N_WAVE_PACKETS 15               /* number of wave packets */
 #define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 #define OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
+#define MU_B 1.0           /* parameter controlling the dimensions of domain */
+#define GAMMA 0.0          /* damping factor in wave equation */
+#define GAMMAB 0.0         /* damping factor in wave equation */
+#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 */
 /* end of constants only used by sub_wave and sub_maze */
 
 

sub_hashgrid.c

@@ -2,6 +2,14 @@
 /* Note: a possible improvement would be to use pointers instead of integers */
 /* when referencing other hashgrid cells or particles */ 
 
+ double module2(double x, double y)   /* Euclidean norm */
+ {
+	double m;
+
+	m = sqrt(x*x + y*y);
+	return(m);
+ }
+
 int bc_grouped(int bc)
 /* regroup boundary conditions by type: rectangular = 0, periodic = 1, Klein = 2, Boy = 3, L shape = 4, ... */
 {
@@ -70,6 +78,41 @@ int hash_cell(double x, double y)
 }
 
 
+void init_hashcell_coordinates(t_hashgrid hashgrid[HASHX*HASHY])
+/* initialise coordinates of corners of hashcells, needed for option COLOR_BACKGROUND */
+{
+    static int first = 1;
+    static double lx, ly, padding, dx, dy;
+    int i, j, n;
+    double x, y;
+    
+    if (first)
+    {
+        if (!NO_WRAP_BC) padding = 0.0;
+        else padding = HASHGRID_PADDING;
+        lx = BCXMAX - BCXMIN + 2.0*padding;
+        ly = BCYMAX - (BCYMIN) + 2.0*padding;
+        dx = 1.0*lx/(double)HASHX;
+        dy = 1.0*ly/(double)HASHY;
+        first = 0;
+    }
+    
+    for (i=0; i<HASHX; i++)
+        for (j=0; j<HASHY; j++)
+        {
+            x = BCXMIN - padding + (double)i*dx;
+            y = BCYMIN - padding + (double)j*dy;
+            n = mhash(i, j);
+            hashgrid[n].x1 = x;
+            hashgrid[n].y1 = y;
+            hashgrid[n].x2 = x + dx;
+            hashgrid[n].y2 = y + dy;
+            hashgrid[n].charge = 0.0;
+        }
+    
+}
+
+
 void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
 /* initialise table of neighbouring cells for each hashgrid cell, depending on boundary condition */
 {
@@ -368,7 +411,10 @@ void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
 //         }
 // //         sleep(1);
 //     }
-    sleep(1);
+
+    if (COLOR_BACKGROUND) init_hashcell_coordinates(hashgrid);
+    
+//     sleep(1);
 }
 
 void update_hashgrid(t_particle* particle, t_hashgrid* hashgrid, int verbose)

sub_lj.c

@@ -157,14 +157,6 @@ void write_text( double x, double y, char *st)
 	return(res);
  }
 
- double module2(double x, double y)   /* Euclidean norm */
- {
-	double m;
-
-	m = sqrt(x*x + y*y);
-	return(m);
- }
-
  double argument(double x, double y)
  {
 	double alph;
@@ -1137,7 +1129,7 @@ void add_obstacle(double x, double y, double radius, t_obstacle obstacle[NMAXOBS
 void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES])
 /* initialise circular obstacle configuration */
 {
-    int i, j, n; 
+    int i, j, n, jmin, jmax; 
     double x, y, dx, dy, width, lpocket, xmid = 0.5*(BCXMIN + BCXMAX), radius;
     
     switch (OBSTACLE_PATTERN) {
@@ -1262,12 +1254,48 @@ void init_obstacle_config(t_obstacle obstacle[NMAXOBSTACLES])
             nobstacles = 4;
             break;
         }
-        
+        case (O_HEX):
+        {
+            dx = (XMAX - XMIN)/(double)NOBSX;
+            dy = (YMAX - YMIN)/(double)NOBSY;
+            n = 0;
+            
+            if ((ADD_FIXED_SEGMENTS)&&(SEGMENT_PATTERN == S_CYLINDER))
+            /* pseudo-cylindrical boundary conditions (e.g. for Hall effect) */
+            {
+                jmin = 1; 
+                jmax = NOBSY-1;
+            }
+            else
+            {
+                jmin = 0;
+                jmax = NOBSY;
+            }
+            
+            for (i=0; i<NOBSX; i++)
+                for (j=jmin; j<jmax; j++)
+                {
+                    obstacle[n].xc = XMIN + ((double)i + 0.25)*dx;
+                    obstacle[n].yc = YMIN + ((double)j + 0.5)*dy;
+                    if (j%2 == 1) obstacle[n].xc += 0.5*dx;
+                    if (obstacle[n].xc > XMAX) obstacle[n].xc += (XMAX - XMIN);
+                    obstacle[n].radius = OBSTACLE_RADIUS;
+                    obstacle[n].active = 1;
+                    n++;
+                }
+            nobstacles = n;
+//             printf("Added %i obstacles\n", nobstacles);
+//             for (n=0; n<nobstacles; n++) printf("Obstacle %i at (%.3f, %.3f)\n", n, obstacle[n].xc, obstacle[n].yc);
+            break;
+        }
         default: 
         {
             printf("Function init_obstacle_config not defined for this pattern \n");
         }
     }
+    
+    if (CHARGE_OBSTACLES)
+        for (n=0; n<nobstacles; n++) obstacle[n].charge = OBSTACLE_CHARGE;
 }
 
 void add_rotated_angle_to_segments(double x1, double y1, double x2, double y2, double width, int center, t_segment segment[NMAXSEGMENTS], int group)
@@ -1406,7 +1434,154 @@ void add_circle_to_segments(double x, double y, double r, int nsegs, double angl
 }
         
         
-
+void add_tree_to_segments(double x, double y, double r, double lfoot[2], double rfoot[2], t_segment segment[NMAXSEGMENTS], int group)
+/* add segments forming a tree to linear obstacle configuration */
+{
+    int i, n = nsegments, nplus, nminus, nadd; 
+    double angle, h, h1, h2, h3, x1, x2, xp1, xp2, dx1, dx2, xt, yt;
+    
+    angle = PI/3.0;
+    h = r*tan(angle);
+    h1 = 0.56*h;
+    h2 = h1;
+    dx1 = h1/tan(angle);
+    dx2 = 0.5*dx1;
+    x1 = r - dx1;
+    xp1 = x1 + dx2;
+    x2 = xp1 - dx1;
+    xp2 = x2 + dx2;
+    h3 = xp2*tan(angle);
+    xt = 0.5*dx2;
+    yt = dx1;
+    
+    nadd = 15;
+    
+    if (nsegments + nadd < NMAXSEGMENTS) 
+    {
+        /* bottom left */
+        segment[n].x1 = x - r;
+        segment[n].y1 = y;
+        segment[n].angle1 = PID + angle;
+        segment[n].angle2 = PI + PID;
+        segment[n].concave = 1;
+        n++;
+        
+        /* level 1 */
+        segment[n].x1 = x - x1;
+        segment[n].y1 = y + h1;
+        segment[n].concave = 0;
+        n++;
+        
+        segment[n].x1 = x - xp1;
+        segment[n].y1 = y + h1;
+        segment[n].angle1 = PID + angle;
+        segment[n].angle2 = PI + PID;
+        segment[n].concave = 1;
+        n++;
+        
+        /* level 2 */
+        segment[n].x1 = x - x2;
+        segment[n].y1 = y + h1 + h2;
+        segment[n].concave = 0;
+        n++;
+        
+        segment[n].x1 = x - xp2;
+        segment[n].y1 = y + h1 + h2;
+        segment[n].angle1 = PID + angle;
+        segment[n].angle2 = PI + PID;
+        segment[n].concave = 1;
+        n++;
+        
+        /* top */
+        segment[n].x1 = x;
+        segment[n].y1 = y + h1 + h2 + h3;
+        segment[n].angle1 = PID - angle;
+        segment[n].angle2 = PID + angle;
+        segment[n].concave = 1;
+        n++;
+        
+        /* level 2 */
+        segment[n].x1 = x + xp2;
+        segment[n].y1 = y + h1 + h2;
+        segment[n].angle1 = -PID;
+        segment[n].angle2 = angle;
+        segment[n].concave = 1;
+        n++;
+        
+        segment[n].x1 = x + x2;
+        segment[n].y1 = y + h1 + h2;
+        segment[n].concave = 0;
+        n++;
+        
+        /* level 1 */
+        segment[n].x1 = x + xp1;
+        segment[n].y1 = y + h1;
+        segment[n].angle1 = -PID;
+        segment[n].angle2 = angle;
+        segment[n].concave = 1;
+        n++;
+        
+        segment[n].x1 = x + x1;
+        segment[n].y1 = y + h1;
+        segment[n].concave = 0;
+        n++;
+        
+        /* bottom right */
+        segment[n].x1 = x + r;
+        segment[n].y1 = y;
+        segment[n].angle1 = -PID;
+        segment[n].angle2 = angle;
+        segment[n].concave = 1;
+        n++;
+        
+        /*stem */
+        segment[n].x1 = x + xt;
+        segment[n].y1 = y;
+        segment[n].concave = 0;
+        n++;
+        
+        segment[n].x1 = x + xt;
+        segment[n].y1 = y - yt;
+        segment[n].angle1 = -PID;
+        segment[n].angle2 = 0.0;
+        segment[n].concave = 1;
+        n++;
+        
+        segment[n].x1 = x - xt;
+        segment[n].y1 = y - yt;
+        segment[n].angle1 = PI;
+        segment[n].angle2 = 3.0*PID;
+        segment[n].concave = 1;
+        n++;
+        
+        segment[n].x1 = x - xt;
+        segment[n].y1 = y;
+        segment[n].concave = 0;
+        n++;
+        
+        rfoot[0] = x + xt;
+        rfoot[1] = y - yt;
+        lfoot[0] = x - xt;
+        lfoot[1] = y - yt;
+        
+        for (i=0; i<nadd-1; i++)
+        {
+            segment[nsegments+i].x2 = segment[nsegments+i+1].x1;
+            segment[nsegments+i].y2 = segment[nsegments+i+1].y1;
+        }
+        segment[nsegments+nadd-1].x2 = segment[nsegments].x1;
+        segment[nsegments+nadd-1].y2 = segment[nsegments].y1;
+        
+        for (i=0; i<nadd; i++)
+        {
+//             segment[nsegments+i].concave = 1;
+            segment[nsegments+i].group = group;
+        }
+    
+        nsegments += nadd;
+    }
+    else printf("Warning: NMAXSEGMENTS too small\n");
+}
 
 double nozzle_width(double x, double width, int nozzle_shape)
 /* width of bell-shaped nozzle */
@@ -1620,6 +1795,7 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
 {
     int i, j, cycle = 0, iminus, iplus, nsides, n, concave = 1; 
     double angle, angle2, dangle, dx, width, height, a, b, length, xmid = 0.5*(BCXMIN + BCXMAX), lpocket, r, x, x1, y1, x2, y2, nozx, nozy, y, yy, dy, ca, sa, padding;
+    double lfoot[NTREES][2], rfoot[NTREES][2];
     
     switch (SEGMENT_PATTERN) {
         case (S_RECTANGLE):
@@ -2478,6 +2654,51 @@ void init_segment_config(t_segment segment[NMAXSEGMENTS])
             
             break;
         }
+        case (S_CYLINDER):
+        {
+            add_rectangle_to_segments(XMAX + 0.1, YMAX - 0.05, XMIN - 0.1, YMAX + 1.0, segment, 0);
+            add_rectangle_to_segments(XMAX + 0.1, YMIN - 1.0, XMIN - 0.1, YMIN + 0.05, segment, 0);
+            
+            cycle = 0;
+            concave = 1;    
+            
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].inactivate = 0;
+            }
+            
+            break;
+        }
+        case (S_TREE):
+        {
+            for (i=0; i<NTREES; i++)
+            {
+                x = XMIN + ((double)i - 0.5 + 0.1*(double)rand()/(double)RAND_MAX)*(XMAX-XMIN)*1.1/(double)NTREES; 
+                y = YMIN + 0.1 + 0.014*(double)i*(YMAX - YMIN);
+                add_tree_to_segments(x, y, LAMBDA*(1.0 +  0.4*(double)rand()/(double)RAND_MAX), lfoot[i], rfoot[i], segment, 0);
+            }
+            /* add ground */
+            for (i=0; i<NTREES-1; i++)
+            {
+                segment[nsegments].x1 = rfoot[i][0];
+                segment[nsegments].y1 = rfoot[i][1];
+                segment[nsegments].x2 = lfoot[i+1][0];
+                segment[nsegments].y2 = lfoot[i+1][1];
+                segment[nsegments].concave = 0;
+                nsegments++;
+            }
+            cycle = 0;
+            concave = 0;
+            ngroups = 1;
+            
+            for (i=0; i<nsegments; i++) 
+            {
+                segment[i].concave = 1;
+                segment[i].inactivate = 0;
+            }
+            
+            break;
+        }
         default: 
         {
             printf("Function init_segment_config not defined for this pattern \n");
@@ -2844,6 +3065,12 @@ int in_segment_region(double x, double y, t_segment segment[NMAXSEGMENTS])
             if (x > x1 + length + 0.1) return(1);
             return(vabs(y - SEGMENTS_Y0 + 0.2*cos(DPI*(x-x1)/length)) > 0.05);
         }
+        case (S_CYLINDER):
+        {
+            if (y > YMAX - 0.05 - MU) return(0);
+            if (y < YMIN + 0.05 + MU) return(0);
+            return(1);
+        }
         case (S_MAZE):
         {
             for (i=0; i<nsegments; i++)
@@ -3872,6 +4099,7 @@ int add_particle(double x, double y, double vx, double vy, double mass, short in
         particle[i].neighb = 0;
         particle[i].diff_neighb = 0;
         particle[i].thermostat = 1;
+        particle[i].charge = CHARGE;
 
         particle[i].energy = 0.0;
         particle[i].emean = 0.0;
@@ -4057,7 +4285,8 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl
         {
             hue = (double)(particle.diff_neighb+1)/(double)(particle.neighb+1);
 //                hue = PARTICLE_HUE_MIN + (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*hue;
-            hue = 180.0*(1.0 + hue);
+//             hue = 180.0*(1.0 + hue);
+            hue = 20.0 + 320.0*hue;
             *radius = particle.radius;
             *width = BOUNDARY_WIDTH;
             break;
@@ -4166,9 +4395,12 @@ void compute_particle_colors(t_particle particle, int plot, double rgb[3], doubl
         }
         case (P_DIFF_NEIGHB):  
         {
-            hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb);
-            hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx);
-            hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby);
+            hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_DIFFNEIGH);
+            hsl_to_rgb_palette(hue, 0.9, 0.5, rgbx, COLOR_PALETTE_DIFFNEIGH);
+            hsl_to_rgb_palette(hue, 0.9, 0.5, rgby, COLOR_PALETTE_DIFFNEIGH);
+//             hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb);
+//             hsl_to_rgb_twilight(hue, 0.9, 0.5, rgbx);
+//             hsl_to_rgb_twilight(hue, 0.9, 0.5, rgby);
             break;
         }
         case (P_INITIAL_POS): 
@@ -4241,7 +4473,9 @@ void set_segment_pressure_color(double pressure, double lum, double rgb[3])
     if (hue > PARTICLE_HUE_MIN) hue = PARTICLE_HUE_MIN;
     else if (hue < PARTICLE_HUE_MAX) hue = PARTICLE_HUE_MAX;
         
-    hsl_to_rgb_turbo(hue, 0.9, lum, rgb);
+//     hsl_to_rgb_turbo(hue, 0.9, lum, rgb);
+    
+    hsl_to_rgb_palette(hue, 0.9, lum, rgb, COLOR_PALETTE_PRESSURE);
     
     glColor3f(rgb[0], rgb[1], rgb[2]);
 }
@@ -4781,7 +5015,8 @@ void draw_one_particle(t_particle particle, double xc, double yc, double radius,
     
     /* draw crosses on particles of second type */
     if ((TWO_TYPES)&&(DRAW_CROSS))
-        if (particle.type == 2)
+//         if (particle.type == 2)
+        if (particle.charge > 0.0)
         {
             if (ROTATION) 
             {
@@ -4922,9 +5157,158 @@ void draw_trajectory(t_tracer trajectory[TRAJECTORY_LENGTH*N_TRACER_PARTICLES],
         }
     }
 }
+
+
+void color_background(t_particle particle[NMAXCIRCLES], int bg_color, t_hashgrid hashgrid[HASHX*HASHY])
+/* color background according to particle properties */
+{
+    int i, j, k, n, p, q, m, nnb, number, avrg_fact;
+    double rgb[3], hue, value, p1, p2, pp1, pp2, oldhue;
+    static int first = 1, counter = 0;
+    
+    p1 = 0.75;
+    p2 = 1.0 - p1;
+    pp1 = 0.95;
+    pp2 = 1.0 - pp1;
+   
+    glBegin(GL_QUADS);
+    for (i=0; i<HASHX; i++)
+        for (j=0; j<HASHY; j++)
+        {
+            n = mhash(i, j);
+            if (first) 
+            {
+                hashgrid[n].hue1 = 180.0;
+                hashgrid[n].hue2 = 180.0;
+            }
+            /* set two old values for option DOUBLE_MOVIE */
+            if (DOUBLE_MOVIE)
+            {
+                if (counter) oldhue = hashgrid[n].hue1;
+                else oldhue = hashgrid[n].hue2;
+            }
+            else oldhue = hashgrid[n].hue1;
+            
+            switch (bg_color) {
+                case (BG_DENSITY):
+                {
+                    nnb = hashgrid[n].nneighb;
+                    number = 0;
+                    for (q=0; q<nnb; q++)
+                    {
+                        m = hashgrid[n].neighbour[q];
+                        number += hashgrid[m].number;
+                    }
+                    number += hashgrid[n].number;
+                    
+//                     hue = 50.0*(double)hashgrid[n].number;
+                    hue = 75.0*(double)number/(double)(nnb + 1);
+                    hue = p1*oldhue + p2*hue;
+                    rgb[0] = hue/360.0;
+                    rgb[1] = hue/360.0;
+                    rgb[2] = hue/360.0;
+                    break;
+                }
+                case (BG_CHARGE):
+                {
+                    avrg_fact = 3;      /* weight of central cell in hashgrid average */
+                    value = 0.0;
+                    nnb = hashgrid[n].nneighb;
+                    for (q=0; q<nnb; q++)
+                    {
+                        m = hashgrid[n].neighbour[q];
+                        for (k=0; k<hashgrid[m].number; k++)
+                        {
+                            p = hashgrid[m].particles[k];
+                            value += particle[p].charge;
+                        }
+                    }
+                    /* hashcell n counts double */
+                    for (k=0; k<hashgrid[n].number; k++)
+                    {
+                        p = hashgrid[n].particles[k];
+                        value += (double)(avrg_fact-1)*particle[p].charge;
+                    }
+                    value *= 1.0/(double)(nnb + avrg_fact);
+                    if (CHARGE_OBSTACLES) value += hashgrid[n].charge;
+                    hue = (-tanh(0.5*value)+1.0)*180.0;
+                    hue = p1*oldhue + p2*hue;
+                    hsl_to_rgb_twilight(hue, 0.9, 0.5, rgb);
+                    break;
+                }
+                case (BG_EKIN):
+                {
+                    value = 0.0;
+                    nnb = hashgrid[n].nneighb;
+                    for (q=0; q<nnb; q++)
+                    {
+                        m = hashgrid[n].neighbour[q];
+                        for (k=0; k<hashgrid[m].number; k++)
+                        {
+                            p = hashgrid[m].particles[k];
+                            value += particle[p].energy;
+                        }
+                    }
+                    /* hashcell n counts double */
+                    for (k=0; k<hashgrid[n].number; k++)
+                    {
+                        p = hashgrid[n].particles[k];
+                        value += particle[p].energy;
+                    }
+                    value *= 1.0/(double)(nnb + 1);
+                    hue = ENERGY_HUE_MIN + (ENERGY_HUE_MAX - ENERGY_HUE_MIN)*value/PARTICLE_EMAX;
+                    if (hue > ENERGY_HUE_MIN) hue = ENERGY_HUE_MIN;
+                    if (hue < ENERGY_HUE_MAX) hue = ENERGY_HUE_MAX;
+                    hue = p1*oldhue + p2*hue;
+                    hsl_to_rgb_palette(hue, 0.9, 0.5, rgb, COLOR_PALETTE_EKIN);
+                    break;
+                }
+                case (BG_FORCE):
+                {
+                    value = 0.0;
+                    nnb = hashgrid[n].nneighb;
+                    for (q=0; q<nnb; q++)
+                    {
+                        m = hashgrid[n].neighbour[q];
+                        for (k=0; k<hashgrid[m].number; k++)
+                        {
+                            p = hashgrid[m].particles[k];
+                            value += module2(particle[p].fx, particle[p].fy);
+                        }
+                    }
+                    /* hashcell n counts double */
+                    for (k=0; k<hashgrid[n].number; k++)
+                    {
+                        p = hashgrid[n].particles[k];
+                        value += module2(particle[p].fx, particle[p].fy);
+                    }
+                    value *= BG_FORCE_SLOPE/(double)(nnb + 1);
+                    hue = (1.0 - tanh(value))*360.0;
+                    hue = pp1*oldhue + pp2*hue;
+                    hsl_to_rgb_turbo(hue, 0.9, 0.5, rgb);
+                    break;
+                }
+            }
+            if (DOUBLE_MOVIE)
+            {
+                if (counter) hashgrid[n].hue1 = hue;
+                else hashgrid[n].hue2 = hue;
+            }
+            else hashgrid[n].hue1 = hue;
+            
+            glColor3f(rgb[0], rgb[1], rgb[2]);
+            glVertex2d(hashgrid[n].x1, hashgrid[n].y1);
+            glVertex2d(hashgrid[n].x2, hashgrid[n].y1);
+            glVertex2d(hashgrid[n].x2, hashgrid[n].y2);
+            glVertex2d(hashgrid[n].x1, hashgrid[n].y2);
+        }
+    glEnd();
+    first = 0;
+    counter = 1 - counter;
+}
     
 
-void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_collision *collisions, int ncollisions)
+void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_collision *collisions, int ncollisions, int bg_color, t_hashgrid hashgrid[HASHX*HASHY])
 {
     int i, j, k, m, width, nnbg, nsides;
     double ej, hue, huex, huey, rgb[3], rgbx[3], rgby[3], radius, x1, y1, x2, y2, angle, ca, sa, length, linkcolor, sign = 1.0, angle1, signy = 1.0, periodx, periody, x, y, lum, logratio;
@@ -4932,22 +5316,33 @@ void draw_particles(t_particle particle[NMAXCIRCLES], int plot, double beta, t_c
     
 //     if (!TRACER_PARTICLE) blank();
     if (plot == P_NOPARTICLE) blank();
+    
+    if ((COLOR_BACKGROUND)&&(bg_color > 0)) color_background(particle, bg_color, hashgrid);
+    else blank();
+    
     glColor3f(1.0, 1.0, 1.0);
     
     /* show region of partial thermostat */
-    if ((PARTIAL_THERMO_COUPLING)&&(PARTIAL_THERMO_REGION == TH_INBOX))
+    if (PARTIAL_THERMO_COUPLING)
     {
-        if (INCREASE_BETA)
+        if (PARTIAL_THERMO_REGION == TH_INBOX) 
         {
-            logratio = log(beta/BETA)/log(0.5*BETA_FACTOR);
-            if (logratio > 1.0) logratio = 1.0;
-            else if (logratio < 0.0) logratio = 0.0;
-            if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio;
-            else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio;
+            if (INCREASE_BETA)
+            {
+                logratio = log(beta/BETA)/log(0.5*BETA_FACTOR);
+                if (logratio > 1.0) logratio = 1.0;
+                else if (logratio < 0.0) logratio = 0.0;
+                if (BETA_FACTOR > 1.0) hue = PARTICLE_HUE_MAX - (PARTICLE_HUE_MAX - PARTICLE_HUE_MIN)*logratio;
+                else hue = PARTICLE_HUE_MIN - (PARTICLE_HUE_MIN - PARTICLE_HUE_MAX)*logratio;
+            }
+            else hue = 0.25*PARTICLE_HUE_MIN + 0.75*PARTICLE_HUE_MAX;
+            erase_area_hsl_turbo(0.0, YMIN, 2.0*PARTIAL_THERMO_WIDTH, PARTIAL_THERMO_HEIGHT*(YMAX - YMIN),  hue, 0.9, 0.15);
+        }
+        else if (PARTIAL_THERMO_REGION == TH_LAYER) 
+        {
+            hue = 0.75*PARTICLE_HUE_MIN + 0.25*PARTICLE_HUE_MAX;
+            erase_area_hsl_turbo(0.0, YMIN, XMAX, PARTIAL_THERMO_HEIGHT - YMIN,  hue, 0.9, 0.15);
         }
-        else hue = 0.25*PARTICLE_HUE_MIN + 0.75*PARTICLE_HUE_MAX;
-        erase_area_hsl_turbo(0.0, YMIN, 2.0*PARTIAL_THERMO_WIDTH, PARTIAL_THERMO_HEIGHT*(YMAX - YMIN),  hue, 0.9, 0.15);
-        
     }
     
     /* draw "altitude lines" */
@@ -5113,19 +5508,30 @@ void draw_container(double xmin, double xmax, t_obstacle obstacle[NMAXOBSTACLES]
 /* draw the container, for certain boundary conditions */
 {
     int i, j;
-    double rgb[3], x, phi, angle, dx, dy, ybin, x1, x2, h;
+    double rgb[3], x, y, r, phi, angle, dx, dy, ybin, x1, x2, h;
     char message[100];
     
     /* draw fixed obstacles */
     if (ADD_FIXED_OBSTACLES)
     {
         glLineWidth(CONTAINER_WIDTH);
-        hsl_to_rgb(300.0, 0.1, 0.5, rgb);
+        if (CHARGE_OBSTACLES) hsl_to_rgb(30.0, 0.1, 0.5, rgb);
+        else hsl_to_rgb(300.0, 0.1, 0.5, rgb);
         for (i = 0; i < nobstacles; i++)
             draw_colored_circle_precomp(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius, rgb);
         glColor3f(1.0, 1.0, 1.0);
         for (i = 0; i < nobstacles; i++)
-            draw_circle_precomp(obstacle[i].xc - xtrack, obstacle[i].yc - ytrack, obstacle[i].radius);
+        {
+            x = obstacle[i].xc - xtrack;
+            y = obstacle[i].yc - ytrack;
+            r = obstacle[i].radius;
+            draw_circle_precomp(x, y, r);
+            if (CHARGE_OBSTACLES) 
+            {
+                draw_line(x - r, y, x + r, y);
+                draw_line(x, y - r, x, y + r);
+            }
+        }
     }
     if (ADD_FIXED_SEGMENTS)
     {
@@ -5433,7 +5839,8 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
             xxtext = XMIN + 0.08*scale;
         }
         xmid = 0.5*(XMIN + XMAX) - 0.1*scale;
-        xmidtext = xmid - 0.3*scale;
+//         xmidtext = xmid - 0.3*scale;
+        xmidtext = xmid - 0.2*scale;
         for (i=0; i<N_P_AVERAGE; i++) pressures[i] = 0.0;
         if (RECORD_PRESSURES) for (j=0; j<N_PRESSURES; j++) 
         {
@@ -5513,12 +5920,12 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
         erase_area_hsl(xbox, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0);
         glColor3f(1.0, 1.0, 1.0);
         sprintf(message, "Density %.3f", density);
-        write_text(xtext, y, message);
+        write_text(xtext + 0.1, y, message);
         
         erase_area_hsl(xmid, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0);
         glColor3f(1.0, 1.0, 1.0);
         sprintf(message, "Temperature %.2f", params.mean_energy);
-        write_text(xmidtext, y, message);
+        write_text(xmidtext - 0.1, y, message);
 
         erase_area_hsl(xxbox, y + 0.025*scale, 0.37*scale, 0.05*scale, 0.0, 0.9, 0.0);
         glColor3f(1.0, 1.0, 1.0);
@@ -5540,6 +5947,28 @@ void print_parameters(t_lj_parameters params, short int left, double pressure[N_
         sprintf(message, "E field = %.2f", 25.0*NVID*DT_PARTICLE*params.efield);
         write_text(xtext + 0.08, y, message);
     }   
+    else if (INCREASE_B)  /* print magnetic field */
+    {
+        erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "B field = %.2f", 25.0*NVID*DT_PARTICLE*params.bfield);
+        write_text(xtext + 0.08, y, message);
+    }   
+    
+    if (PRINT_NPARTICLES)   /* print number of particles */
+    {
+        erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "%i electrons", params.nactive);
+        write_text(xtext + 0.08, y, message);
+    } 
+    if (PRINT_TYPE_PROP)   /* print proportion of types */
+    {
+        erase_area_hsl(xbox, y + 0.025*scale, 0.27*scale, 0.05*scale, 0.0, 0.9, 0.0);
+        glColor3f(1.0, 1.0, 1.0);
+        sprintf(message, "%.1f %% anions", 100.0*params.prop);
+        write_text(xtext + 0.08, y, message);
+    } 
     
     if (RECORD_PRESSURES) 
     {
@@ -6015,10 +6444,10 @@ void print_particles_speeds(t_particle particle[NMAXCIRCLES])
     write_text(xrighttext + 0.1, y, message);
 }
 
-double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall)
+double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacle obstacle[NMAXOBSTACLES], t_segment segment[NMAXSEGMENTS], double xleft, double xright, double *pleft, double *pright, double pressure[N_PRESSURES], int wall, double krepel)
 {
     int i, k;
-    double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin;
+    double xmin, xmax, ymin, ymax, padding, r, rp, r2, cphi, sphi, f, fperp = 0.0, x, y, xtube, distance, dx, dy, width, ybin, angle, x1, x2, h, ytop, norm, dleft, dplus, dminus, tmp_pleft = 0.0, tmp_pright = 0.0, proj, pscal, pvect, pvmin, charge, d2;
     
     /* compute force from fixed circular obstacles */
     if (ADD_FIXED_OBSTACLES) for (i=0; i<nobstacles; i++)
@@ -6037,6 +6466,15 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
             particle[j].fx += f*cphi;
             particle[j].fy += f*sphi;
         }
+        
+        if (CHARGE_OBSTACLES)
+        {
+            charge = obstacle[i].charge*particle[k].charge;
+            d2 = distance*distance;
+            f = KCOULOMB_OBSTACLE*charge/(1.0e-12 + d2);
+            particle[j].fx += f*cphi;
+            particle[j].fy += f*sphi;
+        }
     }
     /* compute force from fixed linear obstacles */
     particle[j].close_to_boundary = 0;
@@ -6550,29 +6988,88 @@ double compute_boundary_force(int j, t_particle particle[NMAXCIRCLES], t_obstacl
 }
 
 
+void compute_partner_force(int j, int n, double eq_distance, double f[2], t_particle particle[NMAXCIRCLES])
+/* compute force of partner particle n on particle j */
+{
+    double dx, dy, r, ca, sa, force;
+    
+    dx = particle[n].xc - particle[j].xc;
+    dy = particle[n].yc - particle[j].yc;
+    
+    if (dx > 0.5*(BCXMAX - BCXMIN)) dx -= (BCXMAX-BCXMIN);
+    else if (dx < -0.5*(BCXMAX - BCXMIN)) dx += (BCXMAX-BCXMIN);
+    if (dy > 0.5*(BCYMAX - BCYMIN)) dy -= (BCYMAX-BCYMIN);
+    else if (dy < -0.5*(BCYMAX - BCYMIN)) dy += (BCYMAX-BCYMIN);
+    
+    r = module2(dx, dy);
+    if (r < 1.0e-10) r = 1.0e-10;
+    if (r > 0.5) r = 0.5;
+    ca = dx/r;
+    sa = dy/r;
+    
+    force = KSPRING_PAIRS*(r - eq_distance);
+    f[0] = force*ca;
+    f[1] = force*sa;
+}
+
+
 void compute_particle_force(int j, double krepel, t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY])
 /* compute force from other particles on particle j */
 {
-    int i0, j0, m0, k, m, q, close;
+    int i0, j0, m0, k, l, m, q, close, n, test;
     double fx = 0.0, fy = 0.0, force[2], torque = 0.0, torque_ij, x, y;
     
     particle[j].neighb = 0;
     if (REACTION_DIFFUSION) particle[j].diff_neighb = 0;
 
-    /* NEW */
     for (k=0; k<particle[j].hash_nneighb; k++) 
-        if (particle[particle[j].hashneighbour[k]].active)
     {
-        close = compute_repelling_force(j, k, force, &torque_ij, particle, krepel);
+        n = particle[j].hashneighbour[k];
+        if (particle[n].active)
+        {
+            if (PAIR_PARTICLES)
+            {
+                /* test whether j is not a partner particle */
+                test = 1;
+                for (l=0; l<particle[n].npartners; l++)
+                    if (particle[n].partner[l] == j) test = 0;
+                
+                if (test)
+                {
+                    close = compute_repelling_force(j, k, force, &torque_ij, particle, krepel);
+                    fx += force[0];
+                    fy += force[1];
+                    torque += torque_ij;
+                    if (close) 
+                    {
+                        particle[j].neighb++;
+                        if (REACTION_DIFFUSION&&(particle[j].type != particle[particle[j].hashneighbour[k]].type)) 
+                            particle[j].diff_neighb++;
+                    }
+                }
+            }
+            else
+            {
+                close = compute_repelling_force(j, k, force, &torque_ij, particle, krepel);
+                fx += force[0];
+                fy += force[1];
+                torque += torque_ij;
+                if (close) 
+                {
+                    particle[j].neighb++;
+                    if (REACTION_DIFFUSION&&(particle[j].type != particle[particle[j].hashneighbour[k]].type)) 
+                        particle[j].diff_neighb++;
+                }
+            }
+        }
+    }
+    
+    if (PAIR_PARTICLES) for (l=0; l<particle[j].npartners; l++)
+    {
+        n = particle[j].partner[l];
+        compute_partner_force(j, n, particle[j].partner_eqd[l], force, particle);
         fx += force[0];
         fy += force[1];
-        torque += torque_ij;
-        if (close) 
-        {
-            particle[j].neighb++;
-            if (REACTION_DIFFUSION&&(particle[j].type != particle[particle[j].hashneighbour[k]].type)) 
-                particle[j].diff_neighb++;
-        }
     }
                 
     particle[j].fx += fx;
@@ -6654,19 +7151,108 @@ int twotype_config(int i, t_particle particle[NMAXCIRCLES])
 }
 
 
+void init_particle_pairs(t_particle particle[NMAXCIRCLES])
+/* initialize data structure for paired particles */
+{
+    int i, k, l, n, counter;
+    double angle, dist;
+    
+    if (ncircles*NPARTNERS > NMAXCIRCLES)
+    {
+        printf("Error: NMAXCIRCLES is too small\n");
+        exit(1); 
+    }
+    
+    for (i=0; i<ncircles; i++) if (particle[i].active)
+    {
+        for (k=0; k<NPARTNERS; k++)
+        {
+            n = ncircles*(k+1) + i;
+            angle = DPI*(double)rand()/(double)RAND_MAX + DPI*(double)k/(double)NPARTNERS;
+            particle[n].type = 2;
+            particle[n].xc = particle[i].xc + (MU + MU_B)*cos(angle);
+            particle[n].yc = particle[i].yc + (MU + MU_B)*sin(angle);
+                
+            particle[n].radius = MU_B;
+            particle[n].mass_inv = 1.0/PARTICLE_MASS_B;
+            particle[n].thermostat = 1;
+            particle[n].charge = CHARGE_B;
+            particle[n].active = 1;
+        
+            particle[n].partner[0] = i;
+            particle[i].partner[k] = n;
+            
+            particle[n].partner_eqd[0] = MU + MU_B;
+            particle[i].partner_eqd[k] = MU + MU_B;
+        }
+        particle[i].npartners = NPARTNERS;
+        
+        switch (PAIRING_TYPE) {
+            case (POLY_STAR):
+            {
+                for (k=0; k<NPARTNERS; k++)
+                {
+                    n = ncircles*(k+1) + i;
+                    particle[n].npartners = 1;
+                }
+                break;
+            }
+            case (POLY_ALL):
+            {
+                for (k=0; k<NPARTNERS; k++)
+                {
+                    n = ncircles*(k+1) + i;
+                    particle[n].npartners = NPARTNERS;
+                    counter = 1;
+                    for (l=0; l<NPARTNERS; l++) if (l != k)
+                    {
+                        particle[n].partner[counter] = ncircles*(l+1) + i;
+                        particle[n].partner_eqd[counter] = MU + MU_B;
+                        counter++;
+                    }
+                }
+                break;
+            }
+            case (POLY_WATER):
+            {
+//                 dist = 2.0*sin(PARTNER_ANGLE*PI/360.0)*(MU + MU_B);
+                for (k=0; k<NPARTNERS; k++)
+                {
+                    n = ncircles*(k+1) + i;
+                    particle[n].npartners = NPARTNERS;
+                    counter = 1;
+                    for (l=0; l<NPARTNERS; l++) if (l != k)
+                    {
+                        particle[n].partner[counter] = ncircles*(l+1) + i;
+                        dist = 2.0*sin(vabs((double)(l-k))*PARTNER_ANGLE*PI/360.0)*(MU + MU_B);
+                        particle[n].partner_eqd[counter] = dist;
+                        counter++;
+                    }
+                }
+                break;
+            }
+            
+        }
+    }
+    ncircles += ncircles*NPARTNERS;
+    printf("ncircles = %i\n", ncircles);
+}
+
+
 int initialize_configuration(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY], t_obstacle obstacle[NMAXOBSTACLES], double px[NMAXCIRCLES], double py[NMAXCIRCLES], double pangle[NMAXCIRCLES], int tracer_n[N_TRACER_PARTICLES],
 t_segment segment[NMAXSEGMENTS])
 /* initialize all particles, obstacles, and the hashgrid */
 {
-    int i, j, k, n, type, nactive = 0;
-    double x, y, h, xx, yy, rnd;
+    int i, j, k, n, type, nactive = 0, hashcell;
+    double x, y, h, xx, yy, rnd, angle;
     
     for (i=0; i < ncircles; i++) 
     {
         /* set particle type */
         particle[i].type = 0;
 //         if ((TWO_TYPES)&&((double)rand()/RAND_MAX > TYPE_PROPORTION)) 
-        if ((TWO_TYPES)&&(twotype_config(i, particle))) 
+//         if ((TWO_TYPES)&&(twotype_config(i, particle))) 
+        if (TWO_TYPES)
         {
             type = twotype_config(i, particle);
             if (type == 1)
@@ -6780,6 +7366,7 @@ t_segment segment[NMAXSEGMENTS])
         else if ((PLOT == P_NUMBER)||(PLOT_B == P_NUMBER))
             particle[i].color_hue = 360.0*(double)(i%N_PARTICLE_COLORS)/(double)N_PARTICLE_COLORS;
     }
+    
     /* initialize dummy values in case particles are added */
     for (i=ncircles; i < NMAXCIRCLES; i++) 
     {
@@ -6878,6 +7465,9 @@ t_segment segment[NMAXSEGMENTS])
             particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
             particle[i].radius = MU_B;
         }
+        
+    /* add copies in case of particle pairing */
+    if (PAIR_PARTICLES) init_particle_pairs(particle);
     
     /* inactivate particles in obstacle */
     printf("Inactivating particles inside obstacles\n");
@@ -6931,8 +7521,13 @@ t_segment segment[NMAXSEGMENTS])
     if (ADD_FIXED_OBSTACLES)
     {
         for (i=0; i< ncircles; i++) for (j=0; j < nobstacles; j++)
+        {
             if (module2(particle[i].xc - obstacle[j].xc, particle[i].yc - obstacle[j].yc) < OBSTACLE_RADIUS + particle[i].radius)
                 particle[i].active = 0;
+        
+            hashcell = hash_cell(obstacle[j].xc, obstacle[j].yc);
+            hashgrid[hashcell].charge = OBSTACLE_CHARGE;
+        }
     }
     
     /* case of segment obstacles */
@@ -7136,16 +7731,17 @@ int add_particles(t_particle particle[NMAXCIRCLES], double px[NMAXCIRCLES], doub
     
     x = ADDXMIN + (ADDXMAX - ADDXMIN)*(double)rand()/(double)RAND_MAX;
     y = ADDYMIN + 0.5*(ADDYMAX - ADDYMIN)*(double)rand()/(double)RAND_MAX;
-    add_particle(x, y, 0.0, 0.0, PARTICLE_MASS, 0, particle);
+    add_particle(x, y, V_INITIAL, 0.0, PARTICLE_MASS, 0, particle);
 //     add_particle(BCXMIN + 0.1, 0.5*(BCYMIN + BCYMAX), 200.0, 0.0, PARTICLE_MASS, 0, particle);
 //     i++;
     
     particle[ncircles - 1].radius = MU;
     particle[ncircles - 1].eq_dist = EQUILIBRIUM_DIST;
-    particle[ncircles - 1].thermostat = 0;
+    particle[ncircles - 1].thermostat = 1;
     px[ncircles - 1] = particle[ncircles - 1].vx;
     py[ncircles - 1] = particle[ncircles - 1].vy;
-    return (nadd_particle + 1);
+//     return (nadd_particle + 1);
+    return (ncircles);
 }
 
 
@@ -7963,6 +8559,67 @@ int chem_catalytic_convert(int i, int type2, int newtype, t_particle particle[NM
     return(ncollisions); 
 }
 
+void change_type_proportion(t_particle particle[NMAXCIRCLES], double prop)
+/* change proportion of particles of types 1 or 2 */
+{
+    int i, n0=0, n1=0, n2=0, ntot, nc=0, nmod, counter=0, cmax=1000;
+    
+    for (i=0; i<ncircles; i++) if (particle[i].active)
+    {
+        if (particle[i].type == 0) n0++;
+        if (particle[i].type == 1) n1++;
+        if (particle[i].type == 2) n2++;
+    }
+    
+    ntot = n0 + n1 + n2;
+    n0 += n1;
+    
+    if ((double)(n0+1) < prop*(double)ntot)
+    {
+        nmod = (int)((double)ntot*prop - (double)n0);
+        if (nmod > 0) while ((nc<nmod)&&(counter<cmax))
+        {   
+            i = rand()%ncircles;
+            if ((particle[i].type == 2)&&(particle[i].active))
+            {
+                particle[i].type = 1;
+                particle[i].radius = MU;
+                particle[i].interaction = INTERACTION;
+                particle[i].eq_dist = EQUILIBRIUM_DIST;
+                particle[i].spin_range = SPIN_RANGE;
+                particle[i].spin_freq = SPIN_INTER_FREQUENCY;
+                particle[i].mass_inv = 1.0/PARTICLE_MASS;
+                particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT;
+                particle[i].charge = CHARGE;
+                nc++;
+            }
+            counter++;
+        }
+    }
+    else if ((double)(n0-1) > prop*(double)ntot)
+    {
+        nmod = (int)((double)n0 - (double)ntot*prop);
+        if (nmod > 0) while ((nc<nmod)&&(counter<cmax))
+        {   
+            i = rand()%ncircles;
+            if ((particle[i].type <= 1)&&(particle[i].active)) 
+            {
+                particle[i].type = 2;
+                particle[i].radius = MU_B;
+                particle[i].interaction = INTERACTION_B;
+                particle[i].eq_dist = EQUILIBRIUM_DIST_B;
+                particle[i].spin_range = SPIN_RANGE_B;
+                particle[i].spin_freq = SPIN_INTER_FREQUENCY_B;
+                particle[i].mass_inv = 1.0/PARTICLE_MASS_B;
+                particle[i].inertia_moment_inv = 1.0/PARTICLE_INERTIA_MOMENT_B;
+                particle[i].charge = CHARGE_B;
+                nc++;
+            }
+            counter++;
+        }
+    }
+}
+
 int update_types(t_particle particle[NMAXCIRCLES], t_collision *collisions, int ncollisions, int *particle_numbers, int time, double *delta_e)
 /* update the types in case of reaction-diffusion equation */
 {

sub_sphere.c

@@ -71,9 +71,15 @@ void init_wave_sphere(t_wave_sphere wsphere[NX*NY])
         }
         
         /* cotangent, taking care of not dividing by zero */
-        for (j=1; j<NY-1; j++) wsphere[i*NY+j].cottheta = wsphere[i*NY+j].ctheta/wsphere[i*NY+j].stheta;
-        wsphere[i*NY].cottheta = wsphere[i*NY+1].cottheta;
-        wsphere[i*NY + NY-1].cottheta = wsphere[i*NY+1].cottheta;
+        /* TODO clean up cottheta range ? */
+        for (j=DPOLE; j<NY-DPOLE; j++) wsphere[i*NY+j].cottheta = wsphere[i*NY+j].ctheta/wsphere[i*NY+j].stheta;
+        for (j=0; j<DPOLE; j++) wsphere[i*NY+j].cottheta = wsphere[i*NY+DPOLE].cottheta;
+        for (j=NY-DPOLE; j<NY; j++) wsphere[i*NY+j].cottheta = wsphere[i*NY+DPOLE-1].cottheta;        
+        
+        /* old version */
+//         for (j=1; j<NY-1; j++) wsphere[i*NY+j].cottheta = wsphere[i*NY+j].ctheta/wsphere[i*NY+j].stheta;
+//         wsphere[i*NY].cottheta = wsphere[i*NY+1].cottheta;
+//         wsphere[i*NY + NY-1].cottheta = wsphere[i*NY+1].cottheta;
     }
 }
 
@@ -229,9 +235,9 @@ void read_negative_dem_values(double *height_values, t_wave_sphere wsphere[NX*NY
 void init_dem(t_wave_sphere wsphere[NX*NY], int dem_number)
 /* init heights from digital elevation map */
 {
-    int i, j, ii, jj, k, nx, ny, maxrgb, nmaxpixels = 4915200, scan, rgbval, diff, sshift, nshift, hmin, hmax, ishift;
+    int i, j, ii, jj, k, nx, ny, maxrgb, nmaxpixels = 4915200, scan, rgbval, diff, sshift, nshift, hmin, hmax, ishift, hsum;
     int *rgb_values;
-    double cratio, rx, ry, cy, dx, dy, pscal, norm, vscale1, vscale2, gradx, grady, deltar, deltai[3], deltaj[3], dphi, dtheta, n[3], hsea;
+    double cratio, rx, ry, cy, dx, dy, pscal, norm, vscale1, vscale2, gradx, grady, deltar, deltai[3], deltaj[3], dphi, dtheta, n[3], hsea, hmean;
     double *height_values, *height_values_tmp;
     FILE *image_file;
     
@@ -259,6 +265,20 @@ void init_dem(t_wave_sphere wsphere[NX*NY], int dem_number)
             hsea = 255.0*PLANET_SEALEVEL/DEM_MAXHEIGHT;
             break;
         }
+        case (DEM_VENUS):
+        {
+            nmaxpixels = 4096*2048;
+            image_file = fopen("Venus_Magellan_Topography_Global_4641m_v02_scaled2.ppm", "r");
+            hsea = 255.0*PLANET_SEALEVEL/DEM_MAXHEIGHT;
+            break;
+        }
+        case (DEM_MERCURY):
+        {
+            nmaxpixels = 1280*380;
+            image_file = fopen("Mercury_Messenger_DEM_Global_665m_1024_1_cropped.ppm", "r");
+            hsea = 255.0*PLANET_SEALEVEL/DEM_MAXHEIGHT;
+            break;
+        }
     }
     
     rgb_values = (int *)malloc(3*nmaxpixels*sizeof(int));
@@ -269,6 +289,9 @@ void init_dem(t_wave_sphere wsphere[NX*NY], int dem_number)
     scan = fscanf(image_file,"%i %i\n", &nx, &ny);
     scan = fscanf(image_file,"%i\n", &maxrgb);
     
+    printf("nx = %i, ny = %i, maxrgb = %i\n", nx, ny, maxrgb);
+    sleep(1);
+    
     hmin = maxrgb;
     hmax = 0;
     
@@ -294,12 +317,29 @@ void init_dem(t_wave_sphere wsphere[NX*NY], int dem_number)
                 {
                     scan = fscanf(image_file,"%i\n", &rgbval);
                     rgb_values[3*(j*nx+i)+k] = rgbval;
-                    if (rgbval < hmin) hmin = rgbval;
+                    if (rgbval < hmin) 
+                    {
+                        if (B_DOMAIN == D_SPHERE_VENUS)
+                        {
+                            if (rgbval > 0) hmin = rgbval;
+                        }
+                        else hmin = rgbval;
+                    }
                     if (rgbval > hmax) hmax = rgbval;
                 }
     
     printf("hmin = %i, hmax = %i, hsea = %i\n", hmin, hmax, (int)hsea);
-                        
+    
+    if (B_DOMAIN == D_SPHERE_VENUS)
+    {
+        hsum = 0;
+        for (j=0; j<ny; j++)
+            for (i=0; i<nx; i++)
+                hsum += rgb_values[3*(j*nx+i)];
+        hmean = (double)hsum/(double)(nx*ny);
+        printf("hmean = %.2f\n", hmean);
+    }
+    
     cratio = 1.0/(double)(hmax-hmin);
     rx = (double)nx/(double)NX;
     ry = (double)ny/(double)(NY - sshift - nshift);
@@ -317,6 +357,13 @@ void init_dem(t_wave_sphere wsphere[NX*NY], int dem_number)
             height_values[i*NY+j] = ((double)rgb_values[3*(jj*nx+ii)]-hsea)*cratio;
             wsphere[i*NY+j].altitude = ((double)rgb_values[3*(jj*nx+ii)]-hsea)*cratio;
             
+            /* take care of black areas (missing data) on venus */
+            if ((B_DOMAIN == D_SPHERE_VENUS)&&(rgb_values[3*(jj*nx+ii)] == 0))
+            {
+                height_values[i*NY+j] = VENUS_NODATA_FACTOR*hmean*cratio;
+                wsphere[i*NY+j].altitude = VENUS_NODATA_FACTOR*hmean*cratio;
+            }
+            
             if (OTHER_PLANET)
                 wsphere[i*NY+j].indomain = (wsphere[i*NY+j].altitude < 0.0);
             
@@ -712,6 +759,22 @@ void init_planet_map(t_wave_sphere wsphere[NX*NY], int planet)
             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);
@@ -833,6 +896,34 @@ int ij_to_sphere(int i, int j, double r, t_wave_sphere wsphere[NX*NY], double xy
     return(pscal/norm_observer > COS_VISIBLE);
 }
 
+
+int ij_to_sphere_r(int i, int j, double r, t_wave_sphere wsphere[NX*NY], double xyz[3])
+/* convert spherical to rectangular coordinates */
+{
+    double pscal, newr;
+    static double norm_observer;
+    static int first = 1;
+    
+    if (first)
+    {
+        norm_observer = sqrt(observer[0]*observer[0] + observer[1]*observer[1] + observer[2]*observer[2]);
+        first = 0;
+    }
+    
+    xyz[0] = wsphere[i*NY+j].x;
+    xyz[1] = wsphere[i*NY+j].y;
+    xyz[2] = wsphere[i*NY+j].z;
+    
+    pscal = xyz[0]*observer[0] + xyz[1]*observer[1] + xyz[2]*observer[2];
+    
+    xyz[0] *= r;
+    xyz[1] *= r;
+    xyz[2] *= r;
+    
+    return(pscal/norm_observer > COS_VISIBLE);
+}
+
+
 int init_circle_sphere(t_circles_sphere *circles, int circle_pattern)
 /* initialise circles on sphere */
 /* for billiard shape D_SPHERE_CIRCLES */
@@ -987,6 +1078,224 @@ int init_circle_sphere(t_circles_sphere *circles, int circle_pattern)
     return(ncircles);
 }
 
+void init_sphere_maze(t_wave_sphere wsphere[NX*NY], int spiral, int wave, int npole)
+/* initialize maze on sphere */
+{
+    int nblocks, block, i, j, k, n, p, q, np, na, inrect;
+    double rmin, rmax, angle, r, dr, phi, dphi, phi1, ww, width = 0.02, theta, theta1, phaseshift; 
+    t_maze* maze;
+    t_rectangle* polyrect;
+    
+    maze = (t_maze *)malloc(NXMAZE*NYMAZE*sizeof(t_maze));
+    polyrect = (t_rectangle *)malloc(NMAXPOLY*sizeof(t_rectangle));
+    
+    printf("Initializing maze\n");
+    init_circular_maze(maze);
+    printf("Maze initialized\n");
+    
+    printf("Building maze walls\n");
+    /* build walls of maze */
+    
+    np = 0;
+    na = 0;
+    
+    /* build walls of maze */
+    nblocks = NYMAZE/NXMAZE;
+    rmin = 0.3;
+    rmax = PID + 0.2;
+    angle = DPI/(double)nblocks;
+        
+    dr = (rmax - rmin)/(double)(NXMAZE);
+    
+    /* add straight walls */
+    for (block = 0; block < nblocks; block++)
+    {
+        printf("Block %i\n", block);
+        dphi = angle;
+        
+        /* first circle */
+        n = nmaze(0, block*NXMAZE);
+        r = rmin - 0.5*width;
+        phi = (double)block*angle;
+            
+        if (maze[n].south)
+        {
+            polyrect[np].x1 = phi - 0.5*width;
+            polyrect[np].y1 = r;
+            polyrect[np].x2 = phi + 0.5*width;
+            polyrect[np].y2 = r+dr+width;
+//             printf("Adding rectangle at (%.3f, %.3f) - (%.3f, %.3f)\n", polyrect[np].x1, polyrect[np].y1, polyrect[np].x2, polyrect[np].y2);
+            np++;
+        }
+                
+        /* second circle */
+        r = rmin + dr - 0.5*width;
+        dphi *= 0.5;
+        for (q=0; q<2; q++)
+        {
+            n = nmaze(1, block*NXMAZE + q);
+            phi = (double)(block)*angle + (double)q*dphi;
+            
+            if (maze[n].south)
+            {
+                polyrect[np].x1 = phi - 0.5*width;
+                polyrect[np].y1 = r;
+                polyrect[np].x2 = phi + 0.5*width;
+                polyrect[np].y2 = r+dr+width;
+//                 printf("Adding rectangle at (%.3f, %.3f) - (%.3f, %.3f)\n", polyrect[np].x1, polyrect[np].y1, polyrect[np].x2, polyrect[np].y2);
+                np++;
+            }
+        }
+                
+        /* other circles */
+        ww = 2;
+        i = 2;
+        while (ww < NXMAZE)
+        {
+            dphi *= 0.5;
+            for (p = 0; p < ww; p++)
+            {
+                r = rmin + (double)i*dr - 0.5*width;
+//                 printf("Segment, i = %i, dphi = %.2lg, r = %.2lg\n", i, dphi, r);
+                for (q = 0; q < 2*ww; q++)
+                {
+                    j = block*NXMAZE + q;
+                    n = nmaze(i,j);
+                    phi = (double)(block)*angle + (double)q*dphi;
+                    
+                    if (maze[n].south)
+                    {
+                        polyrect[np].x1 = phi - 0.5*width;
+                        polyrect[np].y1 = r;
+                        polyrect[np].x2 = phi + 0.5*width;
+                        polyrect[np].y2 = r+dr+width;
+//                         printf("Adding rectangle at (%.3f, %.3f) - (%.3f, %.3f)\n", polyrect[np].x1, polyrect[np].y1, polyrect[np].x2, polyrect[np].y2);
+                        np++;
+                    }
+                }
+                i++;
+            }
+            ww *= 2;
+        }
+                
+    }
+    
+    /* add circular arcs */
+    for (block = 0; block < nblocks; block++)
+    {
+        dphi = angle;
+        
+        /* first circle */
+        n = nmaze(0, block*NXMAZE);
+        r = rmin;
+        phi = (double)block*angle;
+        
+        if ((block > 0)&&(maze[n].west))
+        {
+            polyrect[np].x1 = phi;
+            polyrect[np].y1 = r - 0.5*width;
+            polyrect[np].x2 = phi + dphi + 0.05*width;
+            polyrect[np].y2 = r + 0.5*width;
+            np++;
+        }
+                
+        /* second circle */
+        r = rmin + dr;
+        dphi *= 0.5;
+        for (q=0; q<2; q++)
+        {
+            n = nmaze(1, block*NXMAZE + q);
+            phi = (double)(block)*angle + (double)q*dphi;
+            
+            if (maze[n].west)
+            {
+                polyrect[np].x1 = phi;
+                polyrect[np].y1 = r - 0.5*width;
+                polyrect[np].x2 = phi + dphi + 0.05*width;
+                polyrect[np].y2 = r + 0.5*width;
+                np++;
+            }
+        }
+                
+        /* other circles */
+        ww = 2;
+        i = 2;
+        while (ww < NXMAZE)
+        {
+            dphi *= 0.5;
+            for (p = 0; p < ww; p++)
+            {
+                r = rmin + (double)i*dr;
+                printf("Circle, i = %i, dphi = %.2lg, r = %.2lg\n", i, dphi, r);
+                for (q = 0; q < 2*ww; q++)
+                {
+                    j = block*NXMAZE + q;
+                    n = nmaze(i,j);
+                    phi = (double)(block)*angle + (double)q*dphi;
+                    
+                    if (maze[n].west)
+                    {
+                        polyrect[np].x1 = phi;
+                        polyrect[np].y1 = r - 0.5*width;
+                        polyrect[np].x2 = phi + dphi + 0.05*width;
+                        polyrect[np].y2 = r + 0.5*width;
+                        np++;
+                    }
+                }
+                i++;
+            }
+            ww *= 2;
+        }
+    }
+    
+    /* outer boundary of maze */
+    polyrect[np].x1 = dphi;
+    polyrect[np].y1 = rmax - 0.5*width;
+    polyrect[np].x2 = DPI;
+    polyrect[np].y2 = rmax + 0.5*width;
+    np++;
+    
+
+    printf("Maze walls built\n");
+    
+    phaseshift = ZERO_MERIDIAN*PI/180.0;
+    for (i=0; i<NX; i++)
+    {
+        phi = DPI*(double)i/(double)NX + phaseshift;
+        if (phi > DPI) phi -= DPI;
+        for (j=0; j<NY; j++)
+        {
+            theta = PI*(1.0 - (double)j/(double)NY);
+            if (spiral) 
+            {
+                phi -= PID/(double)NY;
+                if (phi < 0.0) phi += DPI;
+            }
+            if (wave)
+            {
+                theta1 = theta; 
+                theta += 0.03*cos(6.0*phi)*theta;
+                if (theta > PI) theta = PI; 
+                phi += 0.00075*cos(20.0*theta1)*theta1;
+//                 if (phi > DPI) phi -= DPI;
+                if (phi < 0.0) phi += DPI;
+            }
+            inrect = 0;
+            for (n=0; n<np; n++) for (k=-1; k<2; k++)
+            {
+                phi1 = phi + k*DPI;
+                if (((phi1 > polyrect[n].x1)&&(phi1 < polyrect[n].x2)&&(theta > polyrect[n].y1)&&(theta < polyrect[n].y2))) inrect = 1;
+            }
+            wsphere[i*NY+j].indomain = 1-inrect;
+        }
+        if (npole) for (j=NY-DPOLE*3/2; j<NY; j++) wsphere[i*NY+j].indomain = 0;
+    }
+    
+    free(maze);
+    free(polyrect);
+    
+}
+
 int xy_in_billiard_sphere(int i, int j, t_wave_sphere wsphere[NX*NY])
 /* returns 1 if (x,y) represents a point in the billiard */
 {
@@ -1004,6 +1313,9 @@ int xy_in_billiard_sphere(int i, int j, t_wave_sphere wsphere[NX*NY])
             cos_rot = cos(JULIA_ROT*DPI/360.0);
             sin_rot = sin(JULIA_ROT*DPI/360.0);
         }
+        else if (B_DOMAIN == D_SPHERE_MAZE) init_sphere_maze(wsphere, 0, 0, 0);
+        else if (B_DOMAIN == D_SPHERE_MAZE_SPIRAL) init_sphere_maze(wsphere, 1, 0, 1);
+        else if (B_DOMAIN == D_SPHERE_MAZE_WAVE) init_sphere_maze(wsphere, 0, 1, 1);
         first = 0;
     }
     
@@ -1088,22 +1400,9 @@ int xy_in_billiard_sphere(int i, int j, t_wave_sphere wsphere[NX*NY])
             if (u*u + v*v < MANDELLIMIT) return(1);
             return(0);
         }
-        case (D_SPHERE_EARTH):
-        {
-            return(wsphere[i*NY+j].indomain);
-        }
-        case (D_SPHERE_MARS):
-        {
-            return(wsphere[i*NY+j].indomain);
-        }
-        case (D_SPHERE_MOON):
-        {
-            return(wsphere[i*NY+j].indomain);
-        }
         default:
         {
-            printf("Function xy_in_billiard_sphere not defined for this billiard \n");
-            return(0);
+            return(wsphere[i*NY+j].indomain);
         }
     }
 }
@@ -1204,19 +1503,23 @@ void init_circular_wave_sphere(double phi0, double theta0, double phi[NX*NY], do
 /* initialise field with drop - phi is wave height, psi is phi at time t-1 */
 /* phi0 is longidude, theta0 is latitude */
 {
-    int i, j;
+    int i, j, jmin, jmax;
     double xy[2], pscal, dist, dist2, stheta, ctheta, phishift;
 
     ctheta = cos(theta0 + PID);
     stheta = sin(theta0 + PID);
     phishift = ZERO_MERIDIAN*PI/180.0;
+    
+    /* safety distance to avoid singularity at poles */
+    jmin = DPOLE;
+    jmax = NY-DPOLE;
         
     printf("Initializing wave\n"); 
 //     #pragma omp parallel for private(i,j,xy,dist2)
     for (i=0; i<NX; i++)
     {
         if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
-        for (j=0; j<NY; j++)
+        for (j=jmin; j<jmax; j++)
         {
             pscal = cos(wsphere[i*NY+j].phi - phi0 - phishift)*wsphere[i*NY+j].stheta*stheta;
             pscal += wsphere[i*NY+j].ctheta*ctheta;
@@ -1231,6 +1534,16 @@ void init_circular_wave_sphere(double phi0, double theta0, double phi[NX*NY], do
             else phi[i*NY+j] = 0.0;
             psi[i*NY+j] = 0.0;
         }
+        for (j=0; j<jmin; j++)
+        {
+            phi[i*NY+j] = 0.0;
+            psi[i*NY+j] = 0.0;
+        }
+        for (j=jmax; j<NY; j++)
+        {
+            phi[i*NY+j] = 0.0;
+            psi[i*NY+j] = 0.0;
+        }
     }
 }
 
@@ -1238,19 +1551,21 @@ void add_circular_wave_sphere(double amp, double phi0, double theta0, double phi
 /* initialise field with drop - phi is wave height, psi is phi at time t-1 */
 /* phi0 is longidude, theta0 is latitude */
 {
-    int i, j;
+    int i, j, jmin, jmax;
     double xy[2], pscal, dist, dist2, stheta, ctheta, phishift;
 
     ctheta = cos(theta0 + PID);
     stheta = sin(theta0 + PID);
     phishift = ZERO_MERIDIAN*PI/180.0;
     
-    printf("Initializing wave\n"); 
-//     #pragma omp parallel for private(i,j,xy,dist2)
+    /* safety distance to avoid singularity at poles */
+    jmin = DPOLE;
+    jmax = NY-DPOLE;
+    
+    #pragma omp parallel for private(i,j,xy,dist2)
     for (i=0; i<NX; i++)
     {
-        if (i%100 == 0) printf("Initializing column %i of %i\n", i, NX);
-        for (j=0; j<NY; j++)
+        for (j=jmin; j<jmax; j++)
         {
             pscal = cos(wsphere[i*NY+j].phi - phi0 - phishift)*wsphere[i*NY+j].stheta*stheta;
             pscal += wsphere[i*NY+j].ctheta*ctheta;
@@ -1260,8 +1575,6 @@ void add_circular_wave_sphere(double amp, double phi0, double theta0, double phi
             
 	    if ((xy_in[i*NY+j])||(TWOSPEEDS)) 
                 phi[i*NY+j] += amp*INITIAL_AMP*exp(-dist2/INITIAL_VARIANCE)*cos(-sqrt(dist2)/INITIAL_WAVELENGTH);
-            else phi[i*NY+j] = 0.0;
-            psi[i*NY+j] = 0.0;
         }
     }
 }
@@ -1284,13 +1597,23 @@ void compute_light_angle_sphere(short int xy_in[NX*NY], t_wave wave[NX*NY], t_wa
     
     #pragma omp parallel for private(i,j)
     for (i=0; i<NX; i++)
-        for (j=0; j<NY; j++) if (xy_in[i*NY+j])
+    {
+        for (j=0; j<NY - DPOLE; j++) if (xy_in[i*NY+j])
         {
             r = 1.0 + RSCALE*(*wave[i*NY+j].p_zfield[movie]);
             if (r > RMAX) r = RMAX;
+            if (r < RMIN) r = RMIN;
             wsphere[i*NY+j].radius = r;
         }
-    
+        /* TODO ? Avoid artifacts due to singularity at north pole */
+        for (j=NY - DPOLE; j<NY; j++) if (xy_in[i*NY+j]) 
+        {
+            r = 1.0 + RSCALE*(*wave[i*NY+j].p_zfield[movie]);
+            if (r > RMAX) r = RMAX;
+            if (r < 1.0) r = 1.0;
+            wsphere[i*NY+j].radius = r;
+        }
+    }
     
     if (SHADE_WAVE)
     {
@@ -1523,9 +1846,12 @@ void draw_wave_sphere_ij(int i, int iplus, int j, int jplus, int jcolor, int mov
 /* draw wave at simulation grid point (i,j) */
 {
     int k, l;
-    short int draw;
+    short int draw, draw_bc=1;
     double xyz[3], ca;
     
+    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]);
+    
 //     if ((TWOSPEEDS)||(xy_in[i*NY+j]))
 //     if (wsphere[i*NY+j].indomain)
 //     if ((wsphere[i*NY+j].indomain)||(phi[i*NY+j] >= wsphere[i*NY+j].altitude + 0.001))
@@ -1542,22 +1868,25 @@ void draw_wave_sphere_ij(int i, int iplus, int j, int jplus, int jcolor, int mov
             glColor3f(wsphere[i*NY+j].r*ca, wsphere[i*NY+j].g*ca, wsphere[i*NY+j].b*ca);
         else glColor3f(COLOR_OUT_R*ca, COLOR_OUT_G*ca, COLOR_OUT_B*ca);
     }
-    glBegin(GL_TRIANGLE_FAN);
-    if (ij_to_sphere(i, j, *wave[i*NY+j].p_zfield[movie], wsphere, xyz))
-        draw_vertex_sphere(xyz);
-    if (ij_to_sphere(iplus, j, *wave[iplus*NY+j].p_zfield[movie], wsphere, xyz))
-        draw_vertex_sphere(xyz);
-    if (ij_to_sphere(iplus, jplus, *wave[iplus*NY+j+1].p_zfield[movie], wsphere, xyz))
-        draw_vertex_sphere(xyz);
-    if (ij_to_sphere(i, jplus, *wave[i*NY+j+1].p_zfield[movie], wsphere, xyz))
-        draw_vertex_sphere(xyz);
-    glEnd ();
+    if (draw_bc)
+    {
+        glBegin(GL_TRIANGLE_FAN);
+        if (ij_to_sphere(i, j, *wave[i*NY+j].p_zfield[movie], wsphere, xyz))
+            draw_vertex_sphere(xyz);
+        if (ij_to_sphere(iplus, j, *wave[iplus*NY+j].p_zfield[movie], wsphere, xyz))
+            draw_vertex_sphere(xyz);
+        if (ij_to_sphere(iplus, jplus, *wave[iplus*NY+j+1].p_zfield[movie], wsphere, xyz))
+            draw_vertex_sphere(xyz);
+        if (ij_to_sphere(i, jplus, *wave[i*NY+j+1].p_zfield[movie], wsphere, xyz))
+            draw_vertex_sphere(xyz);
+        glEnd ();
+    }
 }
 
 
 void draw_wave_sphere_3D(int movie, double phi[NX*NY], double psi[NX*NY], short int xy_in[NX*NY], t_wave wave[NX*NY], t_wave_sphere wsphere[NX*NY], int zplot, int cplot, int palette, int fade, double fade_value, int refresh)
 {
-    int i, j, imax, imin;
+    int i, j, imax, imin, jmax;
     double observer_angle, angle2;
     
     blank();
@@ -1581,6 +1910,9 @@ void draw_wave_sphere_3D(int movie, double phi[NX*NY], double psi[NX*NY], short
     if (imin >= NX-1) imin = NX-2;
     if (imax >= NX-1) imax = NX-2;
     
+//     jmax = NY-25;
+    jmax = NY-3;
+    
 //     printf("Angle = %.5lg, angle2 = %.5lg, imin = %i, imax = %i\n", observer_angle, angle2, imin, imax);
     
     if (observer[2] > 0.0)
@@ -1588,79 +1920,79 @@ void draw_wave_sphere_3D(int movie, double phi[NX*NY], double psi[NX*NY], short
         if (imin < imax)
         {
             for (i=imax; i>imin; i--)
-                for (j=0; j<NY-2; j++)
+                for (j=0; j<=jmax; j++)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
             for (i=imax+1; i<NX-1; i++)
-                for (j=0; j<NY-2; j++)
+                for (j=0; j<=jmax; j++)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
-            for (j=0; j<NY-2; j++)
+            for (j=0; j<=jmax; j++)
                 draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
             for (i=0; i<=imin; i++)
-                for (j=0; j<NY-2; j++)
+                for (j=0; j<=jmax; j++)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         }
         else
         {
             for (i=imax; i<imin; i++)
-                for (j=0; j<NY-2; j++)
+                for (j=0; j<=jmax; j++)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
             for (i=imax-1; i>=0; i--)
-                for (j=0; j<NY-2; j++)
+                for (j=0; j<=jmax; j++)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
-            for (j=0; j<NY-2; j++)
+            for (j=0; j<=jmax; j++)
                 draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
             for (i=NX-2; i>=imin; i--)
-                for (j=0; j<NY-2; j++)
+                for (j=0; j<=jmax; j++)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         }
     
         /* North pole */
-        for (i=0; i<NX-1; i++) for (j=NY-2; j<NY-1; j++)
-            draw_wave_sphere_ij(i, i+1, j-1, j, j, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        for (i=0; i<NX-1; i++) 
+            draw_wave_sphere_ij(i, i+1, NY-3, NY-2, NY-2, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
-        for (j=NY-2; j<NY-1; j++)
-            draw_wave_sphere_ij(NX-1, 0, j-1, j, NY-DPOLE, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
+        draw_wave_sphere_ij(NX-1, 0, NY-3, NY-2, NY-DPOLE, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
     }
     else
     {
         if (imin < imax)
         {
             for (i=imax; i>imin; i--)
-                for (j=NY-3; j>=0; j--)
+                for (j=jmax; j>=0; j--)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
             for (i=imax+1; i<NX-1; i++)
-                for (j=NY-3; j>=0; j--)
+                for (j=jmax; j>=0; j--)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
-            for (j=NY-3; j>=0; j--)
+            for (j=jmax; j>=0; j--)
                 draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
             for (i=0; i<=imin; i++)
-                for (j=NY-3; j>=0; j--)
+                for (j=jmax; j>=0; j--)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         }
         else
         {
             for (i=imax; i<imin; i++)
-                for (j=NY-3; j>=0; j--)
+                for (j=jmax; j>=0; j--)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
             for (i=imax-1; i>=0; i--)
-                for (j=NY-3; j>=0; j--)
+                for (j=jmax; j>=0; j--)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
-            for (j=NY-3; j>=0; j--)
+            for (j=jmax; j>=0; j--)
                 draw_wave_sphere_ij(NX-1, 0, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         
-            for (i=NX-2; i>=imin; i--)
-                for (j=NY-3; j>=0; j--)
+//             for (i=NX-2; i>=imin; i--)
+            for (i=NX-2; i>=imin-1; i--)
+                for (j=jmax; j>=0; j--)
                     draw_wave_sphere_ij(i, i+1, j, j+1, j+1, movie, phi, psi, xy_in, wave, wsphere, zplot, cplot, palette, fade, fade_value);
         }
     

sub_wave.c

@@ -2444,7 +2444,7 @@ int xy_in_arc(double x, double y, t_arc arc)
 int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_circle *circles)
 /* returns 1 if (x,y) represents a point in the billiard */
 {
-    double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height, ca, sa, l, ht;
+    double l2, r2, r2mu, omega, b, c, angle, z, x1, y1, x2, y2, y3, u, v, u1, v1, dx, dy, width, alpha, s, a, r, height, ca, sa, l, ht, xshift, zz[9][2], x5, x6, f, fp;
     int i, j, k, k1, k2, condition = 0, m;
     static int first = 1, nsides;
     static double h, hh, ra, rb, ll, salpha;
@@ -3025,6 +3025,82 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
             if ((x1 > MU)&&(x1 < MU + width)) return(1);
             return(0);
         }   
+        case (D_WAVEGUIDES_W):
+        {
+            x1 = vabs(x);
+            
+            /* upper fiber */
+            width = LAMBDA - 2.0*MU;
+            height = MU + width;
+            if (y > MU + width + height) return(0);
+            if (y >= height)
+            {
+                r = module2(x1, y-height);
+                if ((r > MU)&&(r < MU + width)) return(1);
+                if ((y < height + MU)&&(x1 > LAMBDA + MU)&&(x1 < LAMBDA + MU + width)) return(1);
+                return(0);
+            }
+            else
+            {
+                r = module2(x1-LAMBDA, y-height);
+                if ((r > MU)&&(r < MU + width)) return(1);
+            }
+            
+            /* lower fiber */
+            height = -MU - width;
+            width = LAMBDA - 2.0*MU_B;
+            if (y < -MU_B - width + height) return(0);
+            if (y >= +height)
+            {
+                r = module2(x1, y-height);
+                if ((r > MU_B)&&(r < MU_B + width)) return(1);
+                if ((y < height + MU_B)&&(x1 > LAMBDA + MU_B)&&(x1 < LAMBDA + MU_B + width)) return(1);
+                return(0);
+            }
+            else
+            {
+                r = module2(x1-LAMBDA, y-height);
+                if ((r > MU_B)&&(r < MU_B + width)) return(1);
+            }
+            
+            return(0);
+        }   
+        case (D_WAVEGUIDES_COUPLED):
+        {
+            x1 = vabs(x);
+            y1 = vabs(y);
+            
+            width = 0.03;
+            b = 1.0*(1.0 - 2.0*MU - width);
+            a = MU + 0.5*width - b;
+            x6 = x1*x1*x1;
+            x5 = x6*x1*x1;
+            x6 = x5*x1;
+            
+            f = a + b*(1.0 + 2.0*x6)/(1.0 + x6);
+            fp = 6.0*b*x5/((1.0 + x6)*(1.0 + x6));
+            
+            return(vabs(f - y1) < MU*sqrt(1.0 + fp*fp));
+        }
+        case (D_WAVEGUIDE_S):
+        {
+            a = 0.5*(1.0 - MU);
+            if (x < 0.0)
+            {
+                x = -x;
+                y = -y;
+            }
+            if (x < LAMBDA)
+            {
+                if (vabs(y - 2.0*a) < MU) return(1);
+                if (vabs(y) < MU) return(1);
+                if (vabs(y + 2.0*a) < MU) return(1);
+                return(0);
+            }
+            r = module2(x-LAMBDA, y-a);
+            return(vabs(r - a) < MU);
+            
+        }
         case (D_MAZE):
         {
             for (i=0; i<npolyrect; i++)
@@ -3137,6 +3213,125 @@ int xy_in_billiard_single_domain(double x, double y, int b_domain, int ncirc, t_
             if (x > 0) return ((x+a)*(x+a) + y*y > LAMBDA*LAMBDA);
             return ((x-a)*(x-a) + y*y > LAMBDA*LAMBDA);
         }
+        case (D_LENS_WALL):
+        {
+            if (vabs(y) > MU) return(1);
+            a = LAMBDA - 0.25;
+            if (x > 0) return ((x+a)*(x+a) + y*y > LAMBDA*LAMBDA);
+            return ((x-a)*(x-a) + y*y > LAMBDA*LAMBDA);
+        }
+        case (D_TWO_LENSES_WALL):
+        {
+            width = 0.2;
+            a = 1.9;
+            b = 1.9 - 2.0*LAMBDA + 2.0*width;
+            x1 = vabs(x);
+            if (module2(x1 - a, y) > LAMBDA) return(1);
+            if (module2(x1 - b, y) > LAMBDA) return(1);
+            return(0);
+        }
+        case (D_TWO_LENSES_OBSTACLE):
+        {
+            width = 0.2;
+            a = 1.9;
+            b = 1.9 - 2.0*LAMBDA + 2.0*width;
+            x1 = vabs(x);
+            if (module2(x1 - a, y) > LAMBDA) return(1);
+            if (module2(x1 - b, y) > LAMBDA) return(1);
+            return(0);
+        }
+        case (D_FRESNEL_ZONE_PLATE):
+        {
+            a = 6.25;
+            b = 0.125;
+            if (vabs(x + LAMBDA) >= MU) return(1);
+            return(cos(DPI*(a*y*y + b)) < 0.0);
+        }
+        case (D_FRESNEL_ZONE_PLATE_INV):
+        {
+            a = 6.25;
+            if (vabs(x + LAMBDA) >= MU) return(1);
+            return(cos(DPI*a*y*y) > 0.0);
+        }
+        case (D_LENS_ROTATED):
+        {
+            ca = cos(APOLY*PID);
+            sa = sin(APOLY*PID);
+            x1 = x*ca + y*sa;
+            y1 = -x*sa + y*ca;
+            if (vabs(y1) > MU) return(1);
+            a = LAMBDA - 0.25;
+            if (x1 > 0) return ((x1+a)*(x1+a) + y1*y1 > LAMBDA*LAMBDA);
+            return ((x1-a)*(x1-a) + y1*y1 > LAMBDA*LAMBDA);
+        }
+        case (D_LENS_CONCAVE):
+        {
+            width = 0.1;
+            x1 = vabs(x);
+            if (vabs(y) > MU) return(1);
+            return((x1 - LAMBDA - width)*(x1 - LAMBDA - width) + y*y < LAMBDA*LAMBDA);
+        }
+        case (D_LENS_CONVEX_CONCAVE):
+        {
+            if (vabs(y) > MU) return(1);
+            xshift = 0.6;
+            if (x > 0.0)
+            {
+                x1 = x - xshift;
+                a = LAMBDA - 0.25;
+                if (x1 > 0) return ((x1+a)*(x1+a) + y*y > LAMBDA*LAMBDA);
+                return ((x1-a)*(x1-a) + y*y > LAMBDA*LAMBDA);
+            }
+            else
+            {
+                width = 0.1;
+                x1 = vabs(x + xshift);
+                return((x1 - LAMBDA - width)*(x1 - LAMBDA - width) + y*y < LAMBDA*LAMBDA);
+            }
+        }
+        case (D_TREE):
+        {
+            /* upper triangle */
+            h = 1.5*LAMBDA;
+            r = 1.5*LAMBDA;
+            x1 = vabs(x);
+            y1 = y - h;
+            zz[0][0] = r*cos(PI/6.0);   zz[0][1] = -0.5*r;
+            zz[1][0] = 0.0;             zz[1][1] = r;
+            zz[2][0] = 0.0;             zz[2][1] = -0.5*r;
+            
+            if (y1 > 0.95*zz[1][1]) return(1);
+            
+            /* middle triangle */
+            h = 0.25*LAMBDA;
+            r = 2.0*LAMBDA;
+            y2 = y - h;
+            zz[3][0] = r*cos(PI/6.0);   zz[3][1] = -0.5*r;
+            zz[4][0] = 0.0;             zz[4][1] = r;
+            zz[5][0] = 0.0;             zz[5][1] = -0.5*r;
+            
+            /* bottom triangle */
+            h = -1.5*LAMBDA;
+            r = 3.0*LAMBDA;
+            y3 = y - h;
+            zz[6][0] = r*cos(PI/6.0);   zz[6][1] = -0.5*r;
+            zz[7][0] = 0.0;             zz[7][1] = r;
+            zz[8][0] = 0.0;             zz[8][1] = -0.5*r;
+            
+            if (xy_in_triangle(x1, y1, zz[0], zz[1], zz[2])) return(0);
+            if (xy_in_triangle(x1, y2, zz[3], zz[4], zz[5])) return(0);
+            if (xy_in_triangle(x1, y3, zz[6], zz[7], zz[8])) return(0);
+            
+            if (xy_in_triangle(0.9*x1, 0.9*y1, zz[0], zz[1], zz[2])) return(1);
+            if (xy_in_triangle(0.9*x1, 0.9*y2, zz[3], zz[4], zz[5])) return(1);
+            if (xy_in_triangle(0.9*x1, 0.9*y3, zz[6], zz[7], zz[8])) return(1);
+            
+            if (xy_in_triangle(0.8*x1, 0.8*y1, zz[0], zz[1], zz[2])) return(2);
+            if (xy_in_triangle(0.8*x1, 0.8*y2, zz[3], zz[4], zz[5])) return(2);
+            if (xy_in_triangle(0.8*x1, 0.8*y3, zz[6], zz[7], zz[8])) return(2);
+
+            return(1);
+        }
         case (D_MENGER):       
         {
             x1 = 0.5*(x+1.0);
@@ -3347,7 +3542,7 @@ void hex_transfo(double u, double v, double *x, double *y)
 
 void draw_billiard(int fade, double fade_value)      /* draws the billiard boundary */
 {
-    double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax, height, xmax, ca, sa;
+    double x0, y0, x, y, x1, y1, x2, y2, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, width, a, b, c, ymax, height, xmax, ca, sa, xshift, x5, x6, f, fp;
     int i, j, k, k1, k2, mr2, ntiles;
     static int first = 1, nsides;
     static double h, hh, sqr3, ll, salpha, arcangle;
@@ -4298,6 +4493,179 @@ void draw_billiard(int fade, double fade_value)      /* draws the billiard bound
             
             break;
         }
+        case (D_WAVEGUIDES_W):
+        {
+            /* upper waveguide */
+            width = LAMBDA - 2.0*MU;
+            height = MU + width;
+
+            draw_circle_arc(0.0, height, MU, 0.0, PI, NSEG);
+            draw_circle_arc(0.0, height, MU + width, 0.0, PI, NSEG);
+            
+            draw_circle_arc(LAMBDA, height, MU, PI, PI, NSEG);
+            draw_circle_arc(LAMBDA, height, MU + width, PI, PI, NSEG);
+            
+            draw_circle_arc(-LAMBDA, height, MU, PI, PI, NSEG);
+            draw_circle_arc(-LAMBDA, height, MU + width, PI, PI, NSEG);
+            
+            draw_line(LAMBDA + MU, height, LAMBDA + MU, height + MU);
+            draw_line(LAMBDA + MU, height + MU, LAMBDA + MU + width, height + MU);
+            draw_line(LAMBDA + MU + width, height + MU, LAMBDA + MU + width, height);
+            
+            draw_line(-LAMBDA - MU, height, -LAMBDA - MU, height + MU);
+            draw_line(-LAMBDA - MU, height + MU, -LAMBDA - MU - width, height + MU);
+            draw_line(-LAMBDA - MU - width, height + MU, -LAMBDA - MU - width, height);
+            
+            /* lower waveguide */
+            height = -MU - width;
+            width = LAMBDA - 2.0*MU_B;
+            
+            draw_circle_arc(0.0, height, MU_B, 0.0, PI, NSEG);
+            draw_circle_arc(0.0, height, MU_B + width, 0.0, PI, NSEG);
+            
+            draw_circle_arc(LAMBDA, height, MU_B, PI, PI, NSEG);
+            draw_circle_arc(LAMBDA, height, MU_B + width, PI, PI, NSEG);
+            
+            draw_circle_arc(-LAMBDA, height, MU_B, PI, PI, NSEG);
+            draw_circle_arc(-LAMBDA, height, MU_B + width, PI, PI, NSEG);
+            
+            draw_line(LAMBDA + MU_B, height, LAMBDA + MU_B, height + MU_B);
+            draw_line(LAMBDA + MU_B, height + MU_B, LAMBDA + MU_B + width, height + MU_B);
+            draw_line(LAMBDA + MU_B + width, height + MU_B, LAMBDA + MU_B + width, height);
+            
+            draw_line(-LAMBDA - MU_B, height, -LAMBDA - MU_B, height + MU_B);
+            draw_line(-LAMBDA - MU_B, height + MU_B, -LAMBDA - MU_B - width, height + MU_B);
+            draw_line(-LAMBDA - MU_B - width, height + MU_B, -LAMBDA - MU_B - width, height);
+            
+            break;
+        }
+        case (D_WAVEGUIDES_COUPLED):
+        {
+            width = 0.03;
+            b = 1.0*(1.0 - 2.0*MU - width);
+            a = MU + 0.5*width - b;
+            
+            if ((DRAW_WAVE_PROFILE)&&(VERTICAL_WAVE_PROFILE))
+                dx = (XMAX - XMIN - 0.06)/(double)NSEG;
+            else dx = (XMAX - XMIN + 0.1)/(double)NSEG;
+            
+            x1 = XMIN - 0.1;
+            y1 = 1.0;
+            for (i=0; i<NSEG; i++)
+            {
+                x2 = XMIN - 0.1 + (double)i*dx;
+                x6 = x2*x2*x2;
+                x5 = x6*x2*x2;
+                x6 = x5*x2;
+            
+                f = a + b*(1.0 + 2.0*x6)/(1.0 + x6);
+                fp = 6.0*b*x5/((1.0 + x6)*(1.0 + x6));
+                
+                y2 = f + MU*sqrt(1.0 + fp*fp);
+                
+                draw_line(x1, y1, x2, y2);
+                x1 = x2;
+                y1 = y2;
+            }
+            
+            x1 = XMIN - 0.1;
+            y1 = 1.0;
+            for (i=0; i<NSEG; i++)
+            {
+                x2 = XMIN - 0.1 + (double)i*dx;
+                x6 = x2*x2*x2;
+                x5 = x6*x2*x2;
+                x6 = x5*x2;
+            
+                f = a + b*(1.0 + 2.0*x6)/(1.0 + x6);
+                fp = 6.0*b*x5/((1.0 + x6)*(1.0 + x6));
+                
+                y2 = f - MU*sqrt(1.0 + fp*fp);
+                
+                draw_line(x1, y1, x2, y2);
+                x1 = x2;
+                y1 = y2;
+            }
+            
+            x1 = XMIN - 0.1;
+            y1 = -1.0;
+            for (i=0; i<NSEG; i++)
+            {
+                x2 = XMIN - 0.1 + (double)i*dx;
+                x6 = x2*x2*x2;
+                x5 = x6*x2*x2;
+                x6 = x5*x2;
+            
+                f = a + b*(1.0 + 2.0*x6)/(1.0 + x6);
+                fp = 6.0*b*x5/((1.0 + x6)*(1.0 + x6));
+                
+                y2 = -(f + MU*sqrt(1.0 + fp*fp));
+                
+                draw_line(x1, y1, x2, y2);
+                x1 = x2;
+                y1 = y2;
+            }
+            
+            x1 = XMIN - 0.1;
+            y1 = -1.0;
+            for (i=0; i<NSEG; i++)
+            {
+                x2 = XMIN - 0.1 + (double)i*dx;
+                x6 = x2*x2*x2;
+                x5 = x6*x2*x2;
+                x6 = x5*x2;
+            
+                f = a + b*(1.0 + 2.0*x6)/(1.0 + x6);
+                fp = 6.0*b*x5/((1.0 + x6)*(1.0 + x6));
+                
+                y2 = -f + MU*sqrt(1.0 + fp*fp);
+                
+                draw_line(x1, y1, x2, y2);
+                x1 = x2;
+                y1 = y2;
+            }
+            
+            x2 = LAMBDA;
+            x6 = x2*x2*x2;
+            x5 = x6*x2*x2;
+            x6 = x5*x2;
+            
+            f = a + b*(1.0 + 2.0*x6)/(1.0 + x6);
+            fp = 6.0*b*x5/((1.0 + x6)*(1.0 + x6));
+                
+            y1 = f - MU*sqrt(1.0 + fp*fp);
+            y2 = f + MU*sqrt(1.0 + fp*fp);
+            
+            draw_line(x2, YMIN, x2, -y2);
+            draw_line(x2, -y1, x2, y1);
+            draw_line(x2, y2, x2, YMAX);
+            draw_line(-x2, YMIN, -x2, -y2);
+            draw_line(-x2, -y1, -x2, y1);
+            draw_line(-x2, y2, -x2, YMAX);
+            
+            break;
+        }
+        case (D_WAVEGUIDE_S):
+        {
+            a = 0.5*(1.0 - MU);
+            
+            draw_circle_arc(LAMBDA, a, a-MU, -PID, PI, NSEG);
+            draw_circle_arc(LAMBDA, a, a+MU, -PID, PI, NSEG);
+            draw_circle_arc(-LAMBDA, -a, a-MU, PID, PI, NSEG);
+            draw_circle_arc(-LAMBDA, -a, a+MU, PID, PI, NSEG);
+            
+            draw_line(LAMBDA, 1.0, -LAMBDA, 1.0);
+            draw_line(-LAMBDA, 1.0, -LAMBDA, 2.0*a-MU);
+            draw_line( -LAMBDA, 2.0*a-MU, LAMBDA, 2.0*a-MU);
+            
+            draw_line(-LAMBDA, MU, LAMBDA, MU);
+            draw_line(-LAMBDA, -MU, LAMBDA, -MU);
+            
+            draw_line(-LAMBDA, -1.0, LAMBDA, -1.0);
+            draw_line(LAMBDA, -1.0, LAMBDA, -2.0*a+MU);
+            draw_line(LAMBDA, -2.0*a+MU, -LAMBDA, -2.0*a+MU);
+            break;
+        }
         case (D_CIRCLE_SEGMENT):
         {
             glLineWidth(BOUNDARY_WIDTH);
@@ -4560,6 +4928,167 @@ void draw_billiard(int fade, double fade_value)      /* draws the billiard bound
             draw_line(-ll, -MU, ll, -MU);
             break;
         }
+        case (D_LENS_WALL):
+        {
+            a = LAMBDA - 0.25;
+            width = 0.05;
+            if (first)
+            {
+                arcangle = asin(MU/LAMBDA);
+                ll = LAMBDA*cos(arcangle) - a;
+                first = 0;
+            }
+            
+            draw_circle_arc(-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_line(ll, MU, -ll, MU);
+            draw_circle_arc(a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_line(-ll, -MU, ll, -MU);
+            
+            draw_rectangle(-width, MU, width, YMAX + 1.0);
+            draw_rectangle(-width, YMIN-1.0, width, -MU);
+            break;
+        }
+        case (D_TWO_LENSES_WALL):
+        {
+            width = 0.2;
+            a = 1.9;
+            b = 1.9 - 2.0*LAMBDA + 2.0*width;
+            if (first)
+            {
+                arcangle = acos(1.0 - width/LAMBDA);
+                first = 0;
+            }
+            draw_circle_arc(a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(b, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(-b, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            
+            width = 0.05;
+            draw_rectangle(-width, MU, width, YMAX + 1.0);
+            draw_rectangle(-width, YMIN-1.0, width, -MU);
+            break;
+        }
+        case (D_FRESNEL_ZONE_PLATE):
+        {
+            a = 6.25;
+            b = 0.125;
+//             y = sqrt(1.0/(6.0*a));
+            y = sqrt((0.25 - b)/a);
+            draw_rectangle(-LAMBDA - MU, -y, -LAMBDA + MU, y);
+            for (i=1; i<(int)(a*a+1); i++)
+            {
+                x = sqrt((-0.25 - b + (double)i)/a);
+                y = sqrt((0.25 - b + (double)i)/a);
+                draw_rectangle(-LAMBDA - MU, x, -LAMBDA + MU, y);
+                draw_rectangle(-LAMBDA - MU, -x, -LAMBDA + MU, -y);
+            }
+            break;
+        }
+        case (D_FRESNEL_ZONE_PLATE_INV):
+        {
+            a = 6.25;
+            for (i=0; i<2*(int)a; i++)
+            {
+                x = sqrt((0.25 + (double)i)/a);
+                y = sqrt((-0.25 + (double)(i+1))/a);
+                draw_rectangle(-LAMBDA - MU, x, -LAMBDA + MU, y);
+                draw_rectangle(-LAMBDA - MU, -x, -LAMBDA + MU, -y);
+            }
+            break;
+        }
+        case (D_TWO_LENSES_OBSTACLE):
+        {
+            width = 0.2;
+            a = 1.9;
+            b = 1.9 - 2.0*LAMBDA + 2.0*width;
+            if (first)
+            {
+                arcangle = acos(1.0 - width/LAMBDA);
+                first = 0;
+            }
+            draw_circle_arc(a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(b, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(-b, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            
+            width = 0.05;
+            draw_rectangle(-width, -MU, width, MU);
+            break;
+        }
+        case (D_LENS_ROTATED):
+        {
+            ca = cos(APOLY*PID);
+            sa = sin(APOLY*PID);
+            a = LAMBDA - 0.25;
+            if (first)
+            {
+                arcangle = asin(MU/LAMBDA);
+                ll = LAMBDA*cos(arcangle) - a;
+                first = 0;
+            }
+            draw_circle_arc(-a*ca, -a*sa, LAMBDA, -arcangle + APOLY*PID, 2.0*arcangle, NSEG);
+            draw_line(ll*ca - MU*sa, ll*sa + MU*ca, -ll*ca - MU*sa, -ll*sa + MU*ca);
+            draw_circle_arc(a*ca, a*sa, LAMBDA, PI-arcangle + APOLY*PID, 2.0*arcangle, NSEG);
+            draw_line(-ll*ca + MU*sa, -ll*sa - MU*ca, ll*ca + MU*sa, ll*sa - MU*ca);
+            break;
+        }
+        case (D_LENS_CONCAVE):
+        {
+            arcangle = asin(MU/LAMBDA);
+            width = 0.1;
+            a = LAMBDA + width - sqrt(LAMBDA*LAMBDA - MU*MU);
+            draw_circle_arc(LAMBDA + width, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_line(a, MU, -a, MU);
+            draw_circle_arc(-LAMBDA - width, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_line(a, -MU, -a, -MU);
+            break;
+        }
+        case (D_LENS_CONVEX_CONCAVE):
+        {
+            xshift = 0.6;
+            
+            a = LAMBDA - 0.25;
+            if (first)
+            {
+                arcangle = asin(MU/LAMBDA);
+                ll = LAMBDA*cos(arcangle) - a;
+                first = 0;
+            }
+            
+            draw_circle_arc(xshift-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_line(xshift+ll, MU, xshift-ll, MU);
+            draw_circle_arc(xshift+a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_line(xshift-ll, -MU, xshift+ll, -MU);
+            
+            width = 0.1;
+            a = LAMBDA + width - sqrt(LAMBDA*LAMBDA - MU*MU);
+            draw_circle_arc(LAMBDA + width - xshift, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_line(a - xshift, MU, -a - xshift, MU);
+            draw_circle_arc(-LAMBDA - width - xshift, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_line(a - xshift, -MU, -a - xshift, -MU);
+            break;
+        }
+        case (D_TREE):
+        {
+//             h = LAMBDA;
+//             r = 2.0*LAMBDA;
+//             
+//             x1 = r*cos(PI/6.0);
+//             y1 = -h;
+//             y2 = r;
+//             
+//             x1 = x1/0.8;
+//             y1 = h + y1/0.8;
+//             x2 = 0.0;
+//             y2 = h + y2/0.8;
+//             
+//             draw_line(x1, y1, x2, y2);
+//             draw_line(x2, y2, -x1, y1);
+//             draw_line(-x1, y1, x1, y1);
+            
+            /* do nothing */
+            break;
+        }
         case (D_MENGER):
         {
             glLineWidth(3);
@@ -5211,6 +5740,12 @@ void draw_color_scheme_palette_fade(double x1, double y1, double x2, double y2,
                 value = min + 1.0*dy*(double)(j - jmin);
                 color_scheme_palette(COLOR_SCHEME, palette, 0.7*value, 1.0, 0, rgb);
                 break;
+            }
+             case (P_3D_LOG_MEAN_ENERGY):
+            {
+                value = LOG_SCALE*dy_e*(double)(j - jmin)*100.0/E_SCALE;
+                color_scheme_asym_palette(COLOR_SCHEME, palette, value, 1.0, 1, rgb);
+                break;
             }
             default:
             {
@@ -5893,13 +6428,13 @@ double oscillating_bc(int time, int j)
     }
 }
 
-void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
+void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double *tgamma_table[NX], double ior_angle)
 /* compute variable index of refraction */
 /* should be at some point merged with 3D version in suv_wave_3d.c */
 {
     int i, j, k, n, inlens, ncircles;
     double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1;
-    double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3];
+    double a, u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3];
     double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY];
     static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat;
     static int first = 1;
@@ -5935,7 +6470,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                         {
 //                             tc[i*NY+j] = COURANT;
                             tcc_table[i][j] = courant2;
-//                             tgamma[i*NY+j] = GAMMA;
+                            tgamma_table[i][j] = GAMMA;
                         }
                         else 
                         {
@@ -5943,7 +6478,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                             if (speed < 0.01) speed = 0.01;
                             tcc_table[i][j] = courant2*speed;
 //                             tc[i*NY+j] = COURANT*sqrt(speed);
-//                             tgamma[i*NY+j] = GAMMA;
+                            tgamma_table[i][j] = GAMMA;
                         }
                     }
                 }
@@ -5971,7 +6506,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                         {
 //                             tc[i*NY+j] = COURANT;
                             tcc_table[i][j] = courant2;
-//                             tgamma[i*NY+j] = GAMMA;
+                            tgamma_table[i][j] = GAMMA;
                         }
                         else 
                         {
@@ -5980,7 +6515,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                             else if (speed > 10.0) speed = 10.0;
                             tcc_table[i][j] = courant2*speed;
 //                             tc[i*NY+j] = COURANT*sqrt(speed);
-//                             tgamma[i*NY+j] = GAMMA;
+                            tgamma_table[i][j] = GAMMA;
                         }
                     }
                 }
@@ -5998,7 +6533,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                         else c = COURANT*(1.3 - 0.9*r2);
 //                         tc[i*NY+j] = c;
                         tcc_table[i][j] = c;
-//                         tgamma[i*NY+j] = GAMMA;
+                        tgamma_table[i][j] = GAMMA;
                     }
                 }
                 break;
@@ -6033,13 +6568,13 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                         else c = COURANT;
 //                         tc[i*NY+j] = c;
                         tcc_table[i][j] = c*c;
-//                         tgamma[i*NY+j] = GAMMA;
+                        tgamma_table[i][j] = GAMMA;
                     }
                     else
                     {
 //                         tc[i*NY+j] = 0.0;
                         tcc_table[i][j] = 0.0;
-//                         tgamma[i*NY+j] = 0.0;
+                        tgamma_table[i][j] = 0.0;
                     }
                 }
                 break;
@@ -6074,6 +6609,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
 //                         sum = tanh(sum);
 //                         printf("%.3lg\n", sum);
                         tcc_table[i][j] = COURANT*sum + COURANTB*(1.0-sum);
+                        tgamma_table[i][j] = GAMMA;
                     }
                 }
                 break;
@@ -6113,6 +6649,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                         sum = tanh(sum);
 //                         printf("%.3lg\n", sum);
                         tcc_table[i][j] = COURANT*sum + COURANTB*(1.0-sum);
+                        tgamma_table[i][j] = GAMMA;
                     }
                 }
                 break;
@@ -6160,6 +6697,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
 //                         sum = tanh(sum);
 //                         printf("%.3lg\n", sum);
                         tcc_table[i][j] = COURANT*sum + COURANTB*(1.0-sum);
+                        tgamma_table[i][j] = GAMMA;
                     }
                 }
                 break;
@@ -6203,6 +6741,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                             sum += exp(-r2/(sigma_stat));
                         }
                         tcc_table[i][j] = COURANT*sum + COURANTB*(1.0-sum);
+                        tgamma_table[i][j] = GAMMA;
                     }
                 }
                 break;
@@ -6238,6 +6777,7 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                         sum = tanh(sum);
 //                         printf("%.3lg\n", sum);
                         tcc_table[i][j] = COURANT*sum + COURANTB*(1.0-sum);
+                        tgamma_table[i][j] = GAMMA;
                     }
                 }
                 
@@ -6274,17 +6814,203 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                         sum = tanh(sum);
 //                         printf("%.3lg\n", sum);
                         tcc_table[i][j] = COURANT*sum + COURANTB*(1.0-sum);
+                        tgamma_table[i][j] = GAMMA;
                     }
                 }
                 
                 break;
-            }default:
+            }
+            case (IOR_LENS_WALL):
+            {
+                printf("Initializing IOR_LENS_WALL\n");
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if ((vabs(x) < 0.05)&&(vabs(y) > MU)) 
+                        {
+                            tcc_table[i][j] = 0.0;
+                            tgamma_table[i][j] = GAMMA;
+                        }
+                        else 
+                        {
+                            if (xy_in[i][j] != 0) 
+                            {
+                                tcc_table[i][j] = courant2;
+                                tgamma_table[i][j] = GAMMA;
+                            }
+                            else 
+                            {
+                                tcc_table[i][j] = courantb2;
+                                tgamma_table[i][j] = GAMMAB;
+                            }
+                        }
+                    }
+                }
+                
+                break;
+            }
+            case (IOR_LENS_OBSTACLE):
+            {
+                printf("Initializing IOR_LENS_OBSTACLE\n");
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if ((vabs(x) < 0.05)&&(vabs(y) < MU))
+                        {
+                            tcc_table[i][j] = 0.0;
+                            tgamma_table[i][j] = GAMMA;
+                        }
+                        else 
+                        {
+                            if (xy_in[i][j] != 0) 
+                            {
+                                tcc_table[i][j] = courant2;
+                                tgamma_table[i][j] = GAMMA;
+                            }
+                            else 
+                            {
+                                tcc_table[i][j] = courantb2;
+                                tgamma_table[i][j] = GAMMAB;
+                            }
+                        }
+                    }
+                }
+                break;
+            }
+            case (IOR_LENS_CONCAVE):
+            {
+                a = LAMBDA + 0.1 - sqrt(LAMBDA*LAMBDA - MU*MU); 
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if ((vabs(x) < a)&&(vabs(y) > MU))
+                        {
+                            tcc_table[i][j] = 0.0;
+                            tgamma_table[i][j] = GAMMA;
+                        }
+                        else 
+                        {
+                            if (xy_in[i][j] != 0) 
+                            {
+                                tcc_table[i][j] = courant2;
+                                tgamma_table[i][j] = GAMMA;
+                            }
+                            else 
+                            {
+                                tcc_table[i][j] = courantb2;
+                                tgamma_table[i][j] = GAMMAB;
+                            }
+                        }
+                    }
+                }
+                break;
+            }
+            case (IOR_LENS_CONVEX_CONCAVE):
+            {
+                a = LAMBDA + 0.1 - sqrt(LAMBDA*LAMBDA - MU*MU) + 0.6; 
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if ((vabs(x) < a)&&(vabs(y) > MU))
+                        {
+                            tcc_table[i][j] = 0.0;
+                            tgamma_table[i][j] = GAMMA;
+                        }
+                        else 
+                        {
+                            if (xy_in[i][j] != 0) 
+                            {
+                                tcc_table[i][j] = courant2;
+                                tgamma_table[i][j] = GAMMA;
+                            }
+                            else 
+                            {
+                                tcc_table[i][j] = courantb2;
+                                tgamma_table[i][j] = GAMMAB;
+                            }
+                        }
+                    }
+                }
+                break;
+            }
+            case (IOR_TREE):
+            {
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        switch (xy_in[i][j]) {
+                            case (0):
+                            {
+                                tcc_table[i][j] = courantb2;
+                                tgamma_table[i][j] = GAMMAB;
+                                break;
+                            }
+                            case (1):
+                            {
+                                tcc_table[i][j] = courant2;
+                                tgamma_table[i][j] = GAMMA;
+                                break;
+                            }
+                            case (2):
+                            {
+                                tcc_table[i][j] = 0.0;
+                                tgamma_table[i][j] = GAMMA;
+                                break;
+                            }
+                            case (3):
+                            {
+                                tcc_table[i][j] = courant2;
+                                tgamma_table[i][j] = GAMMA;
+                                break;
+                            }
+                        }
+                    }
+                }
+                break;
+            }
+            case (IOR_WAVE_GUIDES_COUPLED):
+            {
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if (xy_in[i][j]) 
+                        {
+                            tcc_table[i][j] = courant2;
+                            tgamma_table[i][j] = GAMMA;
+                        }
+                        else if (vabs(x) > LAMBDA)
+                        {
+                            tcc_table[i][j] = 0.0;
+                            tgamma_table[i][j] = 1.0;                            
+                        }
+                        else
+                        {
+                            tcc_table[i][j] = courantb2;
+                            tgamma_table[i][j] = GAMMAB;
+                        }
+                    }
+                }
+                break;
+            }
+            default:
             {
                 for (i=0; i<NX; i++){
                     for (j=0; j<NY; j++){
 //                         tc[i*NY+j] = COURANT;
                         tcc_table[i][j] = COURANT;
-//                         tgamma[i*NY+j] = GAMMA;
+                        tgamma_table[i][j] = GAMMA;
                     }
                 }
             }
@@ -6299,14 +7025,14 @@ void init_ior_2d(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
                 {
 //                     tc[i*NY+j] = COURANT;
                     tcc_table[i][j] = courant2;
-//                     if (xy_in[i*NY+j] == 1) tgamma[i*NY+j] = GAMMA;
-//                     else tgamma[i*NY+j] = GAMMAB;
+                    if (xy_in[i][j] == 1) tgamma_table[i][j] = GAMMA;
+                    else tgamma_table[i][j] = GAMMAB;
                 }
                 else if (TWOSPEEDS)
                 {
 //                     tc[i*NY+j] = COURANTB;
                     tcc_table[i][j] = courantb2;
-//                     tgamma[i*NY+j] = GAMMAB;
+                    tgamma_table[i][j] = GAMMAB;
                 }
             }
         }
@@ -6576,3 +7302,90 @@ void add_wave_packets(double *phi[NX], double *psi[NX], short int * xy_in[NX], t
     if (local) add_wave_packets_locally(phi, psi, packet, time, radius, add_period, type_envelope);
     else add_wave_packets_globally(phi, psi, xy_in, packet, time);
 }
+
+int old_source_schedule(int i)
+{
+    int mod;
+    
+    if (i < 200) return(0);
+    mod = i%(10*OSCILLATING_SOURCE_PERIOD);
+    if ((mod < 2*OSCILLATING_SOURCE_PERIOD)&&(rand()%3 < 2)) return(1);
+    return(0);
+}
+
+void old_init_input_signal()
+{
+    int i;
+ 
+    for (i=0; i<NSTEPS; i++) input_signal[i] = old_source_schedule(i);
+}
+
+
+int init_input_signal()
+{
+    int i, j, k, mod, ldash, ldot, lspace, linter, initial_time;
+    char text[200] = "dddd_dD_dDDd_dDDd_DdDD__Dd_d_dDD__DdDD_d_dD_dDd__ddDDD_DDDDD_ddDDD_ddddD", c;
+    
+    ldash = 12;
+    ldot = 3;
+    linter = 30;
+    lspace = 36;
+    initial_time = 200;
+        
+    for (i=0; i<initial_time; i++) input_signal[i] = 0;
+    
+    j = 0;
+    while (i < NSTEPS)
+    {
+        while (j < strlen(text))
+        {
+            c = text[j];
+            
+            switch(c) {
+                case ('D'):
+                {
+                    for (k=0; k<ldash; k++) 
+                    {
+                        input_signal[i] = 1;
+                        i++;
+                    }
+                    break;
+                }
+                case ('d'):
+                {
+                    for (k=0; k<ldot; k++) 
+                    {
+                        input_signal[i] = 1;
+                        i++;
+                    }
+                    break;
+                }
+                case ('_'):
+                {
+                    for (k=0; k<lspace; k++) 
+                    {
+                        input_signal[i] = 0;
+                        i++;
+                    }
+                    break;
+                }
+            }
+            
+            for (k=0; k<linter; k++) 
+            {
+                input_signal[i] = 0;
+                i++;
+            }
+            
+            j++;
+        }
+        
+        input_signal[i] = 0;
+        i++;
+    }
+    
+    for (i=0; i<NSTEPS; i++) printf("%i ", input_signal[i]);
+    
+    printf("\n\n j = %i, string length = %i\n", j, (int)strlen(text));
+}
+

sub_wave_comp.c

@@ -479,7 +479,7 @@ void init_polygon_config_comp(t_polygon polygons[NMAXCIRCLES])
 int xy_in_billiard_half(double x, double y, int domain, int pattern, int top)
 /* returns 1 if (x,y) represents a point in the billiard */
 {
-    double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy;
+    double l2, r2, r2mu, omega, c, angle, z, x1, y1, x2, y2, u, v, u1, v1, dx, dy, width, a, b;
     int i, j, k, k1, k2, condition, type;
 
     switch (domain) {
@@ -571,6 +571,26 @@ int xy_in_billiard_half(double x, double y, int domain, int pattern, int top)
                 else return(1);                
             }
         }
+        case (D_TWO_LENSES_WALL):
+        {
+            width = 0.2;
+            a = 1.9;
+            b = 1.9 - 2.0*LAMBDA + 2.0*width;
+            x1 = vabs(x);
+            if (module2(x1 - a, y) > LAMBDA) return(1);
+            if (module2(x1 - b, y) > LAMBDA) return(1);
+            return(0);
+        }
+        case (D_TWO_LENSES_OBSTACLE):
+        {
+            width = 0.2;
+            a = 1.9;
+            b = 1.9 - 2.0*LAMBDA + 2.0*width;
+            x1 = vabs(x);
+            if (module2(x1 - a, y) > LAMBDA) return(1);
+            if (module2(x1 - b, y) > LAMBDA) return(1);
+            return(0);
+        }
         default:
         {
             printf("Function ij_in_billiard not defined for this billiard \n");
@@ -599,19 +619,29 @@ int ij_in_billiard_comp(int i, int j)
 }
 
 
-void draw_billiard_half(int domain, int pattern, int top)      /* draws the billiard boundary */
+void draw_billiard_half(int domain, int pattern, int top, int fade, double fade_value)      
+/* draws the billiard boundary */
 /* two domain version implemented for D_CIRCLES */
 {
-    double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, signtop;
+    double x0, x, y, x1, y1, dx, dy, phi, r = 0.01, pos[2], pos1[2], alpha, dphi, omega, z, l, signtop, width, a, b, arcangle;
     int i, j, k, k1, k2, mr2;
+    static int first = 1;
     
     glEnable(GL_SCISSOR_TEST);
     if (top) glScissor(0.0, YMID, NX, YMID);
 //     if (top) glScissor(0.0, YMID, NX, YMAX);
     else glScissor(0.0, 0.0, NX, YMID);
 
-    if (BLACK) glColor3f(1.0, 1.0, 1.0);
-    else glColor3f(0.0, 0.0, 0.0);
+    if (fade)
+    {
+        if (BLACK) glColor3f(fade_value, fade_value, fade_value);
+        else glColor3f(1.0 - fade_value, 1.0 - fade_value, 1.0 - fade_value);        
+    }
+    else
+    {
+        if (BLACK) glColor3f(1.0, 1.0, 1.0);
+        else glColor3f(0.0, 0.0, 0.0);
+    }
     glLineWidth(5);
     
     if (top) signtop = 1.0;
@@ -734,6 +764,45 @@ void draw_billiard_half(int domain, int pattern, int top)      /* draws the bill
             glEnd();
             break;
         }
+        case (D_TWO_LENSES_WALL):
+        {
+            width = 0.2;
+            a = 1.9;
+            b = 1.9 - 2.0*LAMBDA + 2.0*width;
+            if (first)
+            {
+                arcangle = acos(1.0 - width/LAMBDA);
+                first = 0;
+            }
+            draw_circle_arc(a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(b, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(-b, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            
+            width = 0.05;
+            draw_rectangle(-width, MU, width, YMAX + 1.0);
+            draw_rectangle(-width, YMIN-1.0, width, -MUB);
+            break;
+        }
+        case (D_TWO_LENSES_OBSTACLE):
+        {
+            width = 0.2;
+            a = 1.9;
+            b = 1.9 - 2.0*LAMBDA + 2.0*width;
+            if (first)
+            {
+                arcangle = acos(1.0 - width/LAMBDA);
+                first = 0;
+            }
+            draw_circle_arc(a, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(b, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(-a, 0.0, LAMBDA, -arcangle, 2.0*arcangle, NSEG);
+            draw_circle_arc(-b, 0.0, LAMBDA, PI-arcangle, 2.0*arcangle, NSEG);
+            
+            width = 0.05;
+            draw_rectangle(-width, -MUB, width, MU);
+            break;
+        }
         case (D_MANDELBROT):
         {
             /* Do nothing */
@@ -764,10 +833,10 @@ void draw_billiard_half(int domain, int pattern, int top)      /* draws the bill
 }
 
 
-void draw_billiard_comp()      /* draws the billiard boundary */
+void draw_billiard_comp(int fade, double fade_value)      /* draws the billiard boundary */
 {
-    draw_billiard_half(B_DOMAIN, CIRCLE_PATTERN, 1);
-    draw_billiard_half(B_DOMAIN_B, CIRCLE_PATTERN_B, 0);
+    draw_billiard_half(B_DOMAIN, CIRCLE_PATTERN, 1, fade, fade_value);
+    draw_billiard_half(B_DOMAIN_B, CIRCLE_PATTERN_B, 0, fade, fade_value);
 }
 
 
@@ -978,14 +1047,15 @@ void draw_wave_comp_highres_palette(int size, double *phi[NX], double *psi[NX],
                         color_scheme(COLOR_SCHEME, LOG_SHIFT + LOG_SCALE*log(energy), scale, time, rgb);
                         break;
                     }
-//                     case (P_LOG_MEAN_ENERGY):
-//                     {
-//                         energy = compute_energy(phi, psi, xy_in, i, j);
-//                         if (energy == 0.0) energy = 1.0e-20;
-//                         total_energy[i][j] += energy;
-//                         color_scheme(COLOR_SCHEME, LOG_SCALE*log(total_energy[i][j]/(double)(time+1)), scale, time, rgb);
-//                         break;
-//                     }
+                    case (P_LOG_MEAN_ENERGY):
+                    {
+                        energy = compute_energy(phi, psi, xy_in, i, j);
+                        if (energy == 0.0) energy = 1.0e-20;
+                        total_energy[i][j] += energy;
+                        energy = LOG_SHIFT + LOG_SCALE*log(total_energy[i][j]/(double)(time+1));
+                        color_scheme_palette(COLOR_SCHEME, palette, energy, scale, time, rgb);
+                        break;
+                    }
                 }
                     
 //                 if (PLOT == P_AMPLITUDE)
@@ -1014,7 +1084,8 @@ void draw_wave_comp_highres_palette(int size, double *phi[NX], double *psi[NX],
     glEnd ();
     
     /* draw horizontal mid line */
-    glColor3f(1.0, 1.0, 1.0);
+    if (fade) glColor3f(fade_value, fade_value, fade_value);
+    else glColor3f(1.0, 1.0, 1.0);
     glBegin(GL_LINE_STRIP);
     xy_to_pos(XMIN, 0.5*(YMIN+YMAX), pos);
     glVertex2d(pos[0], pos[1]);
@@ -1202,3 +1273,109 @@ void print_energies(double energies[6], double top_energy, double bottom_energy)
     
 }
 
+void init_ior_2d_comp(short int *xy_in[NX], double *tcc_table[NX], double ior_angle)
+/* compute variable index of refraction */
+/* should be at some point merged with 3D version in suv_wave_3d.c */
+{
+    int i, j, k, n, inlens, ncircles;
+    double courant2 = COURANT*COURANT, courantb2 = COURANTB*COURANTB, lambda1, mu1;
+    double u, v, u1, x, y, xy[2], norm2, speed, r2, c, salpha, h, ll, ca, sa, x1, y1, dx, dy, sum, sigma, x0, y0, rgb[3];
+    double xc[NGRIDX*NGRIDY], yc[NGRIDX*NGRIDY], height[NGRIDX*NGRIDY];
+    static double xc_stat[NGRIDX*NGRIDY], yc_stat[NGRIDX*NGRIDY], sigma_stat;
+    static int first = 1;
+    t_circle circles[NMAXCIRCLES];
+    
+    rgb[0] = 1.0;
+    rgb[1] = 1.0;
+    rgb[2] = 1.0;
+    
+    if (VARIABLE_IOR)
+    {
+        switch (IOR) {
+            case (IOR_LENS_WALL):
+            {
+                printf("Initializing IOR_LENS_WALL\n");
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY/2; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if ((vabs(x) < 0.05)&&(vabs(y) > MUB)) tcc_table[i][j] = 0.0;
+                        else 
+                        {
+                            if (xy_in[i][j] != 0) tcc_table[i][j] = courant2;
+                            else tcc_table[i][j] = courantb2;
+                        }
+                    }
+                    for (j=NY/2; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if ((vabs(x) < 0.05)&&(vabs(y) > MU)) tcc_table[i][j] = 0.0;
+                        else 
+                        {
+                            if (xy_in[i][j] != 0) tcc_table[i][j] = courant2;
+                            else tcc_table[i][j] = courantb2;
+                        }
+                    }
+                }
+                
+                break;
+            }
+            case (IOR_LENS_OBSTACLE):
+            {
+                printf("Initializing IOR_LENS_OBSTACLE\n");
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY/2; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if ((vabs(x) < 0.05)&&(vabs(y) < MUB)) tcc_table[i][j] = 0.0;
+                        else 
+                        {
+                            if (xy_in[i][j] != 0) tcc_table[i][j] = courant2;
+                            else tcc_table[i][j] = courantb2;
+                        }
+                    }
+                    for (j=NY/2; j<NY; j++){
+                        ij_to_xy(i, j, xy);
+                        x = xy[0];
+                        y = xy[1];
+                        if ((vabs(x) < 0.05)&&(vabs(y) < MU)) tcc_table[i][j] = 0.0;
+                        else 
+                        {
+                            if (xy_in[i][j] != 0) tcc_table[i][j] = courant2;
+                            else tcc_table[i][j] = courantb2;
+                        }
+                    }
+                }
+                
+                break;
+            }
+            default:
+            {
+                for (i=0; i<NX; i++){
+                    for (j=0; j<NY; j++){
+                        tcc_table[i][j] = COURANT;
+                    }
+                }
+            }
+        }
+    }
+    else
+    {
+        #pragma omp parallel for private(i,j)
+        for (i=0; i<NX; i++){
+            for (j=0; j<NY; j++){
+                if (xy_in[i*NY+j] != 0)
+                {
+                    tcc_table[i][j] = courant2;
+                }
+                else if (TWOSPEEDS)
+                {
+                    tcc_table[i][j] = courantb2;
+                }
+            }
+        }
+    }
+}

wave_3d.c

@@ -291,6 +291,9 @@
 #define POTENTIAL 10
 #define POT_FACT 20.0
 #define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
+#define MU_B 1.0           /* parameter controlling the dimensions of domain */
+#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 */
 /* end of constants only used by sub_wave and sub_maze */
 
 

wave_billiard.c

@@ -49,11 +49,10 @@
 #define NO_EXTRA_BUFFER_SWAP 1    /* some OS require one less buffer swap when recording images */
 
 #define VARIABLE_IOR 0      /* set to 1 for a variable index of refraction */
-#define IOR 9               /* choice of index of refraction, see list in global_pdes.c */
+#define IOR 13               /* 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 */
 #define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
 
-
 /* General geometrical parameters */
 
 #define WINWIDTH 	1920  /* window width */
@@ -63,33 +62,16 @@
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
-#define YMIN -0.397916667
-#define YMAX 1.997916667	/* y interval for 9/16 aspect ratio */
-// #define YMIN -1.197916667
-// #define YMAX 1.197916667	/* y interval for 9/16 aspect ratio */
+#define YMIN -1.197916667
+#define YMAX 1.197916667	/* y interval for 9/16 aspect ratio */
 
 #define HIGHRES 1       /* set to 1 if resolution of grid is double that of displayed image */
 
-// #define WINWIDTH 	1280  /* window width */
-// #define WINHEIGHT 	720   /* window height */
-// 
-// // #define NX 640          /* number of grid points on x axis */
-// // #define NY 360          /* number of grid points on y axis */
-// // #define NX 1280          /* number of grid points on x axis */
-// // #define NY 720          /* number of grid points on y axis */
-// #define NX 2560          /* number of grid points on x axis */
-// #define NY 1440          /* number of grid points on y axis */
-// 
-// #define XMIN -2.0
-// #define XMAX 2.0	/* x interval  */
-// #define YMIN -1.125
-// #define YMAX 1.125	/* y interval for 9/16 aspect ratio */
-
 #define JULIA_SCALE 1.0 /* scaling for Julia sets */
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 20        /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 521        /* choice of domain shape, see list in global_pdes.c */
 
 #define CIRCLE_PATTERN 103   /* pattern of circles or polygons, see list in global_pdes.c */
 
@@ -102,9 +84,10 @@
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 1.0	    /* parameter controlling the dimensions of domain */
-#define MU 0.005            /* parameter controlling the dimensions of domain */
+#define MU 0.48             /* parameter controlling the dimensions of domain */
+#define MU_B 0.42           /* parameter controlling the dimensions of domain */
 #define NPOLY 6             /* number of sides of polygon */
-#define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
+#define APOLY -0.666666666666          /* angle by which to turn polygon, in units of Pi/2 */ 
 #define MDEPTH 6            /* depth of computation of Menger gasket */
 #define MRATIO 3            /* ratio defining Menger gasket */
 #define MANDELLEVEL 1000    /* iteration level for Mandelbrot set */
@@ -112,7 +95,6 @@
 #define FOCI 1              /* set to 1 to draw focal points of ellipse */
 #define NGRIDX 60           /* number of grid point for grid of disks */
 #define NGRIDY 10           /* number of grid point for grid of disks */
-// #define NGRIDY 18            /* number of grid point for grid of disks */
 
 #define X_SHOOTER -0.2
 #define Y_SHOOTER -0.6
@@ -130,38 +112,38 @@
 
 /* Physical parameters of wave equation */
 
-#define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
-#define OSCILLATE_LEFT 1     /* set to 1 to add oscilating boundary condition on the left */
+#define 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 0  /* oscillation schedule, see list in global_pdes.c */
 
-#define OMEGA 0.024        /* frequency of periodic excitation */
+#define OMEGA 0.01         /* frequency of periodic excitation */
 #define AMPLITUDE 1.0      /* amplitude of periodic excitation */ 
 #define ACHIRP 0.25        /* acceleration coefficient in chirp */
 #define DAMPING 0.0        /* damping of periodic excitation */
-#define COURANT 0.1        /* Courant number */
-#define COURANTB 0.01      /* Courant number in medium B */
+#define COURANT 0.06       /* Courant number */
+#define COURANTB 0.18       /* Courant number in medium B */
 #define GAMMA 0.0          /* damping factor in wave equation */
-#define GAMMAB 0.0         /* damping factor in wave equation */
+#define GAMMAB 0.007            /* damping factor in wave equation */
 #define GAMMA_SIDES 1.0e-4      /* damping factor on boundary */
 #define GAMMA_TOPBOT 1.0e-7     /* damping factor on boundary */
 #define KAPPA 0.0           /* "elasticity" term enforcing oscillations */
 #define KAPPA_SIDES 5.0e-4  /* "elasticity" term on absorbing boundary */
 #define KAPPA_TOPBOT 0.0    /* "elasticity" term on absorbing boundary */
-#define OSCIL_LEFT_YSHIFT -400.0   /* y-dependence of left oscillation (for non-horizontal waves) */
+#define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
 /* The physical damping coefficient is given by GAMMA/(DT)^2 */
 /* Increasing COURANT speeds up the simulation, but decreases accuracy */
 /* For similar wave forms, COURANT^2*GAMMA should be kept constant */
 
-#define ADD_OSCILLATING_SOURCE 0        /* set to 1 to add an oscillating wave source */
-#define OSCILLATING_SOURCE_PERIOD 30    /* period of oscillating source */
+#define ADD_OSCILLATING_SOURCE 1        /* set to 1 to add an oscillating wave source */
+#define OSCILLATING_SOURCE_PERIOD 8     /* period of oscillating source */
 #define ALTERNATE_OSCILLATING_SOURCE 1  /* set to 1 to alternate sign of oscillating source */
 
 #define ADD_WAVE_PACKET_SOURCES 0       /* set to 1 to add several sources emitting wave packets */
-#define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
-#define N_WAVE_PACKETS 15               /* number of wave packets */
-#define WAVE_PACKET_RADIUS 20            /* radius of wave packets */
+#define WAVE_PACKET_SOURCE_TYPE 2       /* type of wave packet sources */
+#define N_WAVE_PACKETS 2                /* number of wave packets */
+#define WAVE_PACKET_RADIUS 25            /* radius of wave packets */
 
 /* Boundary conditions, see list in global_pdes.c  */
 
@@ -169,11 +151,10 @@
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 3600       /* number of frames of movie */
-// #define NSTEPS 500       /* number of frames of movie */
-#define NVID 7            /* number of iterations between images displayed on screen */
+#define NSTEPS 2700       /* number of frames of movie */
+#define NVID 20           /* number of iterations between images displayed on screen */
 #define NSEG 1000         /* number of segments of boundary */
-#define INITIAL_TIME 700      /* time after which to start saving frames */
+#define INITIAL_TIME 0      /* time after which to start saving frames */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
 #define PRINT_SPEED 0       /* print speed of moving source */
 #define PRINT_FREQUENCY 0       /* print frequency (for phased array) */
@@ -182,29 +163,26 @@
 #define PSLEEP 1         /* sleep time during pause */
 #define SLEEP1  1        /* initial sleeping time */
 #define SLEEP2  1        /* final sleeping time */
-#define MID_FRAMES 20    /* number of still frames between parts of two-part movie */
-#define END_FRAMES 100    /* number of still frames at end of movie */
+#define MID_FRAMES 100    /* number of still frames between parts of two-part movie */
+#define END_FRAMES 300    /* number of still frames at end of movie */
 #define FADE 1           /* set to 1 to fade at end of movie */
 
 /* Parameters of initial condition */
 
-#define INITIAL_AMP 2.0            /* amplitude of initial condition */
-// #define INITIAL_VARIANCE 0.000015    /* variance of initial condition */
-#define INITIAL_VARIANCE 0.000025    /* variance of initial condition */
-#define INITIAL_WAVELENGTH  0.05  /* wavelength of initial condition */
+#define INITIAL_AMP 0.5            /* amplitude of initial condition */
+#define INITIAL_VARIANCE 0.000012   /* variance of initial condition */
+#define INITIAL_WAVELENGTH  0.05   /* wavelength of initial condition */
 
 /* Plot type, see list in global_pdes.c  */
 
 #define PLOT 0
-// #define PLOT 7
 
 #define PLOT_B 5        /* plot type for second movie */
 
 /* Color schemes */
 
-#define COLOR_PALETTE 18     /* Color palette, see list in global_pdes.c  */
-// #define COLOR_PALETTE 17      /* Color palette, see list in global_pdes.c  */
-#define COLOR_PALETTE_B 12     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 17     /* 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 */
 
@@ -215,9 +193,9 @@
 #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 E_SCALE 50.0     /* scaling factor for energy representation */
+#define E_SCALE 40.0      /* scaling factor for energy representation */
 #define LOG_SCALE 1.0     /* scaling factor for energy log representation */
-#define LOG_SHIFT 1.5     /* shift of colors on log scale */
+#define LOG_SHIFT 1.25     /* shift of colors on log scale */
 #define FLUX_SCALE 5.0e3    /* scaling factor for enegy flux represtnation */
 #define RESCALE_COLOR_IN_CENTER 0   /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
 
@@ -229,14 +207,15 @@
 #define HUEAMP -180.0      /* amplitude of variation of hue for color scheme C_HUE */
 
 #define DRAW_COLOR_SCHEME 1    /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 1.0     /* scale of color scheme bar */
+#define COLORBAR_RANGE 1.5     /* scale of color scheme bar */
 #define COLORBAR_RANGE_B 0.4   /* 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 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 SAVE_TIME_SERIES 0      /* set to 1 to save wave time series at a point */
 
 #define NXMAZE 8      /* width of maze */
@@ -275,7 +254,7 @@ double courant2, courantb2;  /* Courant parameters squared */
 // void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], double *psi_out[NX], 
 //                       short int *xy_in[NX])
 void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX], 
-                      short int *xy_in[NX], double *tcc[NX])
+                      short int *xy_in[NX], double *tcc[NX], double *tgamma[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 /* this version of the function has been rewritten in order to minimize the number of if-branches */
@@ -283,7 +262,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
     int i, j, iplus, iminus, jplus, jminus;
     double delta, x, y, c, cc, gamma, tb_shift;
     static long time = 0;
-    static double tc[NX][NY], tgamma[NX][NY];
+    static double tc[NX][NY];
     static short int first = 1;
     
     time++;
@@ -299,15 +278,21 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
                 if (xy_in[i][j] != 0)
                 {
                     tc[i][j] = COURANT;
-                    if (!VARIABLE_IOR) tcc[i][j] = courant2;
-                    if (xy_in[i][j] == 1) tgamma[i][j] = GAMMA;
-                    else tgamma[i][j] = GAMMAB;
+                    if (!VARIABLE_IOR) 
+                    {
+                        tcc[i][j] = courant2;
+                        if (xy_in[i][j] == 1) tgamma[i][j] = GAMMA;
+                        else tgamma[i][j] = GAMMAB;
+                    }
                 }
                 else if (TWOSPEEDS)
                 {
                     tc[i][j] = COURANTB;
-                    tcc[i][j] = courantb2;
-                    tgamma[i][j] = GAMMAB;
+                    if (!VARIABLE_IOR) 
+                    {
+                        tcc[i][j] = courantb2;
+                        tgamma[i][j] = GAMMAB;
+                    }
                 }
             }
         }
@@ -543,7 +528,7 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 }
 
 
-void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX])
+void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX], double *tgamma_table[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 {
@@ -553,11 +538,11 @@ void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *x
     // At the beginning w[t] is saved in phi, w[t-1] in psi and tmp is space
     // for the next wave state w[t+1]. Take w[t] and w[t-1] to calculate the
     // next wave state. Write this new state in temp
-    evolve_wave_half(phi, psi, tmp, xy_in, tcc_table);
+    evolve_wave_half(phi, psi, tmp, xy_in, tcc_table, tgamma_table);
     // now w[t] is saved in tmp, w[t-1] in phi and the result is written to psi
-    evolve_wave_half(tmp, phi, psi, xy_in, tcc_table);
+    evolve_wave_half(tmp, phi, psi, xy_in, tcc_table, tgamma_table);
     // now w[t] is saved in psi, w[t-1] in tmp and the result is written to phi
-    evolve_wave_half(psi, tmp, phi, xy_in, tcc_table);
+    evolve_wave_half(psi, tmp, phi, xy_in, tcc_table, tgamma_table);
     // now w[t] is saved in phi, w[t-1] in psi and tmp is free again to take
     // the new wave state w[t+1] in the next call to this function, thus
     // matching the given parameter names again
@@ -584,11 +569,10 @@ void draw_color_bar_palette(int plot, double range, int palette, int circular, i
         draw_color_scheme_palette_fade(XMAX - 1.5*width, YMIN + 0.1, XMAX - 0.5*width, YMAX - 0.1, plot, -range, range, palette, fade, fade_value);
 }
 
-
 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, ior_angle = 0.0, omega, phase_shift; 
-    double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX], *total_flux, *tcc_table[NX];
+    double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *color_scale[NX], *total_flux, *tcc_table[NX], *tgamma_table[NX];
     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;
     static int counter = 0;
@@ -612,6 +596,7 @@ void animation()
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
         color_scale[i] = (double *)malloc(NY*sizeof(double));
         tcc_table[i] = (double *)malloc(NX*sizeof(double));
+        tgamma_table[i] = (double *)malloc(NX*sizeof(double));
     }
     
     if (MEAN_FLUX) total_flux = (double *)malloc(4*NX*NY*sizeof(double));
@@ -639,6 +624,8 @@ void animation()
     printf("Rotated rectangles and arcs initialized\n");
     printf("%i rotated rectangles, %i arcs\n", npolyrect_rot, npolyarc);
     
+    if (DRAW_WAVE_TIMESERIES) init_input_signal();
+    
     courant2 = COURANT*COURANT;
     courantb2 = COURANTB*COURANTB;
     c = COURANT*(XMAX - XMIN)/(double)NX;
@@ -667,9 +654,7 @@ void animation()
     if (MEAN_FLUX) 
         for (i=0; i<4*NX*NY; i++)
             total_flux[i] = 0.0;
-    
-    if (VARIABLE_IOR) init_ior_2d(xy_in, tcc_table, ior_angle);
-    
+
     ratio = (XMAX - XMIN)/8.4;  /* for Tokarsky billiard */
     
 //     isospectral_initial_point(0.2, 0.0, startleft, startright);    /* for isospectral billiards */
@@ -682,8 +667,10 @@ void animation()
 //     printf("xleft = (%.3f, %.3f) xright = (%.3f, %.3f)\n", xin_left, yin_left, xin_right, yin_right);
     
     init_wave_flat(phi, psi, xy_in);
+    
+    if (VARIABLE_IOR) init_ior_2d(xy_in, tcc_table, tgamma_table, ior_angle);
 
-//     init_circular_wave(-0.5, 0.0, phi, psi, xy_in);
+//     init_circular_wave(-1.5, 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);
@@ -773,8 +760,7 @@ void animation()
         else draw_wave_epalette(phi, psi, total_energy, total_flux, color_scale, xy_in, scale, i, PLOT, COLOR_PALETTE, 0, 1.0);
         for (j=0; j<NVID; j++) 
         {
-//             evolve_wave(phi, psi, phi_tmp, psi_tmp, xy_in);
-            evolve_wave(phi, psi, tmp, xy_in, tcc_table);
+            evolve_wave(phi, psi, tmp, xy_in, tcc_table, tgamma_table);
             if (SAVE_TIME_SERIES)
             {
                 wave_value = (long int)(phi[sample_left[0]][sample_left[1]]*1.0e16);
@@ -795,8 +781,9 @@ void animation()
         if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
         {
             if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
-            add_circular_wave(sign, -0.5, 0.0, phi, psi, xy_in);
-
+            add_circular_wave(sign, -1.5, 0.8, phi, psi, xy_in);
+            add_circular_wave(sign, -1.5, -0.35, phi, psi, xy_in);
+        
 //             p = phased_array_schedule(i);
             
 //             phase_shift = 0.02 + 0.06*(double)i/(double)NSTEPS;
@@ -873,14 +860,14 @@ void animation()
 //             speed = speed/c;
 //             speed = 120.0*speed/((double)NVID*COURANT);
         }
-        if (ADD_WAVE_PACKET_SOURCES) add_wave_packets(phi, psi, xy_in, packet, i, WAVE_PACKET_RADIUS, 0, 10, 1);
+        if (ADD_WAVE_PACKET_SOURCES) add_wave_packets(phi, psi, xy_in, packet, i, WAVE_PACKET_RADIUS, 1, 4, 2);
         if (PRINT_SPEED) print_speed(speed, 0, 1.0);
         if (PRINT_FREQUENCY) print_frequency(phase_shift, 0, 1.0);
         if ((VARIABLE_IOR)&&(REFRESH_IOR)&&(i%3 == 0))
         {
             ior_angle = ior_angle_schedule(i);
             printf("IOR angle = %.5lg\n", ior_angle); 
-            init_ior_2d(xy_in, tcc_table, ior_angle);
+            init_ior_2d(xy_in, tcc_table, tgamma_table, ior_angle);
             printf("speed = %.5lg\n", tcc_table[3*NX/4][NY/2]);
         }
 
@@ -987,6 +974,7 @@ void animation()
         free(xy_in[i]);
         free(color_scale[i]);
         free(tcc_table[i]);
+        free(tgamma_table[i]);
     }
     
     if (MEAN_FLUX) free(total_flux);

wave_common.c

@@ -894,8 +894,7 @@ double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_
         case (P_LOG_ENERGY):
         {
             energy = compute_energy(phi, psi, xy_in, i, j);
-//                       energy = LOG_SHIFT + LOG_SCALE*log(energy);
-//                         if (energy < 0.0) energy = 0.0;
+            if (energy < 1.0e-20) energy = 1.0e-20;
             value = LOG_SHIFT + LOG_SCALE*log(energy);
             color_scheme_palette(COLOR_SCHEME, palette, value, scale, time, rgb);
             break;
@@ -947,8 +946,8 @@ double wave_value(int i, int j, double *phi[NX], double *psi[NX], double *total_
 }
 
 
-void draw_wave_profile(double *values, int size, int fade, double fade_value)
-/* draw a profile of the wave, if option DRAW_WAVE_PROFILE is active */
+void draw_wave_profile_horizontal(double *values, int size, int fade, double fade_value)
+/* draw a horizontal profile of the wave, if option DRAW_WAVE_PROFILE is active */
 {
     int i, k;
     double vmin, vmax, deltav, a, b, y;
@@ -976,6 +975,7 @@ void draw_wave_profile(double *values, int size, int fade, double fade_value)
         if (values[i*NY+jmin] > vmax) vmax = values[i*NY+jmin];
         if (values[i*NY+jmin] < vmin) vmin = values[i*NY+jmin];        
     }
+//     if (vmax <= vmin) vmax = vmin + 0.01;
     if ((vmin < 0.0)&&(vmax > 0.0))
     {
         if (vmax > -vmin) vmin = -vmax;
@@ -983,7 +983,7 @@ void draw_wave_profile(double *values, int size, int fade, double fade_value)
     }
 //     printf("vmin = %.3lg, vmax = %.3lg\n", vmin, vmax);
     deltav = vmax-vmin;
-    if (deltav == 0.0) deltav = 0.01;
+    if (deltav <= 0.0) deltav = 0.01;
     a = deltaj/deltav;
     b = ymin - a*vmin;
 
@@ -1008,6 +1008,249 @@ void draw_wave_profile(double *values, int size, int fade, double fade_value)
 }
 
 
+void draw_wave_profile_vertical(double *values, int size, int fade, double fade_value)
+/* draw a vertical profile of the wave, if option DRAW_WAVE_PROFILE is active */
+{
+    int j, k;
+    double vmin, vmax, deltav, a, b, x;
+    static int imin, imax, jmin, jmax, imid, d, first = 1;
+    static double deltai, xmin;
+    
+    if (first)
+    {        
+        jmin = 50;
+        jmax = NY - 50;
+        imin = NX - 150;
+        imax = NX - 50;
+        imid = (imin + imax)/2;
+        imid -= (imid%size);
+        d = (imax - imin)/10 + 1;
+        deltai = (double)(imax - imin - 2*d); 
+        xmin = (double)(imin + d);
+        first = 0;
+    }
+    
+    vmin = 1.0e10;
+    vmax = -vmin;
+    for (j=jmin; j<jmax; j+=size)
+    {
+        if (values[imin*NY+j] > vmax) vmax = values[imin*NY+j];
+        if (values[imin*NY+j] < vmin) vmin = values[imin*NY+j];        
+    }
+//     if ((vmin < 0.0)&&(vmax > 0.0))
+    {
+        if (vmax > -vmin) vmin = -vmax;
+        else if (vmax < -vmin) vmax = -vmin;
+    }
+//     else if (vmax > 0.0) vmin = 0.0;
+//     else if (vmin < 0.0) vmax = 0.0;
+//     printf("vmin = %.3lg, vmax = %.3lg\n", vmin, vmax);
+    deltav = vmax-vmin;
+    if (deltav <= 0.0) deltav = 0.01;
+    a = deltai/deltav;
+    b = xmin - a*vmin;
+
+    erase_area_ij(imin, jmin, imax, jmax);
+    
+    if (fade) glColor3f(fade_value, fade_value, fade_value);
+    else glColor3f(1.0, 1.0, 1.0);
+    glBegin(GL_LINE_STRIP);
+    for (j=jmin; j<jmax; j+=size)
+    {
+        x = a*values[imin*NY+j] + b;
+        glVertex2d(x, (double)j);
+    }
+    glEnd();
+    
+    glBegin(GL_LINE_LOOP);
+    glVertex2i(imin, jmin);
+    glVertex2i(imax, jmin);
+    glVertex2i(imax, jmax);
+    glVertex2i(imin, jmax);
+    glEnd();
+}
+
+
+void draw_wave_profile(double *values, int size, int fade, double fade_value)
+/* draw a profile of the wave, if option DRAW_WAVE_PROFILE is active */
+{
+    if (VERTICAL_WAVE_PROFILE) draw_wave_profile_vertical(values, size, fade, fade_value);
+    else draw_wave_profile_horizontal(values, size, fade, fade_value);
+}
+
+void draw_exit_timeseries(double *values, int size, int fade, double fade_value)
+/* draw a profile of the wave, if option DRAW_WAVE_TIMESERIES is active */
+{
+    int t, t1, s, padding = 50, nvals = 200;
+    double a, b, deltav, x, y, wmin, wmax;
+    static int ij[2], ij_test[2], imin, imax, jmin, jmax, jmid, itest, jtest, counter = 0, first  = 1, window = 3;
+    static double timeseries[200], average[200], vmin, vmax, ymin, deltaj, deltai;
+    
+    if (first)
+    {
+        xy_to_ij(LAMBDA, -1.0 + MU, ij);
+//         xy_to_ij(LAMBDA - 0.1*MU, -1.0 + MU, ij_test);
+        xy_to_ij(LAMBDA, -1.0 + MU, ij_test);
+        imin = ij[0];
+        imax = NX - padding;
+        jmin = padding;
+        jmax = 2*ij[1] - jmin;
+        jmid = (jmin + jmax)/2;
+        ymin = (double)(jmin);
+        itest = ij_test[0];
+        jtest = ij_test[1];
+        deltai = (double)(imax - imin)/(double)nvals;
+        deltaj = (double)(jmax - jmin); 
+        vmin = 1.0e10;
+        vmax = -vmin;
+        for (t=0; t<nvals; t++) timeseries[t] = 0.0;
+        first = 0;
+    }
+    
+    timeseries[counter] = values[itest*NY + jtest] + values[(itest-1)*NY + jtest] + values[(itest)*NY + jtest + 1] + values[(itest-1)*NY + jtest + 1];
+    counter++;
+    if (counter == nvals) counter = 0;
+    
+//     for (t=0; t<nvals; t++) printf("val[%i] = %.3lg\n", t, timeseries[t]);
+        
+    for (t=0; t<nvals; t++)
+    {
+        if (timeseries[t] > vmax) vmax = timeseries[t];
+        if (timeseries[t] < vmin) vmin = timeseries[t];        
+    }
+    if (vmax > -vmin) vmin = -vmax;
+    else if (vmax < -vmin) vmax = -vmin;
+    
+    printf("vmin = %.3lg, vmax = %.3lg\n", vmin, vmax);
+    deltav = vmax-vmin;
+    if (deltav <= 0.0) deltav = 0.01;
+    a = deltaj/deltav;
+    b = ymin - a*vmin;
+    
+    /* compute average */
+    for (t=window; t<nvals-window; t++)
+    {
+        average[t] = 0.0;
+        for (s=-window; s<window; s++) 
+        {
+            t1 = counter - t + s;
+            if (t1 < 0) t1 += nvals;
+            if (t1 > nvals-1) t1 -= nvals;
+            average[t] += timeseries[t1]*timeseries[t1];
+        }
+        average[t] = sqrt(average[t]*0.5/(double)window);
+    }
+        
+    erase_area_ij(imin, jmin, imax, jmax);
+    
+    if (fade) glColor3f(fade_value, fade_value, fade_value);
+    else glColor3f(1.0, 1.0, 1.0);
+    
+    glBegin(GL_LINE_STRIP);
+    for (t=0; t<nvals; t++)
+    {
+        t1 = counter - t;
+        if (t1 < 0) t1 += nvals;
+        x = (double)imin + deltai*(double)t;
+        y = a*timeseries[t1] + b;
+        glVertex2d(x, y);
+    }
+    glEnd();
+    
+    /* draw average */
+    if (fade) glColor3f(fade_value, 0.0, 0.0);
+    else glColor3f(1.0, 0.0, 0.0);
+    glBegin(GL_LINE_STRIP);
+    for (t1=window; t1<nvals-window; t1++)
+    {
+        x = (double)imin + deltai*(double)t1;
+        y = a*average[t1] + b;
+        glVertex2d(x, y);
+    }
+    glEnd();
+    glBegin(GL_LINE_STRIP);
+    for (t1=window; t1<nvals-window; t1++)
+    {
+        x = (double)imin + deltai*(double)t1;
+        y = -a*average[t1] + b;
+        glVertex2d(x, y);
+    }
+    glEnd();
+    
+    if (fade) glColor3f(fade_value, fade_value, fade_value);
+    else glColor3f(1.0, 1.0, 1.0);
+    
+    glBegin(GL_LINE_LOOP);
+    glVertex2i(imin, jmin);
+    glVertex2i(imax, jmin);
+    glVertex2i(imax, jmax);
+    glVertex2i(imin, jmax);
+    glEnd();
+}
+
+
+void draw_entrance_timeseries(double *values, int size, int fade, double fade_value)
+/* draw a profile of the wave, if option DRAW_WAVE_TIMESERIES is active */
+{
+    int t, t1, padding = 50, nvals = 200;
+    double vmin, vmax, deltav, x, y;
+    static int ij[2], imin, imax, jmin, jmax, jmid, itest, jtest, counter = 0, first  = 1, 
+    time = 0;
+    static double ymin, deltaj, deltai, a, b;
+    
+    if (first)
+    {
+        xy_to_ij(-LAMBDA, 1.0 - MU, ij);
+        imin = padding;
+        imax = ij[0];
+        jmax = NY - padding;
+        jmin = 2*ij[1] - jmax;
+        jmid = (jmin + jmax)/2;
+        ymin = (double)(jmin);
+        deltai = (double)(imax - imin)/(double)nvals;
+        deltaj = (double)(jmax - jmin); 
+        b = (double)jmin + 0.25*deltaj;
+        a = 0.5*deltaj;
+        
+        first = 0;
+    }
+    
+    counter++;
+        
+    erase_area_ij(imin, jmin, imax, jmax);
+    
+    if (fade) glColor3f(fade_value, fade_value, fade_value);
+    else glColor3f(1.0, 1.0, 1.0);
+    
+//     glEnable(GL_SCISSOR_TEST);
+//     glScissor(imin, jmin, imax-imin, jmax-jmin);
+    glBegin(GL_LINE_STRIP);
+    for (t=0; t<nvals; t++)
+    {
+        t1 = counter - t + nvals;
+        if (t1 >= NSTEPS) t1 = NSTEPS;
+        x = (double)imin + deltai*(double)t;
+        y = a*(double)input_signal[t1] + b;
+        glVertex2d(x, y);
+    }
+    glEnd();
+//     glDisable(GL_SCISSOR_TEST);
+    
+    glBegin(GL_LINE_LOOP);
+    glVertex2i(imin, jmin);
+    glVertex2i(imax, jmin);
+    glVertex2i(imax, jmax);
+    glVertex2i(imin, jmax);
+    glEnd();
+}
+
+void draw_wave_timeseries(double *values, int size, int fade, double fade_value)
+/* draw a profile of the wave, if option DRAW_WAVE_TIMESERIES is active */
+{
+    draw_exit_timeseries(values, size, fade, fade_value);
+    draw_entrance_timeseries(values, size, fade, fade_value);
+}
+
 void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], double *total_energy[NX], double *total_flux, short int *xy_in[NX], double scale, int time, int plot, int palette, int fade, double fade_value)
 /* same as draw_wave_highres, but with color scheme option */
 {
@@ -1017,7 +1260,7 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
     static double dtinverse = ((double)NX)/(COURANT*(XMAX-XMIN)), dx = (XMAX-XMIN)/((double)NX);
     double *values;
     
-    if (DRAW_WAVE_PROFILE) values = (double *)malloc(NX*NY*sizeof(double));
+    if ((DRAW_WAVE_PROFILE)||(DRAW_WAVE_TIMESERIES)) values = (double *)malloc(NX*NY*sizeof(double));
 
     glBegin(GL_QUADS);
     
@@ -1030,7 +1273,7 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
             {
                 value = wave_value(i, j, phi, psi, total_energy, total_flux, xy_in, scale, time, plot, palette, rgb);
                 
-                if (DRAW_WAVE_PROFILE) values[i*NY+j] = value;
+                if ((DRAW_WAVE_PROFILE)||(DRAW_WAVE_TIMESERIES)) values[i*NY+j] = value;
                 
                 if (fade) for (k=0; k<3; k++) rgb[k] *= fade_value;
                 glColor3f(rgb[0], rgb[1], rgb[2]);
@@ -1044,11 +1287,10 @@ void draw_wave_highres_palette(int size, double *phi[NX], double *psi[NX], doubl
 
     glEnd ();
     
-    if (DRAW_WAVE_PROFILE) 
-    {
-        draw_wave_profile(values, size, fade, fade_value);
-        free(values);
-    }
+    if (DRAW_WAVE_TIMESERIES) draw_wave_timeseries(values, size, fade, fade_value);
+    if (DRAW_WAVE_PROFILE) draw_wave_profile(values, size, fade, fade_value);
+    if ((DRAW_WAVE_PROFILE)||(DRAW_WAVE_TIMESERIES)) free(values);
+    
 }
 
 /* modified function for "flattened" wave tables */

wave_comparison.c

@@ -48,6 +48,11 @@
 #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 */
 
+#define VARIABLE_IOR 1      /* set to 1 for a variable index of refraction */
+#define IOR 10               /* 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 */
+#define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
+
 #define TIME_LAPSE 0     /* set to 1 to add a time-lapse movie at the end */
 #define TIME_LAPSE_FACTOR 4    /* factor of time-lapse movie */
 
@@ -80,8 +85,8 @@
 
 /* Choice of the billiard table */
 
-#define B_DOMAIN 999      /* choice of domain shape, see list in global_pdes.c */
-#define B_DOMAIN_B 999    /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN 63      /* choice of domain shape, see list in global_pdes.c */
+#define B_DOMAIN_B 63    /* choice of domain shape, see list in global_pdes.c */
 
 #define CIRCLE_PATTERN 13      /* pattern of circles, see list in global_pdes.c */
 #define CIRCLE_PATTERN_B 13     /* pattern of circles, see list in global_pdes.c */
@@ -97,11 +102,9 @@
 #define HEX_NONUNIF_COMPRESSSION 0.15 /* compression factor for HEX_NONUNIF pattern */
 #define HEX_NONUNIF_COMPRESSSION_B -0.15 /* compression factor for HEX_NONUNIF pattern */
 
-#define LAMBDA -1.1	    /* parameter controlling the dimensions of domain */
-#define MU 0.1   	    /* parameter controlling the dimensions of domain */
-#define MUB 0.1 	    /* parameter controlling the dimensions of domain */
-// #define MU 0.04665361 	    /* parameter controlling the dimensions of domain */
-// #define MUB 0.04665361	    /* parameter controlling the dimensions of domain */
+#define LAMBDA 1.5	    /* parameter controlling the dimensions of domain */
+#define MU 0.7             /* parameter controlling the dimensions of domain */
+#define MUB 0.2 	    /* parameter controlling the dimensions of domain */
 #define NPOLY 3             /* number of sides of polygon */
 #define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define APOLY_B 2.0         /* angle by which to turn polygon, in units of Pi/2 */ 
@@ -133,30 +136,25 @@
 #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 OMEGA 0.004        /* frequency of periodic excitation */
+#define OMEGA 0.024        /* frequency of periodic excitation */
 #define AMPLITUDE 1.0      /* amplitude of periodic excitation */ 
-#define COURANT 0.045       /* Courant number */
-#define COURANTB 0.045     /* Courant number in medium B */
-// #define COURANTB 0.005    /* Courant number in medium B */
-// #define COURANTB 0.008    /* Courant number in medium B */
-#define GAMMA 0.0      /* damping factor in wave equation */
-// #define GAMMA 1.0e-8      /* damping factor in wave equation */
-#define GAMMAB 0.0          /* damping factor in wave equation */
-// #define GAMMAB 1.0e-6           /* damping factor in wave equation */
-// #define GAMMAB 2.0e-4           /* damping factor in wave equation */
-// #define GAMMAB 2.5e-4           /* damping factor in wave equation */
+#define DAMPING 0.0        /* damping of periodic excitation */
+#define COURANT 0.1        /* Courant number */
+#define COURANTB 0.063      /* Courant number in medium B */
+#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 */
-#define GAMMA_TOPBOT 1.0e-6      /* damping factor on boundary */
-#define KAPPA 0.0       /* "elasticity" term enforcing oscillations */
-#define KAPPA_SIDES 5.0e-4       /* "elasticity" term on absorbing boundary */
-#define KAPPA_TOPBOT 0.0       /* "elasticity" term on absorbing boundary */
+#define GAMMA_TOPBOT 1.0e-7     /* damping factor on boundary */
+#define KAPPA 0.0           /* "elasticity" term enforcing oscillations */
+#define KAPPA_SIDES 5.0e-4  /* "elasticity" term on absorbing boundary */
+#define KAPPA_TOPBOT 0.0    /* "elasticity" term on absorbing boundary */
 /* The Courant number is given by c*DT/DX, where DT is the time step and DX the lattice spacing */
 /* The physical damping coefficient is given by GAMMA/(DT)^2 */
 /* Increasing COURANT speeds up the simulation, but decreases accuracy */
 /* 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 2     /* period of oscillating source */
+#define OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
 #define ALTERNATE_OSCILLATING_SOURCE 1  /* set to 1 to alternate sign of oscillating source */
 #define NSOURCES 48         /* number of sources */
 
@@ -168,9 +166,9 @@
 
 // #define NSTEPS 500      /* number of frames of movie */
 #define NSTEPS 2600      /* number of frames of movie */
-#define NVID 25          /* number of iterations between images displayed on screen */
+#define NVID 7            /* number of iterations between images displayed on screen */
 #define NSEG 100         /* number of segments of boundary */
-#define INITIAL_TIME 0    /* time after which to start saving frames */
+#define INITIAL_TIME 100    /* time after which to start saving frames */
 #define COMPUTE_ENERGIES 0  /* set to 1 to compute and print energies */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
 
@@ -178,26 +176,25 @@
 #define PSLEEP 1         /* sleep time during pause */
 #define SLEEP1  1        /* initial sleeping time */
 #define SLEEP2  1   /* final sleeping time */
-#define MID_FRAMES 20   /* number of still frames between movies */
-#define END_FRAMES 100   /* number of still frames at end of movie */
+#define MID_FRAMES 50   /* number of still frames between movies */
+#define END_FRAMES 300   /* number of still frames at end of movie */
 #define FADE 1           /* set to 1 to fade at end of movie */
 
 /* Parameters of initial condition */
 
-#define INITIAL_AMP 0.002          /* amplitude of initial condition */
-// #define INITIAL_AMP 0.002          /* amplitude of initial condition */
-#define INITIAL_VARIANCE 0.0003    /* variance of initial condition */
-#define INITIAL_WAVELENGTH  0.015  /* wavelength of initial condition */
+#define INITIAL_AMP 1.0            /* amplitude of initial condition */
+#define INITIAL_VARIANCE 0.00001    /* variance of initial condition */
+#define INITIAL_WAVELENGTH  0.025  /* wavelength of initial condition */
 
 /* Plot type, see list in global_pdes.c  */
 
 #define PLOT 0
 
-#define PLOT_B 3
+#define PLOT_B 5
 
 /* Color schemes */
 
-#define COLOR_PALETTE 18     /* Color palette, see list in global_pdes.c  */
+#define COLOR_PALETTE 17     /* 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 */
@@ -207,14 +204,13 @@
 
 #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 SLOPE 0.75        /* 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 ATTENUATION 0.0  /* exponential attenuation coefficient of contrast with time */
-#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 */
-#define FLUX_SCALE 1.0e4    /* scaling factor for enegy flux represtnation */
+#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 */
+#define FLUX_SCALE 5.0e3    /* scaling factor for enegy flux represtnation */
 #define RESCALE_COLOR_IN_CENTER 0   /* set to 1 to decrease color intentiy in the center (for wave escaping ring) */
 
 #define COLORHUE 260     /* initial hue of water color for scheme C_LUM */
@@ -224,7 +220,7 @@
 #define HUEMEAN 220.0    /* mean value of hue for color scheme C_HUE */
 #define HUEAMP -220.0      /* amplitude of variation of hue for color scheme C_HUE */
 
-#define DRAW_COLOR_SCHEME 1     /* set to 1 to plot the color scheme */
+#define DRAW_COLOR_SCHEME 0     /* set to 1 to plot the color scheme */
 #define COLORBAR_RANGE 1.5    /* scale of color scheme bar */
 #define COLORBAR_RANGE_B 12.5    /* scale of color scheme bar for 2nd part */
 #define ROTATE_COLOR_SCHEME 0   /* set to 1 to draw color scheme horizontally */
@@ -251,20 +247,13 @@
 #define POT_MAZE 7
 #define POTENTIAL 0
 #define MAZE_WIDTH 0.02     /* half width of maze walls */
-#define VARIABLE_IOR 1      /* 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 */
-#define MANDEL_IOR_SCALE -0.05   /* parameter controlling dependence of IoR on Mandelbrot escape speed */
-// #define COURANT 0.04       /* Courant number */
-// #define COURANTB 0.0       /* Courant number in medium B */
-// #define INITIAL_AMP 0.5            /* amplitude of initial condition */
-// #define INITIAL_VARIANCE 0.0003    /* variance of initial condition */
-// #define INITIAL_WAVELENGTH  0.015  /* wavelength of initial condition */
-// #define TWOSPEEDS 0          /* set to 1 to replace hardcore boundary by medium with different speed */
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
 #define N_WAVE_PACKETS 15               /* number of wave packets */
 #define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 #define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
+#define MU_B 1.0           /* parameter controlling the dimensions of domain */
+#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 */
 /* end of constants only used by sub_wave and sub_maze */
 
 
@@ -283,14 +272,14 @@ FILE *monitor_sources;
 /*********************/
 
 void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX],
-                      short int *xy_in[NX])
+                      short int *xy_in[NX], double *tcc[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 {
     int i, j, iplus, iminus, jplus, jminus, jmid = NY/2;
     double delta, x, y, c, cc, gamma;
     static long time = 0;
-    static double tc[NX][NY], tcc[NX][NY], tgamma[NX][NY];
+    static double tc[NX][NY], tgamma[NX][NY];
     static short int first = 1;
     
     time++;
@@ -300,16 +289,17 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
     {
         for (i=0; i<NX; i++){
             for (j=0; j<NY; j++){
-                if (xy_in[i][j])
+                if (xy_in[i][j] != 0)
                 {
                     tc[i][j] = COURANT;
-                    tcc[i][j] = courant2;
-                    tgamma[i][j] = GAMMA;
+                    if (!VARIABLE_IOR) tcc[i][j] = courant2;
+                    if (xy_in[i][j] == 1) tgamma[i][j] = GAMMA;
+                    else tgamma[i][j] = GAMMAB;
                 }
                 else if (TWOSPEEDS)
                 {
                     tc[i][j] = COURANTB;
-                    tcc[i][j] = courantb2;
+                    if (!VARIABLE_IOR) tcc[i][j] = courantb2;
                     tgamma[i][j] = GAMMAB;
                 }
             }
@@ -683,13 +673,13 @@ void evolve_wave_half(double *phi_in[NX], double *psi_in[NX], double *phi_out[NX
 }
 
 
-void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX])
+void evolve_wave(double *phi[NX], double *psi[NX], double *tmp[NX], short int *xy_in[NX], double *tcc_table[NX])
 /* time step of field evolution */
 /* phi is value of field at time t, psi at time t-1 */
 {
-    evolve_wave_half(phi, psi, tmp, xy_in);
-    evolve_wave_half(tmp, phi, psi, xy_in);
-    evolve_wave_half(psi, tmp, phi, xy_in);
+    evolve_wave_half(phi, psi, tmp, xy_in, tcc_table);
+    evolve_wave_half(tmp, phi, psi, xy_in, tcc_table);
+    evolve_wave_half(psi, tmp, phi, xy_in, tcc_table);
 }
 
 
@@ -712,8 +702,8 @@ void draw_color_bar_palette(int plot, double range, int palette, int fade, doubl
 
 void animation()
 {
-    double time, scale, energies[6], top_energy, bottom_energy, omega, angle, fade_value;
-    double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX];
+    double time, scale, energies[6], top_energy, bottom_energy, omega, angle, fade_value, sign = 1.0, ior_angle = 0.0;
+    double *phi[NX], *psi[NX], *tmp[NX], *total_energy[NX], *tcc_table[NX];
     short int *xy_in[NX];
     int i, j, s, counter = 0, k, first_source = 1, fade, resol = HIGHRES + 1;
     t_wave_source wave_source[NSOURCES];
@@ -728,6 +718,7 @@ void animation()
         tmp[i] = (double *)malloc(NY*sizeof(double));
         total_energy[i] = (double *)malloc(NY*sizeof(double));
         xy_in[i] = (short int *)malloc(NY*sizeof(short int));
+        tcc_table[i] = (double *)malloc(NX*sizeof(double));
     }
     
     /* initialise positions and radii of circles */
@@ -764,6 +755,7 @@ void animation()
 //     add_drop_to_wave(1.0, -0.7, 0.0, phi, psi);
 //     add_drop_to_wave(1.0, 0.0, -0.7, phi, psi);
 
+    if (VARIABLE_IOR) init_ior_2d_comp(xy_in, tcc_table, ior_angle);
     
     /* initialize energies */
     if (COMPUTE_ENERGIES) 
@@ -781,7 +773,7 @@ void animation()
     draw_wave_comp(phi, psi, xy_in, 1.0, 0, PLOT);
 
     printf("drawing billiard\n");
-    draw_billiard_comp();
+    draw_billiard_comp(0, 1.0);
 
     if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, fade, fade_value);
 
@@ -808,62 +800,59 @@ void animation()
         
        for (j=0; j<NVID; j++) 
         {
-            evolve_wave(phi, psi, tmp, xy_in);
+            evolve_wave(phi, psi, tmp, xy_in, tcc_table);
 //             if (i % 10 == 9) oscillate_linear_wave(0.2*scale, 0.15*(double)(i*NVID + j), -1.5, YMIN, -1.5, YMAX, phi, psi);
         }
         
         /* add oscillating waves */
         if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
         {
-//             if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
-//             add_circular_wave(sign, 0.0, 0.0, phi, psi, xy_in);
-
-//             p = phased_array_schedule(i);
+            if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
+            add_circular_wave(sign, -1.75, 0.16, phi, psi, xy_in);
+            add_circular_wave(sign, -1.75, -0.16, phi, psi, xy_in);
             
-//             phase_shift = 0.02 + 0.06*(double)i/(double)NSTEPS;
-            
-            if (first_source) for (k=0; k<NSOURCES; k++)
-            {
-                omega = DPI/(double)NSOURCES;
-                wave_source[k].xc = 0.05*cos(((double)k + 0.5)*omega);;
-                wave_source[k].yc = 0.05*sin(((double)k + 0.5)*omega);
-                if (wave_source[k].yc < 0.0)
-                {
-                    wave_source[k].xc *= 0.5;
-                    wave_source[k].yc *= 0.5;
-                }
-//                 wave_source[k].phase = 0.99 - 1.4*sin(0.35*(1.0 + wave_source[k].xc/0.1));
-//                 wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.1));
-//                 wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05));
-                wave_source[k].phase = 0.99 - 1.4*sin(0.35*(1.0 + wave_source[k].xc/0.05));
-                wave_source[k].amp = 1.0;
-//                 if (wave_source[k].phase) wave_source[k].sign = 1;
-//                 else wave_source[k].sign = -1;
-                wave_source[k].sign = 1;
-                
-                first_source = 0;
-            }
-            
-            fprintf(monitor_sources, "Frame %i\n\n", i);
-            
-            for (k=0; k<NSOURCES; k++) /*if (wave_source[k].xc > 0.0) */
-            {
-                wave_source[k].phase += 0.06;
-//                 if (k==1) printf("x = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].phase);
-//                 fprintf(monitor_sources, "x = %.3lg, y = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].yc, wave_source[k].phase);
-                if (wave_source[k].phase > 1.0)
-                {
-                    add_circular_wave_comp((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in, (wave_source[k].yc > 0));
-                    fprintf(monitor_sources, "x = %.3lg, y = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].yc, wave_source[k].phase);
-                    printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
-                    fprintf(monitor_sources, "Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
-                    wave_source[k].phase -= 1.0;
-                    wave_source[k].sign *= -1;
-                }
-            }
+//             if (first_source) for (k=0; k<NSOURCES; k++)
+//             {
+//                 omega = DPI/(double)NSOURCES;
+//                 wave_source[k].xc = 0.05*cos(((double)k + 0.5)*omega);;
+//                 wave_source[k].yc = 0.05*sin(((double)k + 0.5)*omega);
+//                 if (wave_source[k].yc < 0.0)
+//                 {
+//                     wave_source[k].xc *= 0.5;
+//                     wave_source[k].yc *= 0.5;
+//                 }
+// //                 wave_source[k].phase = 0.99 - 1.4*sin(0.35*(1.0 + wave_source[k].xc/0.1));
+// //                 wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.1));
+// //                 wave_source[k].phase = 0.99 - 1.4*sin(0.7*(1.0 + wave_source[k].xc/0.05));
+//                 wave_source[k].phase = 0.99 - 1.4*sin(0.35*(1.0 + wave_source[k].xc/0.05));
+//                 wave_source[k].amp = 1.0;
+// //                 if (wave_source[k].phase) wave_source[k].sign = 1;
+// //                 else wave_source[k].sign = -1;
+//                 wave_source[k].sign = 1;
+//                 
+//                 first_source = 0;
+//             }
+//             
+//             fprintf(monitor_sources, "Frame %i\n\n", i);
+//             
+//             for (k=0; k<NSOURCES; k++) /*if (wave_source[k].xc > 0.0) */
+//             {
+//                 wave_source[k].phase += 0.06;
+// //                 if (k==1) printf("x = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].phase);
+// //                 fprintf(monitor_sources, "x = %.3lg, y = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].yc, wave_source[k].phase);
+//                 if (wave_source[k].phase > 1.0)
+//                 {
+//                     add_circular_wave_comp((double)wave_source[k].sign*wave_source[k].amp, wave_source[k].xc, wave_source[k].yc, phi, psi, xy_in, (wave_source[k].yc > 0));
+//                     fprintf(monitor_sources, "x = %.3lg, y = %.3lg, phase = %.3lg\n", wave_source[k].xc, wave_source[k].yc, wave_source[k].phase);
+//                     printf("Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
+//                     fprintf(monitor_sources, "Adding wave at (%.2lg, %.2lg)\n", wave_source[k].xc, wave_source[k].yc);
+//                     wave_source[k].phase -= 1.0;
+//                     wave_source[k].sign *= -1;
+//                 }
+//             }
         }
 
-        draw_billiard_comp();
+        draw_billiard_comp(0, 1.0);
         
 	if (COMPUTE_ENERGIES) 
         {
@@ -899,7 +888,7 @@ void animation()
                 {
 //                     draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B);
                     draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, i, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
-                    draw_billiard_comp();
+                    draw_billiard_comp(0, 1.0);
                     if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0);
                     glutSwapBuffers();
                     save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter);
@@ -930,7 +919,7 @@ void animation()
         {
 //             draw_wave_comp(phi, psi, xy_in, scale, i, PLOT);
             draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 0, 1.0);
-            draw_billiard_comp();
+            draw_billiard_comp(0, 1.0);
             if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 0, 1.0);   
             glutSwapBuffers();
         }
@@ -939,7 +928,7 @@ void animation()
         {
             fade_value = 1.0 - (double)i/(double)MID_FRAMES;
             draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT, COLOR_PALETTE, 1, fade_value);
-            draw_billiard_comp();
+            draw_billiard_comp(1, fade_value);
             if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT, COLORBAR_RANGE, COLOR_PALETTE, 1, fade_value); 
             if (!NO_EXTRA_BUFFER_SWAP) glutSwapBuffers();
             save_frame_counter(NSTEPS + i + 1);
@@ -948,7 +937,7 @@ void animation()
         {
 //             draw_wave_comp(phi, psi, xy_in, scale, i, PLOT_B);
             draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 0, 1.0);
-            draw_billiard_comp();
+            draw_billiard_comp(0, 1.0);
             if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 0, 1.0); 
             glutSwapBuffers();
             
@@ -957,7 +946,7 @@ void animation()
             {
                 fade_value = 1.0 - (double)i/(double)END_FRAMES;
                 draw_wave_comp_highres_palette(resol, phi, psi, total_energy, xy_in, scale, NSTEPS, PLOT_B, COLOR_PALETTE_B, 1, fade_value);
-                draw_billiard_comp();
+                draw_billiard_comp(1, fade_value);
                 if (DRAW_COLOR_SCHEME) draw_color_bar_palette(PLOT_B, COLORBAR_RANGE_B, COLOR_PALETTE_B, 1, fade_value); 
                 glutSwapBuffers();
                 save_frame_counter(NSTEPS + MID_FRAMES + 1 + counter + i);
@@ -975,6 +964,7 @@ void animation()
         free(tmp[i]);
         free(total_energy[i]);
         free(xy_in[i]);
+        free(tcc_table[i]);
     }
 
 }

wave_energy.c

@@ -224,6 +224,10 @@
 #define OSCIL_LEFT_YSHIFT 0.0   /* y-dependence of left oscillation (for non-horizontal waves) */
 #define DRAW_WAVE_PROFILE 0     /* set to 1 to draw a profile of the wave */
 #define OSCILLATING_SOURCE_PERIOD 20    /* period of oscillating source */
+#define MU_B 1.0           /* parameter controlling the dimensions of domain */
+#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 DRAW_WAVE_TIMESERIES 0  /* set to 1 to draw a time series of the wave */
 /* end of constants only used by sub_wave and sub_maze */
 
 #include "global_pdes.c"        /* constants and global variables */
@@ -847,7 +851,7 @@ void animation()
     draw_wave_energy(phi, psi, xy_in, 1.0, 0);
 
     printf("drawing billiard\n");
-    draw_billiard_half(B_DOMAIN, CIRCLE_PATTERN, 0);
+    draw_billiard_half(B_DOMAIN, CIRCLE_PATTERN, 0, 0, 1.0);
 
     glutSwapBuffers();
 
@@ -869,7 +873,7 @@ void animation()
 
         draw_wave_energy(phi, psi, xy_in, scale, i);
         
-        draw_billiard_half(B_DOMAIN, CIRCLE_PATTERN, 0);
+        draw_billiard_half(B_DOMAIN, CIRCLE_PATTERN, 0, 0, 1.0);
 
 
         

wave_sphere.c

@@ -48,12 +48,22 @@
 
 #define WINWIDTH 	1920  /* window width */
 #define WINHEIGHT 	1150  /* window height */
+// // // #define NX 2500          /* number of grid points on x axis */
+// // // #define NY 1250          /* number of grid points on y axis */
 #define NX 2560          /* number of grid points on x axis */
 #define NY 1280           /* number of grid points on y axis */
+// #define NX 1280          /* number of grid points on x axis */
+// #define NY 640           /* number of grid points on y axis */
+// #define NX 2048         /* number of grid points on x axis */
+// #define NY 1024         /* number of grid points on y axis */
+// #define NX 3064         /* number of grid points on x axis */
+// #define NY 1536         /* number of grid points on y axis */
+// #define NX 1024         /* number of grid points on x axis */
+// #define NY 512         /* number of grid points on y axis */
 
 #define DPOLE 20         /* safety distance to poles */
 #define SMOOTHPOLE 0.1     /* smoothing coefficient at poles */
-#define ZERO_MERIDIAN 0.0     /* choice of zero meridian (will be at left/right boundary of 2d plot) */
+#define ZERO_MERIDIAN 180.0     /* choice of zero meridian (will be at left/right boundary of 2d plot) */
 
 #define XMIN -2.0
 #define XMAX 2.0	/* x interval */
@@ -66,6 +76,7 @@
 #define JULIA_ROT -20.0       /* rotation of Julia set, in degrees */
 #define JULIA_RE 0.5    
 #define JULIA_IM 0.462    /* parameters for Julia sets */
+// #define JULIA_IM 0.45    /* parameters for Julia sets */
 
 /* Choice of the billiard table */
 
@@ -87,7 +98,8 @@
 #define RANDOM_POLY_ANGLE 1 /* set to 1 to randomize angle of polygons */
 
 #define LAMBDA 0.75	    /* parameter controlling the dimensions of domain */
-#define MU 0.1             /* parameter controlling the dimensions of domain */
+#define MU 0.1              /* parameter controlling the dimensions of domain */
+#define MU_B 1.0            /* parameter controlling the dimensions of domain */
 #define NPOLY 6             /* number of sides of polygon */
 #define APOLY 0.0           /* angle by which to turn polygon, in units of Pi/2 */ 
 #define MDEPTH 7            /* depth of computation of Menger gasket */
@@ -127,7 +139,7 @@
 #define COURANT 0.05       /* Courant number */
 #define COURANTB 0.002    /* Courant number in medium B */
 #define GAMMA 0.0          /* damping factor in wave equation */
-#define GAMMAB 1.0e-6         /* damping factor in wave equation */
+#define GAMMAB 1.0e-2         /* damping factor in wave equation */
 #define GAMMA_SIDES 1.0e-4      /* damping factor on boundary */
 #define GAMMA_TOPBOT 1.0e-7     /* damping factor on boundary */
 #define KAPPA 0.0           /* "elasticity" term enforcing oscillations */
@@ -140,8 +152,8 @@
 /* For similar wave forms, COURANT^2*GAMMA should be kept constant */
 
 #define ADD_OSCILLATING_SOURCE 0        /* set to 1 to add an oscillating wave source */
-#define OSCILLATING_SOURCE_PERIOD 30    /* period of oscillating source */
-#define ALTERNATE_OSCILLATING_SOURCE 1  /* set to 1 to alternate sign of oscillating source */
+#define OSCILLATING_SOURCE_PERIOD 25    /* period of oscillating source */
+#define ALTERNATE_OSCILLATING_SOURCE 0  /* set to 1 to alternate sign of oscillating source */
 
 #define ADD_WAVE_PACKET_SOURCES 0       /* set to 1 to add several sources emitting wave packets */
 #define WAVE_PACKET_SOURCE_TYPE 1       /* type of wave packet sources */
@@ -150,14 +162,17 @@
 
 /* Boundary conditions, see list in global_pdes.c  */
 
+// #define B_COND 1
 #define B_COND 2
 
 #define PRECOMPUTE_BC 0     /* set to 1 to compute neighbours for Laplacian in advance */
 
 /* Parameters for length and speed of simulation */
 
-#define NSTEPS 2500       /* number of frames of movie */
+#define NSTEPS 2400       /* number of frames of movie */
+// #define NSTEPS 500         /* number of frames of movie */
 #define NVID 4            /* number of iterations between images displayed on screen */
+// #define NVID 10            /* number of iterations between images displayed on screen */
 #define NSEG 1000          /* number of segments of boundary */
 #define INITIAL_TIME 0      /* time after which to start saving frames */
 #define BOUNDARY_WIDTH 2    /* width of billiard boundary */
@@ -173,9 +188,9 @@
 
 /* Parameters of initial condition */
 
-#define INITIAL_AMP 0.5            /* amplitude of initial condition */
-#define INITIAL_VARIANCE 0.001  /* variance of initial condition */
-#define INITIAL_WAVELENGTH  0.05  /* wavelength of initial condition */
+#define INITIAL_AMP 1.0            /* amplitude of initial condition */
+#define INITIAL_VARIANCE 0.00005    /* variance of initial condition */
+#define INITIAL_WAVELENGTH  0.002  /* wavelength of initial condition */
 
 /* Plot type, see list in global_pdes.c  */
 
@@ -188,8 +203,8 @@
 #define CHANGE_LUMINOSITY 1     /* set to 1 to let luminosity depend on energy flux intensity */
 #define FLUX_WINDOW 30          /* size of averaging window of flux intensity */
 #define AMPLITUDE_HIGH_RES 1    /* set to 1 to increase resolution of plot */
-#define SHADE_3D 0              /* set to 1 to change luminosity according to normal vector */
-#define SHADE_2D 1              /* set to 1 to change luminosity according to normal vector to plane */
+#define SHADE_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 to plane */
 #define SHADE_WAVE 1            /* set to 1 to have luminosity depend on wave height */
 #define NON_DIRICHLET_BC 0      /* set to 1 to draw only facets in domain, if field is not zero on boundary */
 #define FLOOR_ZCOORD 1          /* set to 1 to draw only facets with z not too negative */
@@ -208,7 +223,7 @@
 /* 3D representation */
 
 #define REPRESENTATION_3D 1     /* choice of 3D representation */ 
-#define PLOT_2D 1               /* switch to 2D representation, equirectangular projection */
+#define PLOT_2D 0               /* switch to 2D representation, equirectangular projection */
 #define PHISHIFT 0.0            /* shift of phi in 2D plot (in degrees) */
 #define FLOODING 1              /* set to 1 to draw waves above altitude (for Earth representations) */
 
@@ -216,10 +231,11 @@
 #define REP_PROJ_3D 1       /* projection on plane orthogonal to observer line of sight */
 
 #define ROTATE_VIEW 1       /* set to 1 to rotate position of observer */
-#define ROTATE_ANGLE 75.0   /* total angle of rotation during simulation */
+#define ROTATE_ANGLE 720.0   /* total angle of rotation during simulation */
+// #define ROTATE_ANGLE -50.0   /* total angle of rotation during simulation */
 #define ROTATE_VIEW_WHILE_FADE 1    /* set to 1 to keep rotating viewpoint during fade */
 
-#define VIEWPOINT_TRAJ 2    /* type of viewpoint trajectory */
+#define VIEWPOINT_TRAJ 1    /* type of viewpoint trajectory */
 #define MAX_LATITUDE 45.0   /* maximal latitude for viewpoint trajectory VP_ORBIT2 */
 
 /* Color schemes */
@@ -235,14 +251,14 @@
 #define COLOR_SCHEME 3   /* choice of color scheme, see list in global_pdes.c  */
 
 #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 SLOPE 0.5        /* sensitivity of color on wave amplitude */
 #define VSCALE_AMPLITUDE 1.5   /* additional scaling factor for color scheme P_3D_AMPLITUDE */
-#define VSCALE_ENERGY 4.0     /* additional scaling factor for color scheme P_3D_ENERGY */
+#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 */
 #define ATTENUATION 0.0    /* exponential attenuation coefficient of contrast with time */
-#define E_SCALE 150.0      /* scaling factor for energy representation */
-#define LOG_SCALE 0.75     /* scaling factor for energy log representation */
+#define E_SCALE 200.0      /* scaling factor for energy representation */
+#define LOG_SCALE 0.25      /* scaling factor for energy log representation */
 #define LOG_SHIFT 0.5      /* shift of colors on log scale */
 #define LOG_ENERGY_FLOOR -10.0    /* floor value for log of (total) energy */
 #define LOG_MEAN_ENERGY_SHIFT 1.0   /* additional shift for log of mean energy */
@@ -260,18 +276,20 @@
 #define NXMAZE 8      /* width of maze */
 #define NYMAZE 32      /* height of maze */
 #define MAZE_MAX_NGBH 5     /* max number of neighbours of maze cell */
-#define RAND_SHIFT 5        /* seed of random number generator */
+#define RAND_SHIFT 1        /* seed of random number generator */
 #define MAZE_XSHIFT 0.0     /* horizontal shift of maze */
 #define MAZE_WIDTH 0.02     /* half width of maze walls */
 
-#define DRAW_COLOR_SCHEME 0       /* set to 1 to plot the color scheme */
-#define COLORBAR_RANGE 2.0      /* scale of color scheme bar */
-#define COLORBAR_RANGE_B 3.0    /* scale of color scheme bar for 2nd part */
+#define DRAW_COLOR_SCHEME 1       /* set to 1 to plot the color scheme */
+#define COLORBAR_RANGE 6.0      /* scale of color scheme bar */
+#define COLORBAR_RANGE_B 6.0    /* scale of color scheme bar for 2nd part */
 #define ROTATE_COLOR_SCHEME 0     /* set to 1 to draw color scheme horizontally */
 #define CIRC_COLORBAR 0         /* set to 1 to draw circular color scheme */
 #define CIRC_COLORBAR_B 0       /* set to 1 to draw circular color scheme */
 
 #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 DRAW_WAVE_TIMESERIES 0  /* set to 1 to draw a time series of the wave */
 #define SAVE_TIME_SERIES 0      /* set to 1 to save wave time series at a point */
 
 #define ADD_POTENTIAL 0         /* set to 1 to add potential to z coordinate */
@@ -288,27 +306,30 @@
 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.816496581, 0.40824829, 0.40824829};      /* vector of "light" direction for P_3D_ANGLE color scheme */
-double observer[3] = {-5.0, 8.0, -7.0};    /* location of observer for REP_PROJ_3D representation */ 
+double light[3] = {-0.816496581, -0.40824829, 0.40824829};      /* vector of "light" direction for P_3D_ANGLE color scheme */
+double observer[3] = {-6.0, -3.0, 4.5};    /* location of observer for REP_PROJ_3D representation */ 
 int reset_view = 0;         /* switch to reset 3D view parameters (for option ROTATE_VIEW) */
 
 #define ADD_DEM 1               /* add DEM (digital elevation model) */
 #define ADD_NEGATIVE_DEM 0      /* add DEM with bathymetric data */
-#define RSCALE_DEM 0.1          /* scaling factor of radial component for DEM */
+#define RSCALE_DEM 0.075          /* scaling factor of radial component for DEM */
 #define SMOOTH_DEM 0            /* set to 1 to smoothen DEM (to make altitude less constant) */
-#define DEM_SMOOTH_STEPS 10     /* number of smoothening steps */
+#define DEM_SMOOTH_STEPS 5      /* number of smoothening steps */
 #define DEM_SMOOTH_HEIGHT 0.5   /* relative height below which to smoothen */
-#define DEM_MAXHEIGHT 8000      /* max height of DEM (estimated from Everest) */
-#define PLANET_SEALEVEL 2500.0      /* sea level for flooded planet */
+#define DEM_MAXHEIGHT 9000.0     /* max height of DEM (estimated from Everest/Olympus Mons) */
+#define PLANET_SEALEVEL 0.0      /* sea level for flooded planet */
+// #define PLANET_SEALEVEL 3850.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 RSCALE 0.01             /* scaling factor of radial component */
-#define RMAX 10.0               /* max value of radial component */
-#define Z_SCALING_FACTOR 0.85   /* overall scaling factor of z axis for REP_PROJ_3D representation */
-#define XY_SCALING_FACTOR 2.0   /* overall scaling factor for on-screen (x,y) coordinates after projection */
+#define RSCALE 0.025             /* scaling factor of radial component */
+#define RMAX 1.5               /* max value of radial component */
+#define RMIN 0.5               /* min value of radial component */
+#define Z_SCALING_FACTOR 0.8   /* overall scaling factor of z axis for REP_PROJ_3D representation */
+#define XY_SCALING_FACTOR 2.2   /* overall scaling factor for on-screen (x,y) coordinates after projection */
 #define ZMAX_FACTOR 1.0         /* max value of z coordinate for REP_PROJ_3D representation */
 #define XSHIFT_3D 0.0           /* overall x shift for REP_PROJ_3D representation */
 #define YSHIFT_3D 0.0           /* overall y shift for REP_PROJ_3D representation */
-#define COS_VISIBLE -0.5        /* limit on cosine of normal to shown facets */
+#define COS_VISIBLE -0.7        /* limit on cosine of normal to shown facets */
 
 #include "global_pdes.c"        /* constants and global variables */
 #include "global_3d.c"          /* constants and global variables */
@@ -712,8 +733,51 @@ void animation()
 //     for (j=1; j<NPOLY; j++)
 //         add_circular_wave_mod(1.0, lambda1*cos(((double)j+0.5)*angle), lambda1*sin(((double)j+0.5)*angle), phi, psi, xy_in);
         
+    /* Kilauea, Hawaii, USA */
+    init_circular_wave_sphere(155.2867*PI/180.0, 19.42109*PI/180.0, phi, psi, xy_in, wsphere);
+    
+    /* Mount Etna, Sicily, Italy */
+    add_circular_wave_sphere(1.0, -14.9950*PI/180.0, 37.7550*PI/180.0, phi, psi, xy_in, wsphere);
+    
+    /* Eyjafjallajökull, Iceland */
+    add_circular_wave_sphere(1.0, 19.61333*PI/180.0, 63.62000*PI/180.0, phi, psi, xy_in, wsphere);
+
+    /* Mount Vesuvius, Italy */
+    add_circular_wave_sphere(1.0, -14.433*PI/180.0, 40.817*PI/180.0, phi, psi, xy_in, wsphere);
+
+    /* Sakurajima, Japan */
+    add_circular_wave_sphere(1.0, -130.65806*PI/180.0, 31.58056*PI/180.0, phi, psi, xy_in, wsphere);
+
+    /* Mount Merapi, Indonesia */
+    add_circular_wave_sphere(1.0, -110.44611*PI/180.0, -7.54139*PI/180.0, phi, psi, xy_in, wsphere);
+    
+    /* Ulawun, Papua New Guinea */
+    add_circular_wave_sphere(1.0, -151.33333*PI/180.0, -5.05000*PI/180.0, phi, psi, xy_in, wsphere);
+    
+    /* Santa María, Guatemala */
+    add_circular_wave_sphere(1.0, 91.55167*PI/180.0, 14.75556*PI/180.0, phi, psi, xy_in, wsphere);
+    
+    /*  Mount Erebus, Antarctica */
+    add_circular_wave_sphere(1.0, -167.15333*PI/180.0, -77.52972*PI/180.0, phi, psi, xy_in, wsphere);
+    
+    /*  Piton de la Fournaise, Reunion island, Indian Ocean */
+    add_circular_wave_sphere(0.75, -55.70889*PI/180.0, -21.24250*PI/180.0, phi, psi, xy_in, wsphere);
+
+    /*  Soufrière Hills, Montserrat */
+    add_circular_wave_sphere(0.75, 62.1773*PI/180.0, 16.7103*PI/180.0, phi, psi, xy_in, wsphere);
+    
+    /*  Axial Seamount, Pacific Ocean */
+    add_circular_wave_sphere(0.5, 130.00*PI/180.0, 45.95*PI/180.0, phi, psi, xy_in, wsphere);
+    
+    /*  Havre Seamount, New Zealand */
+    add_circular_wave_sphere(0.5, 179.968611*PI/180.0, -31.120278*PI/180.0, phi, psi, xy_in, wsphere);
+
+    /*  Boomerang Seamount, Indian Ocean */
+    add_circular_wave_sphere(0.5, -77.825*PI/180.0, -37.721*PI/180.0, phi, psi, xy_in, wsphere);
+
 //     init_wave_flat_sphere(phi, psi, xy_in, wsphere);
-    init_circular_wave_sphere(123.3*PI/180.0, -48.8*PI/180.0, phi, psi, xy_in, wsphere);
+//     for (j=0; j<6; j++)
+//         add_circular_wave_sphere(1.0, (double)j*PI/3.0, 0.925*PID, phi, psi, xy_in, wsphere);
     
 //     add_circular_wave_sphere(1.0, 1.0 - PI/9.0 + DPI/3.0, 0.15, phi, psi, xy_in, wsphere);
 //     add_circular_wave_sphere(1.0, 1.0 - PI/9.0 + 2.0*DPI/3.0, 0.15, phi, psi, xy_in, wsphere);
@@ -801,7 +865,8 @@ void animation()
         if ((ADD_OSCILLATING_SOURCE)&&(i%OSCILLATING_SOURCE_PERIOD == 1))
         {
             if (ALTERNATE_OSCILLATING_SOURCE) sign = -sign;
-            add_circular_wave_mod(sign, -0.5, 0.0, phi, psi, xy_in);
+            for (j=0; j<6; j++)
+                add_circular_wave_sphere(sign, (double)j*PI/3.0, 0.925*PID, phi, psi, xy_in, wsphere);
         }
         if (PRINT_SPEED) print_speed_3d(speed, 0, 1.0);